Project
Loading...
Searching...
No Matches
genEvents.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#include <iostream>
16#include <fstream>
17#include <cstring>
18#include <sys/stat.h>
19
20#include "Rtypes.h"
21#include "TRandom.h"
22#include "TH1F.h"
23#include "TFile.h"
24#include "TCanvas.h"
25#include "TPad.h"
26
27#include <iostream>
28#include <iomanip>
29#include <limits>
30
31#include "genEvents.h"
32#include "GPUTPCClusterData.h"
33#include "GPUTPCMCInfo.h"
35#include "GPUParam.inc"
37#include "GPUTPCGMPropagator.h"
38#include "GPUTPCGMMerger.h"
39#include "GPUChainTracking.h"
40#include "GPUConstantMem.h"
41#include "GPUSettings.h"
42
43using namespace o2::gpu;
44using namespace std;
45namespace o2::gpu
46{
47extern GPUSettingsStandalone configStandalone;
48}
49
50int32_t genEvents::GetSector(double GlobalPhi)
51{
52 double phi = GlobalPhi;
53 // std::cout<<" GetSector: phi = "<<phi<<std::endl;
54
55 if (phi >= mTwoPi) {
56 phi -= mTwoPi;
57 }
58 if (phi < 0) {
59 phi += mTwoPi;
60 }
61 return (int32_t)(phi / mSectorDAngle);
62}
63
64int32_t genEvents::GetDSector(double LocalPhi) { return GetSector(LocalPhi + mSectorAngleOffset); }
65
66double genEvents::GetSectorAngle(int32_t iSector) { return mSectorAngleOffset + iSector * mSectorDAngle; }
67
68int32_t genEvents::RecalculateSector(GPUTPCGMPhysicalTrackModel& t, int32_t& iSector)
69{
70 double phi = atan2(t.GetY(), t.GetX());
71 // std::cout<<" recalculate: phi = "<<phi<<std::endl;
72 int32_t dSector = GetDSector(phi);
73
74 if (dSector == 0) {
75 return 0; // nothing to do
76 }
77 // std::cout<<" dSector = "<<dSector<<std::endl;
78 double dAlpha = dSector * mSectorDAngle;
79 // rotate track on angle dAlpha
80
81 t.Rotate(dAlpha);
82
83 iSector += dSector;
84 if (iSector >= 18) {
85 iSector -= 18;
86 }
87 return 1;
88}
89
90double genEvents::GetGaus(double sigma)
91{
92 double x = 0;
93 do {
94 x = gRandom->Gaus(0., sigma);
95 if (fabs(x) <= 3.5 * sigma) {
96 break;
97 }
98 } while (1);
99 return x;
100}
101
103{
104 const char* rows[3] = {"0-63", "128-159", "64-127"};
105 for (int32_t i = 0; i < 3; i++) {
106 for (int32_t j = 0; j < 2; j++) {
107 char name[1024], title[1024];
108
109 snprintf(name, 1024, "clError%s%d", (j == 0 ? "Y" : "Z"), i);
110
111 snprintf(title, 1024, "Cluster %s Error for rows %s", (j == 0 ? "Y" : "Z"), rows[i]);
112
113 mClusterError[i][j] = new TH1F(name, title, 1000, 0., .7);
114 mClusterError[i][j]->GetXaxis()->SetTitle("Cluster Error [cm]");
115 }
116 }
117}
118
120{
121 TFile* tout = new TFile("generator.root", "RECREATE");
122 TCanvas* c = new TCanvas("ClusterErrors", "Cluste rErrors", 0, 0, 700, 700. * 2. / 3.);
123 c->Divide(3, 2);
124 int32_t ipad = 1;
125 for (int32_t j = 0; j < 2; j++) {
126 for (int32_t i = 0; i < 3; i++) {
127 c->cd(ipad++);
128 int32_t k = i;
129 if (i == 1) {
130 k = 2;
131 }
132 if (i == 2) {
133 k = 1;
134 }
135 if (tout) {
136 mClusterError[k][j]->Write();
137 }
138 gPad->SetLogy();
139 mClusterError[k][j]->Draw();
140 // delete mClusterError[i][j];
141 // mClusterError[i][j]=0;
142 }
143 }
144 c->Print("plots/clusterErrors.pdf");
145 delete c;
146 if (tout) {
147 tout->Close();
148 delete tout;
149 }
150}
151
152int32_t genEvents::GenerateEvent(const GPUParam& param, const char* filename)
153{
154 mRec->ClearIOPointers();
155 static int32_t iEvent = -1;
156 iEvent++;
157 if (iEvent == 0) {
158 gRandom->SetSeed(configStandalone.seed);
159 }
160
161 int32_t nTracks = configStandalone.EG.numberOfTracks; // Number of MC tracks, must be at least as large as the largest fMCID assigned above
162 cout << "NTracks " << nTracks << endl;
163 std::vector<GPUTPCMCInfo> mcInfo(nTracks);
164 memset(mcInfo.data(), 0, nTracks * sizeof(mcInfo[0]));
165
166 // double Bz = param.ConstBz();
167 // std::cout<<"Bz[kG] = "<<param.BzkG()<<std::endl;
168
170 {
171 const GPUTPCGMMerger& merger = mRec->GetProcessors()->tpcMerger;
172 prop.SetPolynomialField(&merger.Param().polynomialField);
173 }
174
175 // Bz*=o2::gpu::gpu_common_constants::kCLight;
176
177 std::vector<GenCluster> vClusters;
178 int32_t clusterId = 0; // Here we count up the cluster ids we fill (must be unique).
179 // gRandom->SetSeed(0);
180 // uint32_t seed = gRandom->GetSeed();
181
182 for (int32_t itr = 0; itr < nTracks; itr++) {
183 // std::cout<<"Track "<<itr<<":"<<std::endl;
184 // gRandom->SetSeed(seed);
185
186 mcInfo[itr].pid = -100; //-100: Unknown / other, 0: Electron, 1, Muon, 2: Pion, 3: Kaon, 4: Proton
187 mcInfo[itr].charge = 1;
188 mcInfo[itr].prim = 1; // Primary particle
189 mcInfo[itr].primDaughters = 0; // Primary particle with daughters in the TPC
190 mcInfo[itr].x = 0; // Position of MC track at entry of TPC / first hit in the TPC
191 mcInfo[itr].y = 0;
192 mcInfo[itr].z = 0;
193 mcInfo[itr].pX = 0; // Momentum of MC track at that position
194 mcInfo[itr].pY = 0;
195 mcInfo[itr].pZ = 0;
196
198 double dphi = mTwoPi / nTracks;
199 double phi = mSectorAngleOffset + dphi * itr;
200 double eta = gRandom->Uniform(-1.5, 1.5);
201
202 double theta = 2 * std::atan(1. / exp(eta));
203 double lambda = theta - M_PI / 2;
204 // double theta = gRandom->Uniform(-60,60)*M_PI/180.;
205 double pt = .08 * std::pow(10, gRandom->Uniform(0, 2.2));
206
207 double q = 1.;
208 int32_t iSector = GetSector(phi);
209 phi = phi - GetSectorAngle(iSector);
210
211 // std::cout<<"phi = "<<phi<<std::endl;
212 double x0 = cosf(phi);
213 double y0 = sinf(phi);
214 double z0 = tanf(lambda);
215 t.Set(x0, y0, z0, pt * x0, pt * y0, pt * z0, q);
216
217 if (RecalculateSector(t, iSector) != 0) {
218 std::cout << "Initial sector wrong!!!" << std::endl;
219 // exit(0);
220 }
221
222 for (int32_t iRow = 0; iRow < GPUCA_ROW_COUNT; iRow++) {
223 // if( iRow>=50 ) break; //SG!!!
224 float xRow = GPUTPCGeometry::Row2X(iRow);
225 // transport to row
226 int32_t err = 0;
227 for (int32_t itry = 0; itry < 1; itry++) {
228 float B[3];
229 prop.GetBxByBz(GetSectorAngle(iSector), t.GetX(), t.GetY(), t.GetZ(), B);
230 float dLp = 0;
231 err = t.PropagateToXBxByBz(xRow, B[0], B[1], B[2], dLp);
232 if (err) {
233 std::cout << "Can not propagate to x = " << xRow << std::endl;
234 t.Print();
235 break;
236 }
237 if (fabsf(t.GetZ()) >= GPUTPCGeometry::TPCLength()) {
238 std::cout << "Can not propagate to x = " << xRow << ": Z outside the volume" << std::endl;
239 t.Print();
240 err = -1;
241 break;
242 }
243 // rotate track coordinate system to current sector
244 int32_t isNewSector = RecalculateSector(t, iSector);
245 if (!isNewSector) {
246 break;
247 } else {
248 std::cout << "track " << itr << ": new sector " << iSector << " at row " << iRow << std::endl;
249 }
250 }
251 if (err) {
252 break;
253 }
254 // std::cout<<" track "<<itr<<": Sector "<<iSector<<" row "<<iRow<<" params :"<<std::endl;
255 // t.Print();
256 // track at row iRow, sector iSector
257 if (iRow == 0) { // store MC track at first row
258 // std::cout<<std::setprecision( 20 );
259 // std::cout<<"track "<<itr<<": x "<<t.X()<<" y "<<t.Y()<<" z "<<t.Z()<<std::endl;
260 GPUTPCGMPhysicalTrackModel tg(t); // global coordinates
261 tg.Rotate(-GetSectorAngle(iSector));
262
263 mcInfo[itr].pid = 2; // pion
264 mcInfo[itr].charge = 3 * q;
265 mcInfo[itr].x = tg.GetX(); // Position of MC track at entry of TPC / first hit in the TPC
266 mcInfo[itr].y = tg.GetY();
267 mcInfo[itr].z = tg.GetZ();
268 mcInfo[itr].pX = tg.GetPx(); // Momentum of MC track at that position
269 mcInfo[itr].pY = tg.GetPy();
270 mcInfo[itr].pZ = tg.GetPz();
271 // std::cout<<" mc Z = "<<tg.GetZ()<<std::endl;
272 }
273
274 // create cluster
275 GenCluster c;
276 float sigmaY = 0.3;
277 float sigmaZ = 0.5;
278 const int32_t rowType = iRow < 64 ? 0 : iRow < 128 ? 2 : 1;
279 t.UpdateValues();
280 param.GetClusterErrors2(iSector, rowType, t.GetZ(), t.GetSinPhi(), t.GetDzDs(), -1.f, 0.f, 0.f, sigmaY, sigmaZ);
281 sigmaY = std::sqrt(sigmaY);
282 sigmaZ = std::sqrt(sigmaZ);
283 mClusterError[rowType][0]->Fill(sigmaY);
284 mClusterError[rowType][1]->Fill(sigmaZ);
285 // std::cout<<sigmaY<<" "<<sigmaY<<std::endl;
286 // if( sigmaY > 0.5 ) sigmaY = 0.5;
287 // if( sigmaZ > 0.5 ) sigmaZ = 0.5;
288 c.sector = (t.GetZ() >= 0.) ? iSector : iSector + 18;
289 c.row = iRow;
290 c.mcID = itr;
291 c.x = t.GetX();
292 c.y = t.GetY() + GetGaus(sigmaY);
293 c.z = t.GetZ() + GetGaus(sigmaZ);
294 c.id = clusterId++;
295 vClusters.push_back(c);
296 } // iRow
297 } // itr
298
299 std::vector<AliHLTTPCClusterMCLabel> labels;
300
301 std::unique_ptr<GPUTPCClusterData> clSectors[GPUChainTracking::NSECTORS];
302
303 for (int32_t iSector = 0; iSector < (int32_t)GPUChainTracking::NSECTORS; iSector++) // HLT Sector numbering, sectors go from 0 to 35, all spanning all rows from 0 to 158.
304 {
305 int32_t nNumberOfHits = 0;
306 for (uint32_t i = 0; i < vClusters.size(); i++) {
307 if (vClusters[i].sector == iSector) {
308 nNumberOfHits++;
309 }
310 }
311 // For every sector we first have to fill the number of hits in this sector to the file
312 mRec->mIOPtrs.nClusterData[iSector] = nNumberOfHits;
313
314 GPUTPCClusterData* clusters = new GPUTPCClusterData[nNumberOfHits];
315 clSectors[iSector].reset(clusters);
316 int32_t icl = 0;
317 for (uint32_t i = 0; i < vClusters.size(); i++) {
318 GenCluster& c = vClusters[i];
319 if (c.sector == iSector) {
320 clusters[icl].id = c.id;
321 clusters[icl].row = c.row; // We fill one hit per TPC row
322 clusters[icl].x = c.x;
323 clusters[icl].y = c.y;
324 clusters[icl].z = c.z;
325 clusters[icl].amp = 100; // Arbitrary amplitude
326 icl++;
327 AliHLTTPCClusterMCLabel clusterLabel;
328 for (int32_t j = 0; j < 3; j++) {
329 clusterLabel.fClusterID[j].fMCID = -1;
330 clusterLabel.fClusterID[j].fWeight = 0;
331 }
332 clusterLabel.fClusterID[0].fMCID = c.mcID;
333 clusterLabel.fClusterID[0].fWeight = 1;
334 labels.push_back(clusterLabel);
335 }
336 }
337 mRec->mIOPtrs.clusterData[iSector] = clusters;
338 }
339
340 // Create vector with cluster MC labels, clusters are counter from 0 to clusterId in the order they have been written above. No separation in sectors.
341
342 mRec->mIOPtrs.nMCLabelsTPC = labels.size();
343 mRec->mIOPtrs.mcLabelsTPC = labels.data();
344
345 mRec->mIOPtrs.nMCInfosTPC = mcInfo.size();
346 mRec->mIOPtrs.mcInfosTPC = mcInfo.data();
347 static const GPUTPCMCInfoCol mcColInfo = {0, (uint32_t)mcInfo.size()};
348 mRec->mIOPtrs.mcInfosTPCCol = &mcColInfo;
349 mRec->mIOPtrs.nMCInfosTPCCol = 1;
350
351 mRec->DumpData(filename);
352 labels.clear();
353 mcInfo.clear();
354 return (0);
355}
356
357void genEvents::RunEventGenerator(GPUChainTracking* rec, const std::string& dir)
358{
359 std::unique_ptr<genEvents> gen(new genEvents(rec));
360 mkdir(dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
361 rec->DumpSettings(dir.c_str());
362
363 gen->InitEventGenerator();
364
365 for (int32_t i = 0; i < (configStandalone.nEvents == -1 ? 10 : configStandalone.nEvents); i++) {
366 GPUInfo("Generating event %d/%d", i, configStandalone.nEvents == -1 ? 10 : configStandalone.nEvents);
367 gen->GenerateEvent(rec->GetParam(), (dir + GPUCA_EVDUMP_FILE "." + std::to_string(i) + ".dump").c_str());
368 }
369 gen->FinishEventGenerator();
370}
default_random_engine gen(dev())
uint64_t exp(uint64_t base, uint8_t exp) noexcept
int32_t i
#define GPUCA_EVDUMP_FILE
Definition GPUDef.h:40
#define GPUCA_ROW_COUNT
uint32_t j
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
Definition B.h:16
static constexpr int32_t NSECTORS
Definition GPUChain.h:58
const GPUParam & GetParam() const
void DumpSettings(const char *dir="")
float const GPUParam float float float float float int32_t int16_t int8_t float float float charge
void InitEventGenerator()
Definition genEvents.h:30
int32_t GenerateEvent(const GPUParam &sectorParam, const char *filename)
Definition genEvents.h:31
void FinishEventGenerator()
Definition genEvents.h:32
static void RunEventGenerator(GPUChainTracking *rec, const std::string &dir)
Definition genEvents.h:34
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint const GLchar * name
Definition glcorearb.h:781
GLuint GLfloat x0
Definition glcorearb.h:5034
GLenum GLfloat param
Definition glcorearb.h:271
GLuint GLfloat GLfloat y0
Definition glcorearb.h:5034
std::ostream & cout()
GPUSettingsStandalone configStandalone
Definition genEvents.cxx:47
Defining DataPointCompositeObject explicitly as copiable.
std::string to_string(gsl::span< T, Size > span)
Definition common.h:52
std::string filename()
GPUReconstruction * rec
AliHLTTPCClusterMCWeight fClusterID[3]
const int nEvents
Definition test_Fifo.cxx:27
std::vector< Cluster > clusters
std::vector< ReadoutWindowData > rows