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