]>
Commit | Line | Data |
---|---|---|
a7a1dd76 | 1 | /// \file AliTPCCmpNG.C |
3a4edebe | 2 | /// |
a7a1dd76 | 3 | /// version: 1.0 |
4 | /// description: | |
5 | /// define a class TPCGenTrack | |
6 | /// save TPC related properties of tracks into a single tree | |
3a4edebe | 7 | /// |
a7a1dd76 | 8 | /// input: |
9 | /// Int_t nEvents ... nr of events to process | |
10 | /// Int_t firstEventNr ... first event number (starts from 0) | |
11 | /// char* fnRecTracks .. name of file with reconstructed tracks | |
12 | /// char* fnHits ... name of file with hits and Kine Tree | |
13 | /// char* fnDigits ... name of file with digits | |
14 | /// char* fnTracks .. output file name, default genTracks.root | |
3a4edebe | 15 | /// |
a7a1dd76 | 16 | /// How to use: |
17 | /// Typical usage: | |
18 | /// .L AliTPCCmpNG.C+ | |
19 | /// TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root") | |
20 | /// t->Exec(); | |
21 | /// .q | |
22 | /// aliroot | |
23 | /// .L AliTPCCmpNG.C+ | |
24 | /// TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root"); | |
25 | /// t2->Exec(); | |
3a4edebe | 26 | /// |
a7a1dd76 | 27 | /// Details: |
3a4edebe | 28 | /// |
a7a1dd76 | 29 | /// Step 1 - summurize information from simulation |
3a4edebe | 30 | /// |
a7a1dd76 | 31 | /// Compile macro with ACLIC: |
32 | /// .L AliTPCCmpNG.C+ | |
33 | /// create an object TPCFindGenTracks, which processes information | |
34 | /// from simulations. As input it needs: | |
35 | /// object gAlice: to get magnetic field | |
36 | /// TreeK: to get parameters of generated particles | |
37 | /// TreeTR: to get track parameters at the TPC entry | |
38 | /// TreeD: to get number of digits and digits pattern | |
39 | /// for a given track in TPC | |
40 | /// These input objects can be in different files, gAlice, TreeK and | |
41 | /// TreeTR are in the file fnHits, TreeD in the file fnDigits (can be | |
42 | /// the same as fnHits. Output is written to the file fnRes | |
43 | /// ("genTracks.root" by default). Use can specify number of | |
44 | /// events to process and the first event number: | |
45 | /// TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root","genTracks.root",1,0) | |
46 | /// The filenames in the example on previous line are defaults, user can | |
47 | /// specify just the file name with gAlice object (and TreeTR and TreeK), | |
48 | /// so equivalent shorter initialization is: | |
49 | /// TPCFindGenTracks *t = new TPCFindGenTracks("galice.root") | |
50 | /// The task is done by calling Exec() method: | |
51 | /// t->Exec(); | |
52 | /// User can set different debug levels by invoking: | |
53 | /// t->SetDebug(debugLevel) | |
54 | /// Number of events to process and the first event number can be | |
55 | /// specified as parameters to Exec: | |
56 | /// t->Exec(nEvents, firstEvent) | |
57 | /// Then you have to quit root to get rid of problems with deleting gAlice | |
58 | /// object (it is not deleted, but read again in the following step): | |
3a4edebe | 59 | /// |
a7a1dd76 | 60 | /// Step 2 - compare reconstructed tracks with simulated |
3a4edebe | 61 | /// |
a7a1dd76 | 62 | /// Load (and compile) the macro: |
63 | /// .L AliTPCCmpNG.C+ | |
64 | /// Create object TPCCmpTr, which does the comparison. As input it requires | |
65 | /// name of the file with reconstructed TPC tracks. You can specify | |
66 | /// name of the file with genTrack tree (summarized info about simulation), | |
67 | /// file with gAlice object, output file name, number of events to process | |
68 | /// and first event number: | |
69 | /// TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root","galice.root",1,0); | |
70 | /// The interface is quite similar to the TPCFindGenTracks class. | |
71 | /// Then just invoke Exec() method: | |
72 | /// t2->Exec(); | |
3a4edebe | 73 | /// |
a7a1dd76 | 74 | /// Step 3 - study the results |
3a4edebe | 75 | /// |
a7a1dd76 | 76 | /// Load the outoput TTree and you can do Draw(), Scan() or other |
77 | /// usual things to do with TTree: | |
78 | /// TFile *f = new TFile("cmpTracks.root") | |
79 | /// TTree *t = (TTree*)f->Get("TPCcmpTracks") | |
80 | /// t->Draw("fVDist[3]","fReconstructed") | |
3a4edebe | 81 | /// |
a7a1dd76 | 82 | /// History: |
3a4edebe | 83 | /// |
a7a1dd76 | 84 | /// 24.09.02 - first version |
85 | /// 24.01.03 - v7, before change from TPC Special Hits to TrackReferences | |
86 | /// 26.01.03 - change from TPC Special Hits to TrackReferences | |
87 | /// (loop over TreeTR instead of TreeH) | |
88 | /// 28.01.03 - v8 last version before removing TPC special point | |
89 | /// 28.01.03 - remove TPC special point, loop over TreeH | |
90 | /// store TParticle and AliTrack | |
91 | /// 29.01.03 - v9 last version before moving the loop over rec. tracks | |
92 | /// into separate step | |
93 | /// 03.02.03 - rename to AliTPCCmpNG.C, remove the part with rec. tracks | |
94 | /// (will be addded in a macro AliTPCCmpTr.C | |
3a4edebe | 95 | /// |
a7a1dd76 | 96 | /// \author Jiri Chudoba |
97 | /// \date 24.09.2002 | |
e8558587 | 98 | |
99 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
100 | #include "iostream.h" | |
101 | #include <stdio.h> | |
102 | #include <string.h> | |
103 | #include "Rtypes.h" | |
104 | #include "TFile.h" | |
105 | #include "TTree.h" | |
106 | #include "TString.h" | |
107 | #include "TBenchmark.h" | |
108 | #include "TStopwatch.h" | |
109 | #include "TParticle.h" | |
110 | #include "AliRun.h" | |
111 | #include "AliStack.h" | |
112 | #include "AliTPCtrack.h" | |
113 | #include "AliSimDigits.h" | |
114 | #include "AliTPCParam.h" | |
115 | #include "TParticle.h" | |
116 | #include "AliTPC.h" | |
117 | #include "AliDetector.h" | |
118 | #include "AliTrackReference.h" | |
119 | #include "TSystem.h" | |
120 | #include "TTimer.h" | |
121 | // #include "AliConst.h" | |
122 | #endif | |
123 | ||
124 | #include "AliTPCTracking.C" | |
125 | ||
126 | // include ML.C: | |
127 | //////////////////////////////////////////////////////////////////////// | |
128 | // | |
129 | // name: ML.C (Macro Library - collection of functions used in diff macros | |
130 | // date: 08.05.2002 | |
131 | // last update: 08.05.2002 | |
132 | // author: Jiri Chudoba | |
133 | // version: 1.0 | |
134 | // description: | |
135 | // | |
136 | // History: | |
137 | // | |
138 | // 08.05.02 - first version | |
139 | // | |
140 | //////////////////////////////////////////////////////////////////////// | |
141 | ||
142 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
143 | #include "iostream.h" | |
144 | #include "Rtypes.h" | |
145 | #include "TSystem.h" | |
146 | #include "TTimer.h" | |
147 | #include "Getline.h" | |
148 | #include "TChain.h" | |
149 | #include "TString.h" | |
150 | #include "TFile.h" | |
151 | #include "AliRun.h" | |
152 | ||
153 | void WaitForReturn(); | |
154 | Bool_t ImportgAlice(TFile *file); | |
155 | TFile* OpenAliceFile(char *fn, Bool_t importgAlice = kFALSE, char *mode = "read"); | |
156 | Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain); | |
157 | ||
158 | #endif | |
159 | ||
160 | //////////////////////////////////////////////////////////////////////// | |
161 | void WaitForReturn() | |
162 | { | |
3a4edebe | 163 | /// wait until user press return; |
164 | ||
e8558587 | 165 | char *input; |
166 | Bool_t done = kFALSE; | |
167 | TTimer *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE); | |
168 | ||
169 | do { | |
170 | timer->TurnOn(); | |
171 | timer->Reset(); | |
172 | // Now let's read the input, we can use here any | |
173 | // stdio or iostream reading methods. like std::cin >> myinputl; | |
174 | input = Getline("Type <return> to continue: "); | |
175 | timer->TurnOff(); | |
176 | if (input) done = kTRUE; | |
177 | } while (!done); | |
178 | } | |
179 | ||
180 | //////////////////////////////////////////////////////////////////////// | |
181 | ||
182 | Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain) { | |
3a4edebe | 183 | /// chain all files fn in the subdirectories which match subDirName |
184 | /// return number of chained files | |
e8558587 | 185 | |
186 | // open baseDir, loop over subdirs | |
187 | ||
188 | void *dirp = gSystem->OpenDirectory(baseDir); | |
189 | if (!dirp) { | |
190 | cerr<<"Could not open base directory "<<baseDir.Data()<<endl; | |
191 | return 0; | |
192 | } | |
193 | const char *afile; | |
194 | Long_t id, size, flag, modTime; | |
195 | Int_t rc = 0; | |
196 | char relName[100]; | |
197 | while((afile = gSystem->GetDirEntry(dirp))) { | |
198 | // printf("file: %s\n",afile); | |
199 | if (strstr(afile,subDirNameMask.Data())) { | |
200 | sprintf(relName,"%s/%s/%s",baseDir.Data(),afile,fn.Data()); | |
201 | // cerr<<"relName = "<<relName<<endl; | |
202 | if(!gSystem->GetPathInfo(relName, &id, &size, &flag, &modTime)) { | |
203 | rc += chain.Add(relName); | |
204 | } | |
205 | } | |
206 | } | |
207 | gSystem->FreeDirectory(dirp); | |
208 | // cerr<<"rc = "<<rc<<endl; | |
209 | return rc; | |
210 | } | |
211 | //////////////////////////////////////////////////////////////////////// | |
212 | Bool_t ImportgAlice(TFile *file) { | |
3a4edebe | 213 | /// read in gAlice object from the file |
214 | ||
e8558587 | 215 | gAlice = (AliRun*)file->Get("gAlice"); |
216 | if (!gAlice) return kFALSE; | |
217 | return kTRUE; | |
218 | } | |
219 | //////////////////////////////////////////////////////////////////////// | |
220 | TFile* OpenAliceFile(char *fn, Bool_t importgAlice, char *mode) { | |
221 | TFile *file = TFile::Open(fn,mode); | |
222 | if (!file->IsOpen()) { | |
223 | cerr<<"OpenAliceFile: cannot open file "<<fn<<" in mode " | |
224 | <<mode<<endl; | |
225 | return 0; | |
226 | } | |
227 | if (!importgAlice) return file; | |
228 | if (ImportgAlice(file)) return file; | |
229 | return 0; | |
230 | } | |
231 | //////////////////////////////////////////////////////////////////////// | |
232 | ||
233 | ||
234 | //////////////////////////////////////////////////////////////////////// | |
235 | // | |
236 | // Start of implementation of the class TPCGenInfo | |
237 | // | |
238 | //////////////////////////////////////////////////////////////////////// | |
239 | class TPCGenInfo: public TObject { | |
240 | ||
241 | public: | |
242 | TPCGenInfo(); | |
a7a1dd76 | 243 | AliTrackReference *fTrackRef; ///< track reference saved in the output tree |
244 | TParticle *fParticle; ///< generated particle | |
245 | Int_t fLabel; ///< track label | |
e8558587 | 246 | |
a7a1dd76 | 247 | Int_t fRowsWithDigitsInn; ///< number of rows with digits in the inner sectors |
248 | Int_t fRowsWithDigits; ///< number of rows with digits in the outer sectors | |
249 | Int_t fRowsTrackLength; ///< last - first row with digit | |
250 | Int_t fDigitsInSeed; ///< digits in the default seed rows | |
e8558587 | 251 | |
252 | ClassDef(TPCGenInfo,1) // container for | |
253 | }; | |
254 | ClassImp(TPCGenInfo) | |
255 | //////////////////////////////////////////////////////////////////////// | |
256 | TPCGenInfo::TPCGenInfo() | |
257 | { | |
258 | fTrackRef = 0; | |
259 | fParticle = 0; | |
260 | } | |
261 | //////////////////////////////////////////////////////////////////////// | |
262 | ||
263 | //////////////////////////////////////////////////////////////////////// | |
264 | // | |
265 | // End of implementation of the class TPCGenInfo | |
266 | // | |
267 | //////////////////////////////////////////////////////////////////////// | |
268 | ||
269 | ||
270 | ||
271 | //////////////////////////////////////////////////////////////////////// | |
272 | // | |
273 | // Start of implementation of the class digitRow | |
274 | // | |
275 | //////////////////////////////////////////////////////////////////////// | |
276 | const Int_t kgRowBytes = 32; | |
277 | ||
278 | class digitRow: public TObject { | |
279 | ||
280 | public: | |
281 | digitRow(); | |
282 | // digitRow(){;} | |
283 | virtual ~digitRow(){;} | |
284 | void SetRow(Int_t row); | |
285 | Bool_t TestRow(Int_t row); | |
286 | digitRow & operator=(const digitRow &digOld); | |
287 | Int_t RowsOn(Int_t upto=8*kgRowBytes); | |
288 | Int_t Last(); | |
289 | Int_t First(); | |
290 | void Reset(); | |
291 | ||
292 | //private: | |
293 | UChar_t fDig[kgRowBytes]; | |
294 | ||
295 | ClassDef(digitRow,1) // container for digit pattern | |
296 | }; | |
297 | ClassImp(digitRow) | |
298 | //////////////////////////////////////////////////////////////////////// | |
299 | digitRow::digitRow() | |
300 | { | |
301 | Reset(); | |
302 | } | |
303 | //////////////////////////////////////////////////////////////////////// | |
304 | digitRow & digitRow::operator=(const digitRow &digOld) | |
305 | { | |
306 | for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i]; | |
307 | return (*this); | |
308 | } | |
309 | //////////////////////////////////////////////////////////////////////// | |
310 | void digitRow::SetRow(Int_t row) | |
311 | { | |
312 | if (row >= 8*kgRowBytes) { | |
313 | cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl; | |
314 | return; | |
315 | } | |
316 | Int_t iC = row/8; | |
317 | Int_t iB = row%8; | |
318 | SETBIT(fDig[iC],iB); | |
319 | } | |
320 | ||
321 | //////////////////////////////////////////////////////////////////////// | |
322 | Bool_t digitRow::TestRow(Int_t row) | |
323 | { | |
a7a1dd76 | 324 | /// return kTRUE if row is on |
325 | ||
e8558587 | 326 | Int_t iC = row/8; |
327 | Int_t iB = row%8; | |
328 | return TESTBIT(fDig[iC],iB); | |
329 | } | |
330 | //////////////////////////////////////////////////////////////////////// | |
331 | Int_t digitRow::RowsOn(Int_t upto) | |
332 | { | |
a7a1dd76 | 333 | /// returns number of rows with a digit |
334 | /// count only rows less equal row number upto | |
335 | ||
e8558587 | 336 | Int_t total = 0; |
337 | for (Int_t i = 0; i<kgRowBytes; i++) { | |
338 | for (Int_t j = 0; j < 8; j++) { | |
339 | if (i*8+j > upto) return total; | |
340 | if (TESTBIT(fDig[i],j)) total++; | |
341 | } | |
342 | } | |
343 | return total; | |
344 | } | |
345 | //////////////////////////////////////////////////////////////////////// | |
346 | void digitRow::Reset() | |
347 | { | |
a7a1dd76 | 348 | /// resets all rows to zero |
349 | ||
e8558587 | 350 | for (Int_t i = 0; i<kgRowBytes; i++) { |
351 | fDig[i] <<= 8; | |
352 | } | |
353 | } | |
354 | //////////////////////////////////////////////////////////////////////// | |
355 | Int_t digitRow::Last() | |
356 | { | |
a7a1dd76 | 357 | /// returns the last row number with a digit |
358 | /// returns -1 if now digits | |
359 | ||
e8558587 | 360 | for (Int_t i = kgRowBytes-1; i>=0; i--) { |
361 | for (Int_t j = 7; j >= 0; j--) { | |
362 | if TESTBIT(fDig[i],j) return i*8+j; | |
363 | } | |
364 | } | |
365 | return -1; | |
366 | } | |
367 | //////////////////////////////////////////////////////////////////////// | |
368 | Int_t digitRow::First() | |
369 | { | |
a7a1dd76 | 370 | /// returns the first row number with a digit |
371 | /// returns -1 if now digits | |
372 | ||
e8558587 | 373 | for (Int_t i = 0; i<kgRowBytes; i++) { |
374 | for (Int_t j = 0; j < 8; j++) { | |
375 | if (TESTBIT(fDig[i],j)) return i*8+j; | |
376 | } | |
377 | } | |
378 | return -1; | |
379 | } | |
380 | ||
381 | //////////////////////////////////////////////////////////////////////// | |
382 | // | |
383 | // end of implementation of a class digitRow | |
384 | // | |
385 | //////////////////////////////////////////////////////////////////////// | |
386 | //////////////////////////////////////////////////////////////////////// | |
387 | // | |
388 | // Start of implementation of the class TPCFindGenTracks | |
389 | // | |
390 | //////////////////////////////////////////////////////////////////////// | |
391 | ||
392 | class TPCFindGenTracks { | |
393 | ||
394 | public: | |
395 | TPCFindGenTracks(); | |
396 | TPCFindGenTracks(char* fnHits, | |
397 | char* fnDigits ="tpc.digits.root", | |
398 | char* fnRes ="genTracks.root", | |
399 | Int_t nEvents=1, Int_t firstEvent=0); | |
400 | virtual ~TPCFindGenTracks(); | |
401 | void Reset(); | |
402 | Int_t Exec(); | |
403 | Int_t Exec(Int_t nEvents, Int_t firstEventNr); | |
404 | void CreateTreeGenTracks(); | |
405 | void CloseOutputFile(); | |
406 | Int_t TreeKLoop(); | |
407 | Int_t TreeTRLoop(); | |
408 | Int_t TreeDLoop(); | |
409 | void SetFirstEventNr(Int_t i) {fFirstEventNr = i;} | |
410 | void SetNEvents(Int_t i) {fNEvents = i;} | |
411 | void SetDebug(Int_t level) {fDebug = level;} | |
412 | ||
413 | Float_t TR2LocalX(AliTrackReference *trackRef, | |
414 | AliTPCParam *paramTPC); | |
415 | ||
416 | public: | |
a7a1dd76 | 417 | Int_t fDebug; //!< debug flag |
418 | Int_t fEventNr; //!< current event number | |
419 | Int_t fLabel; //!< track label | |
420 | Int_t fNEvents; //!< number of events to process | |
421 | Int_t fFirstEventNr; //!< first event to process | |
422 | Int_t fNParticles; //!< number of particles in TreeK | |
423 | TTree *fTreeGenTracks; //!< output tree with generated tracks | |
424 | char *fFnRes; //!< output file name with stored tracks | |
425 | char *fFnHits; //!< input file name with hits | |
426 | char *fFnDigits; //!< input file name with digits | |
427 | TFile *fFileGenTracks; //!< output file with stored fTreeGenTracks | |
428 | TFile *fFileHits; //!< input file with hits | |
429 | TFile *fFileTreeD; //!< input file with digits | |
430 | digitRow *fDigitRow; //!< pointer to the object saved in Branch | |
431 | digitRow *fContainerDigitRow; //!< big container for partial information | |
432 | AliTrackReference *fTrackRef; //!< track reference saved in the output tree | |
433 | AliTrackReference *fContainerTR;//!< big container for partial information | |
434 | Int_t *fIndexTR; //!< index of particle label in the fContainerTR | |
435 | Int_t fLastIndexTR; //!< last used index in fIndexTR | |
436 | ||
437 | AliTPCParam* fParamTPC; //!< AliTPCParam | |
438 | ||
439 | Double_t fVPrim[3]; //!< primary vertex position | |
440 | Double_t fVDist[4]; //!< distance of the particle vertex from primary vertex | |
e8558587 | 441 | // the fVDist[3] contains size of the 3-vector |
a7a1dd76 | 442 | TParticle *fParticle; //!< generated particle |
e8558587 | 443 | |
a7a1dd76 | 444 | Int_t fRowsWithDigitsInn; //!< number of rows with digits in the inner sectors |
445 | Int_t fRowsWithDigits; //!< number of rows with digits in the outer sectors | |
446 | Int_t fRowsTrackLength; //!< last - first row with digit | |
447 | Int_t fDigitsInSeed; //!< digits in the default seed rows | |
e8558587 | 448 | |
449 | private: | |
450 | ||
451 | // some constants for the original non-pareller tracking (by Y.Belikov) | |
a7a1dd76 | 452 | static const Int_t seedRow11 = 158; ///< nRowUp - 1 |
453 | static const Int_t seedRow12 = 139; ///< nRowUp - 1 - (Int_t) 0.125*nRowUp | |
454 | static const Int_t seedRow21 = 149; ///< seedRow11 - shift | |
455 | static const Int_t seedRow22 = 130; ///< seedRow12 - shift | |
d8d3b5b8 | 456 | static const Double_t kRaddeg = 180./TMath::Pi(); |
e8558587 | 457 | |
a7a1dd76 | 458 | static const Int_t fgMaxIndexTR = 50000; ///< maximum number of tracks with a track ref |
459 | static const Int_t fgMaxParticles = 2000000; ///< maximum number of generated particles | |
460 | static const Double_t fgPtCut = .001; ///< do not store particles with generated pT less than this | |
e8558587 | 461 | static const Float_t fgTrackRefLocalXMax = 82.95; |
462 | static const Float_t fgTrackRefLocalXMaxDelta = 500.; | |
463 | ||
464 | ClassDef(TPCFindGenTracks,1) // class which creates and fills tree with TPCGenTrack objects | |
465 | }; | |
466 | ClassImp(TPCFindGenTracks) | |
467 | ||
468 | //////////////////////////////////////////////////////////////////////// | |
469 | TPCFindGenTracks::TPCFindGenTracks() | |
470 | { | |
471 | Reset(); | |
472 | } | |
473 | ||
474 | //////////////////////////////////////////////////////////////////////// | |
475 | TPCFindGenTracks::TPCFindGenTracks(char* fnHits, char* fnDigits, | |
476 | char* fnRes, | |
477 | Int_t nEvents, Int_t firstEvent) | |
478 | { | |
479 | Reset(); | |
480 | fFirstEventNr = firstEvent; | |
481 | fEventNr = firstEvent; | |
482 | fNEvents = nEvents; | |
483 | fFnRes = fnRes; | |
484 | fFnHits = fnHits; | |
485 | fFnDigits = fnDigits; | |
486 | } | |
487 | //////////////////////////////////////////////////////////////////////// | |
488 | void TPCFindGenTracks::Reset() | |
489 | { | |
490 | fDigitRow = 0; | |
491 | fEventNr = 0; | |
492 | fNEvents = 0; | |
493 | fTreeGenTracks = 0; | |
494 | fFnRes = "genTracks.root"; | |
495 | fFnHits = "rfio:galice.root"; | |
496 | fFnDigits = "rfio:its.tpc.trd.digits.root"; | |
497 | fFileGenTracks = 0; | |
498 | fFileHits =0; | |
499 | fFileTreeD =0; | |
500 | fContainerDigitRow = 0; | |
501 | fContainerTR = 0; | |
502 | fIndexTR = 0; | |
503 | fLastIndexTR = -1; | |
504 | fParticle = 0; | |
505 | fTrackRef = 0; | |
506 | fDebug = 0; | |
507 | fVPrim[0] = -1000.; | |
508 | fVPrim[1] = -1000.; | |
509 | fVPrim[2] = -1000.; | |
510 | fParamTPC = 0; | |
511 | ||
512 | } | |
513 | //////////////////////////////////////////////////////////////////////// | |
514 | TPCFindGenTracks::~TPCFindGenTracks() | |
515 | { | |
516 | ; | |
517 | } | |
518 | //////////////////////////////////////////////////////////////////////// | |
519 | Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr) | |
520 | { | |
521 | fNEvents = nEvents; | |
522 | fFirstEventNr = firstEventNr; | |
523 | return Exec(); | |
524 | } | |
525 | ||
526 | //////////////////////////////////////////////////////////////////////// | |
527 | Int_t TPCFindGenTracks::Exec() | |
528 | { | |
529 | TStopwatch timer; | |
530 | timer.Start(); | |
531 | ||
532 | fDigitRow = new digitRow(); | |
533 | CreateTreeGenTracks(); | |
534 | if (!fTreeGenTracks) return 1; | |
535 | fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here | |
536 | if (!fFileHits) return 1; | |
537 | fFileHits->cd(); | |
538 | SetFieldFactor(); | |
539 | ||
540 | fFileTreeD = TFile::Open(fFnDigits,"read"); | |
541 | if (!fFileTreeD->IsOpen()) { | |
542 | cerr<<"Cannot open file "<<fFnDigits<<endl; | |
543 | return 1; | |
544 | } | |
545 | ||
546 | fParamTPC = LoadTPCParam(fFileTreeD); | |
547 | if (!fParamTPC) { | |
548 | cerr<<"TPC parameters not found and could not be created"<<endl; | |
549 | return 1; | |
550 | } | |
551 | ||
552 | for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents; | |
553 | fEventNr++) { | |
554 | fNParticles = gAlice->GetEvent(fEventNr); | |
555 | fContainerDigitRow = new digitRow[fNParticles]; | |
556 | fContainerTR = new AliTrackReference[fgMaxIndexTR]; | |
557 | fIndexTR = new Int_t[fNParticles]; | |
558 | for (Int_t i = 0; i<fNParticles; i++) { | |
559 | fIndexTR[i] = -1; | |
560 | } | |
561 | ||
562 | cout<<"Start to process event "<<fEventNr<<endl; | |
563 | cout<<"\tfNParticles = "<<fNParticles<<endl; | |
564 | if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl; | |
565 | if (TreeDLoop()>0) return 1; | |
566 | if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl; | |
567 | if (TreeTRLoop()>0) return 1; | |
568 | if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl; | |
569 | if (TreeKLoop()>0) return 1; | |
570 | if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl; | |
571 | ||
572 | delete [] fContainerDigitRow; | |
573 | delete [] fContainerTR; | |
574 | delete [] fIndexTR; | |
575 | } | |
576 | ||
577 | CloseOutputFile(); | |
578 | ||
579 | cerr<<"Exec finished"<<endl; | |
580 | fFileHits->cd(); | |
581 | delete gAlice; | |
582 | fFileHits->Close(); | |
583 | delete fFileHits; | |
584 | ||
585 | timer.Stop(); | |
586 | timer.Print(); | |
587 | return 0; | |
588 | } | |
589 | //////////////////////////////////////////////////////////////////////// | |
590 | void TPCFindGenTracks::CreateTreeGenTracks() | |
591 | { | |
592 | fFileGenTracks = TFile::Open(fFnRes,"RECREATE"); | |
593 | if (!fFileGenTracks) { | |
594 | cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl; | |
595 | return; | |
596 | } | |
597 | fTreeGenTracks = new TTree("genTracksTree","genTracksTree"); | |
598 | TBranch *branchBits = fTreeGenTracks->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0); | |
599 | if (!branchBits) { | |
600 | cerr<<"Error in CreateTreeGenTracks: cannot create branch."<<endl; | |
601 | return; | |
602 | } | |
603 | fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I"); | |
604 | fTreeGenTracks->Branch("fLabel",&fLabel,"fLabel/I"); | |
605 | fTreeGenTracks->Branch("fRowsWithDigitsInn",&fRowsWithDigitsInn,"fRowsWithDigitsInn/I"); | |
606 | fTreeGenTracks->Branch("fRowsWithDigits",&fRowsWithDigits,"fRowsWithDigits/I"); | |
607 | fTreeGenTracks->Branch("fRowsTrackLength",&fRowsTrackLength,"fRowsTrackLength/I"); | |
608 | fTreeGenTracks->Branch("fDigitsInSeed",&fDigitsInSeed,"fDigitsInSeed/I"); | |
609 | ||
610 | fTreeGenTracks->Branch("Particle","TParticle",&fParticle); | |
611 | fTreeGenTracks->Branch("fVDist",&fVDist,"fVDist[4]/D"); | |
612 | fTreeGenTracks->Branch("TR","AliTrackReference",&fTrackRef); | |
613 | ||
614 | fTreeGenTracks->AutoSave(); | |
615 | } | |
616 | //////////////////////////////////////////////////////////////////////// | |
617 | void TPCFindGenTracks::CloseOutputFile() | |
618 | { | |
619 | if (!fFileGenTracks) { | |
620 | cerr<<"File "<<fFnRes<<" not found as an open file."<<endl; | |
621 | return; | |
622 | } | |
623 | fFileGenTracks->cd(); | |
624 | fTreeGenTracks->Write(); | |
625 | delete fTreeGenTracks; | |
626 | fFileGenTracks->Close(); | |
627 | delete fFileGenTracks; | |
628 | return; | |
629 | } | |
630 | ||
631 | //////////////////////////////////////////////////////////////////////// | |
632 | Int_t TPCFindGenTracks::TreeKLoop() | |
633 | { | |
a7a1dd76 | 634 | /// open the file with treeK |
635 | /// loop over all entries there and save information about some tracks | |
636 | ||
e8558587 | 637 | fFileHits->cd(); |
638 | if (fDebug > 0) { | |
639 | cout<<"There are "<<fNParticles<<" primary and secondary particles in event " | |
640 | <<fEventNr<<endl; | |
641 | } | |
642 | AliStack * stack = gAlice->Stack(); | |
643 | if (!stack) {cerr<<"Stack was not found!\n"; return 1;} | |
644 | ||
645 | // not all generators give primary vertex position. Take the vertex | |
646 | // of the particle 0 as primary vertex. | |
647 | fParticle = stack->ParticleFromTreeK(0); | |
648 | fVPrim[0] = fParticle->Vx(); | |
649 | fVPrim[1] = fParticle->Vy(); | |
650 | fVPrim[2] = fParticle->Vz(); | |
651 | ||
652 | for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) { | |
653 | // for (Int_t iParticle = 0; iParticle < fDebug; iParticle++) { | |
654 | ||
655 | // load only particles with TR | |
656 | if (fIndexTR[iParticle] < 0) continue; | |
657 | fParticle = stack->ParticleFromTreeK(iParticle); | |
658 | if (fDebug > 3 && iParticle < 10) { | |
659 | cout<<"processing particle "<<iParticle<<" "; | |
660 | fParticle->Print(); | |
661 | } | |
662 | ||
663 | // fill the tree | |
664 | ||
665 | fLabel = iParticle; | |
666 | fVDist[0] = fParticle->Vx()-fVPrim[0]; | |
667 | fVDist[1] = fParticle->Vy()-fVPrim[1]; | |
668 | fVDist[2] = fParticle->Vz()-fVPrim[2]; | |
669 | fVDist[3] = TMath::Sqrt(fVDist[0]*fVDist[0]+fVDist[1]*fVDist[1]+fVDist[2]*fVDist[2]); | |
670 | fDigitRow = &(fContainerDigitRow[iParticle]); | |
671 | fRowsWithDigitsInn = fDigitRow->RowsOn(63); // 63 = number of inner rows | |
672 | fRowsWithDigits = fDigitRow->RowsOn(); | |
673 | fRowsTrackLength = fDigitRow->Last() - fDigitRow->First(); | |
674 | if (fDebug > 2 && iParticle < 10) { | |
675 | cerr<<"Fill track with a label "<<iParticle<<endl; | |
676 | } | |
677 | fDigitsInSeed = 0; | |
678 | if (fDigitRow->TestRow(seedRow11) && fDigitRow->TestRow(seedRow12)) | |
679 | fDigitsInSeed = 1; | |
680 | if (fDigitRow->TestRow(seedRow21) && fDigitRow->TestRow(seedRow22)) | |
681 | fDigitsInSeed += 10; | |
682 | ||
683 | if (fIndexTR[iParticle] >= 0) { | |
684 | fTrackRef = &(fContainerTR[fIndexTR[iParticle]]); | |
685 | // cerr<<"Debug: fTrackRef->X() = "<<fTrackRef->X()<<endl; | |
686 | } else { | |
687 | fTrackRef->SetTrack(-1); | |
688 | } | |
689 | ||
690 | fTreeGenTracks->Fill(); | |
691 | ||
692 | } | |
693 | fTreeGenTracks->AutoSave(); | |
694 | // delete gAlice; gAlice = 0; | |
695 | // fFileHits->Close(); | |
696 | ||
697 | if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl; | |
698 | ||
699 | return 0; | |
700 | } | |
701 | //////////////////////////////////////////////////////////////////////// | |
702 | Int_t TPCFindGenTracks::TreeDLoop() | |
703 | { | |
a7a1dd76 | 704 | /// open the file with treeD |
705 | /// loop over all entries there and save information about some tracks | |
e8558587 | 706 | |
707 | ||
708 | // Int_t nrow_up=fParamTPC->GetNRowUp(); | |
709 | // Int_t nrows=fParamTPC->GetNRowLow()+nrow_up; | |
710 | Int_t nInnerSector = fParamTPC->GetNInnerSector(); | |
711 | Int_t rowShift = 0; | |
712 | Int_t zero=fParamTPC->GetZeroSup(); | |
713 | // Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap); | |
714 | ||
715 | char treeDName[100]; | |
716 | sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr); | |
717 | TTree *treeD=(TTree*)fFileTreeD->Get(treeDName); | |
718 | AliSimDigits digitsAddress, *digits=&digitsAddress; | |
719 | treeD->GetBranch("Segment")->SetAddress(&digits); | |
720 | ||
721 | Int_t sectorsByRows=(Int_t)treeD->GetEntries(); | |
722 | if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl; | |
723 | for (Int_t i=0; i<sectorsByRows; i++) { | |
724 | // for (Int_t i=5720; i<sectorsByRows; i++) { | |
725 | if (!treeD->GetEvent(i)) continue; | |
726 | Int_t sec,row; | |
727 | fParamTPC->AdjustSectorRow(digits->GetID(),sec,row); | |
728 | if (fDebug > 1) cout<<sec<<' '<<row<<" \r"; | |
729 | // cerr<<sec<<' '<<row<<endl; | |
730 | ||
731 | // here I expect that upper sectors follow lower sectors | |
732 | if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow(); | |
733 | digits->First(); | |
734 | do { | |
735 | Int_t iRow=digits->CurrentRow(); | |
736 | Int_t iColumn=digits->CurrentColumn(); | |
737 | Short_t digitValue = digits->CurrentDigit(); | |
738 | // cout<<"Inner loop: sector, iRow, iColumn " | |
739 | // <<sec<<" "<<iRow<<" "<<iColumn<<endl; | |
740 | if (digitValue >= zero) { | |
741 | Int_t label; | |
742 | for (Int_t j = 0; j<3; j++) { | |
743 | label = digits->GetTrackID(iRow,iColumn,j); | |
744 | if (label >= fNParticles) { | |
745 | cerr<<"particle label too big: fNParticles, label " | |
746 | <<fNParticles<<" "<<label<<endl; | |
747 | } | |
748 | if (label >= 0 && label <= fNParticles) { | |
749 | // if (label >= 0 && label <= fDebug) { | |
750 | if (fDebug > 6 ) { | |
751 | cout<<"Inner loop: sector, iRow, iColumn, label, value, row " | |
752 | <<sec<<" " | |
753 | <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue | |
754 | <<" "<<row<<endl; | |
755 | } | |
756 | fContainerDigitRow[label].SetRow(row+rowShift); | |
757 | } | |
758 | } | |
759 | } | |
760 | } while (digits->Next()); | |
761 | } | |
762 | ||
763 | if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl; | |
764 | ||
765 | return 0; | |
766 | } | |
767 | ||
768 | //////////////////////////////////////////////////////////////////////// | |
769 | Int_t TPCFindGenTracks::TreeTRLoop() | |
770 | { | |
a7a1dd76 | 771 | /// loop over TrackReferences and store the first one for each track |
772 | ||
e8558587 | 773 | TTree *treeTR=gAlice->TreeTR(); |
774 | if (!treeTR) { | |
775 | cerr<<"TreeTR not found"<<endl; | |
776 | return 1; | |
777 | } | |
778 | Int_t nPrimaries = (Int_t) treeTR->GetEntries(); | |
779 | if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl; | |
780 | TBranch *TPCBranchTR = treeTR->GetBranch("TPC"); | |
781 | if (!TPCBranchTR) { | |
782 | cerr<<"TPC branch in TR not found"<<endl; | |
783 | return 1; | |
784 | } | |
785 | TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference"); | |
786 | TPCBranchTR->SetAddress(&TPCArrayTR); | |
787 | for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) { | |
788 | TPCBranchTR->GetEntry(iPrimPart); | |
789 | for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) { | |
790 | AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef); | |
791 | Int_t label = trackRef->GetTrack(); | |
792 | ||
793 | // save info in the fContainerTR | |
794 | if (label<0 || label > fNParticles) { | |
795 | cerr<<"Wrong label: "<<label<<endl; | |
796 | continue; | |
797 | // return 1; | |
798 | } | |
799 | ||
800 | // store only if we do not have any track reference yet for this label | |
801 | if (fIndexTR[label] >= 0) continue; | |
802 | ||
803 | // store only references with localX < fgTrackRefLocalXMax +- fgTrackRefLocalXMaxDelta | |
804 | // and the pT > fgPtCut | |
805 | Float_t localX = TR2LocalX(trackRef,fParamTPC); | |
806 | if (TMath::Abs(localX-fgTrackRefLocalXMax)>fgTrackRefLocalXMaxDelta) continue; | |
807 | if (trackRef->Pt() < fgPtCut) continue; | |
808 | ||
809 | ||
810 | // cout<<"label, xg "<<label<<" "<<xg<<endl; | |
811 | if (fLastIndexTR >= fgMaxIndexTR-1) { | |
812 | cerr<<"Too many tracks with track reference. Increase the constant" | |
813 | <<" fgMaxIndexTR"<<endl; | |
814 | return 1; | |
815 | } | |
816 | fIndexTR[label] = ++fLastIndexTR; | |
817 | ||
818 | // someone was lazy to implement copy ctor in AliTrackReference, so we have to do | |
819 | // it here "manually" | |
820 | fContainerTR[fIndexTR[label]].SetPosition(trackRef->X(),trackRef->Y(),trackRef->Z()); | |
821 | ||
822 | fContainerTR[fIndexTR[label]].SetMomentum(trackRef->Px(),trackRef->Py(),trackRef->Pz()); | |
823 | fContainerTR[fIndexTR[label]].SetTrack(trackRef->GetTrack()); | |
824 | fContainerTR[fIndexTR[label]].SetLength(trackRef->GetLength()); | |
825 | // cerr<<"Debug: trackRef->X(), stored: "<<trackRef->X()<<" " | |
826 | // << fContainerTR[fIndexTR[label]].X()<<endl; | |
827 | ||
828 | } | |
829 | } | |
830 | return 0; | |
831 | } | |
832 | //////////////////////////////////////////////////////////////////////// | |
833 | Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef, | |
834 | AliTPCParam *paramTPC) { | |
835 | ||
836 | Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; | |
837 | Int_t index[4]; | |
838 | paramTPC->Transform0to1(x,index); | |
839 | paramTPC->Transform1to2(x,index); | |
840 | return x[0]; | |
841 | } | |
842 | //////////////////////////////////////////////////////////////////////// | |
843 | ||
844 | ||
845 | ||
846 | //////////////////////////////////////////////////////////////////////// | |
847 | // | |
848 | // Start of implementation of the class TPCCmpTr | |
849 | // | |
850 | //////////////////////////////////////////////////////////////////////// | |
851 | ||
852 | class TPCCmpTr { | |
853 | ||
854 | public: | |
855 | TPCCmpTr(); | |
856 | TPCCmpTr(char* fnRecTracks, | |
857 | char* fnGenTracks ="genTracks.root", | |
858 | char* fnCmpRes ="cmpTracks.root", | |
859 | char* fnGalice ="galice.root", | |
860 | Int_t nEvents=1, Int_t firstEvent=0); | |
861 | virtual ~TPCCmpTr(); | |
862 | void Reset(); | |
863 | Int_t Exec(); | |
864 | Int_t Exec(Int_t nEvents, Int_t firstEventNr); | |
865 | void CreateTreeCmp(); | |
866 | void CloseOutputFile(); | |
867 | Bool_t ConnectGenTree(); | |
868 | Int_t TreeGenLoop(Int_t eventNr); | |
869 | Int_t TreeTLoop(Int_t eventNr); | |
870 | void SetFirstEventNr(Int_t i) {fFirstEventNr = i;} | |
871 | void SetNEvents(Int_t i) {fNEvents = i;} | |
872 | void SetDebug(Int_t level) {fDebug = level;} | |
873 | ||
874 | // tmp method, should go to TrackReferenceTPC | |
875 | Float_t TR2LocalX(AliTrackReference *trackRef, | |
876 | AliTPCParam *paramTPC); | |
877 | ||
878 | private: | |
a7a1dd76 | 879 | digitRow *fDigitRow; //!< pointer to the object saved in Branch |
880 | Int_t fEventNr; //!< current event number | |
881 | Int_t fLabel; //!< track label | |
882 | Int_t fNEvents; //!< number of events to process | |
883 | Int_t fFirstEventNr; //!< first event to process | |
e8558587 | 884 | |
a7a1dd76 | 885 | char *fFnCmp; //!< output file name with cmp tracks |
886 | TFile *fFileCmp; //!< output file with cmp tracks | |
887 | TTree *fTreeCmp; //!< output tree with cmp tracks | |
e8558587 | 888 | |
a7a1dd76 | 889 | char *fFnGenTracks; //!< input file name with gen tracks |
e8558587 | 890 | TFile *fFileGenTracks; |
891 | TTree *fTreeGenTracks; | |
892 | ||
a7a1dd76 | 893 | char *fFnHits; //!< input file name with gAlice object (needed for B) |
894 | TFile *fFileHits; //!< input file with gAlice | |
e8558587 | 895 | |
a7a1dd76 | 896 | char *fFnRecTracks; //!< input file name with tpc rec. tracks |
897 | TFile *fFileRecTracks; //!< input file with reconstructed tracks | |
898 | TTree *fTreeRecTracks; //!< tree with reconstructed tracks | |
e8558587 | 899 | |
a7a1dd76 | 900 | AliTPCtrack *fTPCTrack; //!< pointer to TPC track to connect branch |
901 | Int_t *fIndexRecTracks; //!< index of particle label in the TreeT_TPC | |
e8558587 | 902 | |
a7a1dd76 | 903 | Int_t fRowsWithDigitsInn; //!< number of rows with digits in the inner sectors |
904 | Int_t fRowsWithDigits; //!< number of rows with digits in the outer sectors | |
905 | Int_t fRowsTrackLength; //!< last - first row with digit | |
906 | Int_t fDigitsInSeed; //!< digits in the default seed rows | |
907 | TParticle *fParticle; //!< generated particle | |
908 | Double_t fVDist[4]; //!< distance of the particle vertex from primary vertex | |
e8558587 | 909 | // the fVDist[3] contains size of the 3-vector |
a7a1dd76 | 910 | AliTrackReference *fTrackRef; //!< track reference saved in the output tree |
911 | Int_t fReconstructed; //!< flag if track was reconstructed | |
912 | AliTPCParam* fParamTPC; //!< AliTPCParam | |
e8558587 | 913 | |
a7a1dd76 | 914 | Int_t fNParticles; //!< number of particles in the input tree genTracks |
915 | Int_t fDebug; //!< debug flag | |
e8558587 | 916 | |
a7a1dd76 | 917 | Int_t fNextTreeGenEntryToRead; //!< last entry already read from genTracks tree |
918 | TPCGenInfo *fGenInfo; //!< container for all the details | |
e8558587 | 919 | |
a7a1dd76 | 920 | Double_t fRecPhi; ///< reconstructed phi angle (0;2*kPI) |
921 | Double_t fLambda; ///< reconstructed | |
922 | Double_t fRecPt_1; ///< reconstructed | |
923 | Float_t fdEdx; ///< reconstructed dEdx | |
e8558587 | 924 | |
925 | ||
926 | ClassDef(TPCCmpTr,1) // class which creates and fills tree with TPCGenTrack objects | |
927 | }; | |
3a4edebe | 928 | /// \cond CLASSIMP |
e8558587 | 929 | ClassImp(TPCCmpTr) |
3a4edebe | 930 | /// \endcond |
e8558587 | 931 | |
932 | //////////////////////////////////////////////////////////////////////// | |
933 | TPCCmpTr::TPCCmpTr() | |
934 | { | |
935 | Reset(); | |
936 | } | |
937 | ||
938 | //////////////////////////////////////////////////////////////////////// | |
939 | TPCCmpTr::TPCCmpTr(char* fnRecTracks, | |
940 | char* fnGenTracks, | |
941 | char* fnCmp, | |
942 | char* fnGalice, | |
943 | Int_t nEvents, Int_t firstEvent) | |
944 | { | |
945 | Reset(); | |
946 | fFnRecTracks = fnRecTracks; | |
947 | fFnGenTracks = fnGenTracks; | |
948 | fFnCmp = fnCmp; | |
949 | fFnHits = fnGalice; | |
950 | fFirstEventNr = firstEvent; | |
951 | fEventNr = firstEvent; | |
952 | fNEvents = nEvents; | |
953 | } | |
954 | //////////////////////////////////////////////////////////////////////// | |
955 | void TPCCmpTr::Reset() | |
956 | { | |
957 | fDigitRow = 0; | |
958 | fEventNr = 0; | |
959 | fNEvents = 0; | |
960 | fTreeCmp = 0; | |
961 | fFnCmp = "cmpTracks.root"; | |
962 | fFnHits = "galice.root"; | |
963 | fFileGenTracks = 0; | |
964 | fFileHits =0; | |
965 | fParticle = 0; | |
966 | fTrackRef = 0; | |
967 | ||
968 | fRowsWithDigitsInn = 0; | |
969 | fRowsWithDigits = 0; | |
970 | fRowsTrackLength = 0; | |
971 | fDigitsInSeed = 0; | |
972 | ||
973 | fDebug = 0; | |
974 | ||
975 | fParamTPC = 0; | |
976 | fFnRecTracks = "tpc.tracks.root"; | |
977 | fTreeRecTracks = 0; | |
978 | fFileRecTracks = 0; | |
979 | fTPCTrack = 0; | |
980 | fGenInfo = 0; | |
981 | } | |
982 | //////////////////////////////////////////////////////////////////////// | |
983 | TPCCmpTr::~TPCCmpTr() | |
984 | { | |
985 | ; | |
986 | } | |
987 | //////////////////////////////////////////////////////////////////////// | |
988 | Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr) | |
989 | { | |
990 | fNEvents = nEvents; | |
991 | fFirstEventNr = firstEventNr; | |
992 | return Exec(); | |
993 | } | |
994 | ||
995 | //////////////////////////////////////////////////////////////////////// | |
996 | Int_t TPCCmpTr::Exec() | |
997 | { | |
998 | TStopwatch timer; | |
999 | timer.Start(); | |
1000 | ||
1001 | fDigitRow = new digitRow(); | |
1002 | CreateTreeCmp(); | |
1003 | if (!fTreeCmp) { | |
1004 | cerr<<"output tree not created"<<endl; | |
1005 | return 1; | |
1006 | } | |
1007 | fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here | |
1008 | if (!fFileHits) { | |
1009 | cerr<<"Cannot open file with gAlice object "<<fFnHits<<endl; | |
1010 | return 1; | |
1011 | } | |
1012 | fFileHits->cd(); | |
1013 | SetFieldFactor(); | |
1014 | ||
1015 | fParamTPC = LoadTPCParam(fFileHits); | |
1016 | if (!fParamTPC) { | |
1017 | cerr<<"TPC parameters not found and could not be created"<<endl; | |
1018 | return 1; | |
1019 | } | |
1020 | ||
1021 | if (!ConnectGenTree()) { | |
1022 | cerr<<"Cannot connect tree with generated tracks"<<endl; | |
1023 | return 1; | |
1024 | } | |
1025 | fFileHits->cd(); | |
1026 | fNextTreeGenEntryToRead = 0; | |
1027 | cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl; | |
1028 | for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents; | |
1029 | eventNr++) { | |
1030 | fNParticles = gAlice->GetEvent(fEventNr); | |
1031 | fIndexRecTracks = new Int_t[fNParticles]; | |
1032 | for (Int_t i = 0; i<fNParticles; i++) { | |
1033 | fIndexRecTracks[i] = -1; | |
1034 | } | |
1035 | ||
1036 | cout<<"Start to process event "<<fEventNr<<endl; | |
1037 | cout<<"\tfNParticles = "<<fNParticles<<endl; | |
1038 | if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl; | |
1039 | if (TreeTLoop(eventNr)>0) return 1; | |
1040 | ||
1041 | if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl; | |
1042 | if (TreeGenLoop(eventNr)>0) return 1; | |
1043 | if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl; | |
1044 | ||
1045 | delete fTreeRecTracks; | |
1046 | ||
1047 | delete [] fIndexRecTracks; | |
1048 | } | |
1049 | ||
1050 | CloseOutputFile(); | |
1051 | ||
1052 | cerr<<"Exec finished"<<endl; | |
1053 | fFileHits->cd(); | |
1054 | delete gAlice; | |
1055 | fFileHits->Close(); | |
1056 | delete fFileHits; | |
1057 | ||
1058 | timer.Stop(); | |
1059 | timer.Print(); | |
1060 | return 0; | |
1061 | } | |
1062 | //////////////////////////////////////////////////////////////////////// | |
1063 | Bool_t TPCCmpTr::ConnectGenTree() | |
1064 | { | |
a7a1dd76 | 1065 | /// connect all branches from the genTracksTree |
1066 | /// use the same variables as for the new cmp tree, it may work | |
1067 | ||
e8558587 | 1068 | fFileGenTracks = TFile::Open(fFnGenTracks,"READ"); |
1069 | if (!fFileGenTracks) { | |
1070 | cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl; | |
1071 | return kFALSE; | |
1072 | } | |
1073 | fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree"); | |
1074 | if (!fTreeGenTracks) { | |
1075 | cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file " | |
1076 | <<fFnGenTracks<<endl; | |
1077 | return kFALSE; | |
1078 | } | |
1079 | fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr); | |
1080 | fTreeGenTracks->SetBranchAddress("fLabel",&fLabel); | |
1081 | fTreeGenTracks->SetBranchAddress("fRowsWithDigitsInn",&fRowsWithDigitsInn); | |
1082 | fTreeGenTracks->SetBranchAddress("fRowsWithDigits",&fRowsWithDigits); | |
1083 | fTreeGenTracks->SetBranchAddress("fRowsTrackLength",&fRowsTrackLength); | |
1084 | fTreeGenTracks->SetBranchAddress("fDigitsInSeed",&fDigitsInSeed); | |
1085 | fTreeGenTracks->SetBranchAddress("Particle",&fParticle); | |
1086 | fTreeGenTracks->SetBranchAddress("fVDist",fVDist); | |
1087 | fTreeGenTracks->SetBranchAddress("TR",&fTrackRef); | |
1088 | ||
1089 | if (fDebug > 1) { | |
1090 | cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl; | |
1091 | } | |
1092 | return kTRUE; | |
1093 | } | |
1094 | ||
1095 | ||
1096 | //////////////////////////////////////////////////////////////////////// | |
1097 | void TPCCmpTr::CreateTreeCmp() | |
1098 | { | |
1099 | fFileCmp = TFile::Open(fFnCmp,"RECREATE"); | |
1100 | if (!fFileCmp) { | |
1101 | cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl; | |
1102 | return; | |
1103 | } | |
1104 | fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks"); | |
1105 | TBranch *branchBits = fTreeCmp->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0); | |
1106 | if (!branchBits) { | |
1107 | cerr<<"Error in CreateTreeCmp: cannot create branch."<<endl; | |
1108 | return; | |
1109 | } | |
1110 | fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I"); | |
1111 | fTreeCmp->Branch("fLabel",&fLabel,"fLabel/I"); | |
1112 | fTreeCmp->Branch("fRowsWithDigitsInn",&fRowsWithDigitsInn,"fRowsWithDigitsInn/I"); | |
1113 | fTreeCmp->Branch("fRowsWithDigits",&fRowsWithDigits,"fRowsWithDigits/I"); | |
1114 | fTreeCmp->Branch("fRowsTrackLength",&fRowsTrackLength,"fRowsTrackLength/I"); | |
1115 | fTreeCmp->Branch("fDigitsInSeed",&fDigitsInSeed,"fDigitsInSeed/I"); | |
1116 | ||
1117 | fTreeCmp->Branch("fReconstructed",&fReconstructed,"fReconstructed/I"); | |
1118 | fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack); | |
1119 | ||
1120 | fTreeCmp->Branch("Particle","TParticle",&fParticle); | |
1121 | fTreeCmp->Branch("fVDist",&fVDist,"fVDist[4]/D"); | |
1122 | fTreeCmp->Branch("TR","AliTrackReference",&fTrackRef); | |
1123 | ||
1124 | fTreeCmp->AutoSave(); | |
1125 | } | |
1126 | //////////////////////////////////////////////////////////////////////// | |
1127 | void TPCCmpTr::CloseOutputFile() | |
1128 | { | |
1129 | if (!fFileCmp) { | |
1130 | cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl; | |
1131 | return; | |
1132 | } | |
1133 | fFileCmp->cd(); | |
1134 | fTreeCmp->Write(); | |
1135 | delete fTreeCmp; | |
1136 | fFileCmp->Close(); | |
1137 | delete fFileCmp; | |
1138 | return; | |
1139 | } | |
1140 | //////////////////////////////////////////////////////////////////////// | |
1141 | ||
1142 | Float_t TPCCmpTr::TR2LocalX(AliTrackReference *trackRef, | |
1143 | AliTPCParam *paramTPC) { | |
1144 | ||
1145 | Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; | |
1146 | Int_t index[4]; | |
1147 | paramTPC->Transform0to1(x,index); | |
1148 | paramTPC->Transform1to2(x,index); | |
1149 | return x[0]; | |
1150 | } | |
1151 | //////////////////////////////////////////////////////////////////////// | |
1152 | ||
1153 | Int_t TPCCmpTr::TreeTLoop(Int_t eventNr) | |
1154 | { | |
a7a1dd76 | 1155 | /// loop over all TPC reconstructed tracks and store info in memory |
e8558587 | 1156 | |
1157 | ||
1158 | if (!fFileRecTracks) fFileRecTracks = TFile::Open(fFnRecTracks,"read"); | |
1159 | if (!fFileRecTracks->IsOpen()) { | |
1160 | cerr<<"Cannot open file "<<fFnRecTracks<<endl; | |
1161 | return 1; | |
1162 | } | |
1163 | ||
1164 | char treeNameBase[11] = "TreeT_TPC_"; | |
1165 | char treeName[20]; | |
1166 | sprintf(treeName,"%s%d",treeNameBase,eventNr); | |
1167 | ||
1168 | fTreeRecTracks=(TTree*)fFileRecTracks->Get(treeName); | |
1169 | if (!fTreeRecTracks) { | |
1170 | cerr<<"Can't get a tree with TPC rec. tracks named "<<treeName<<endl; | |
1171 | return 1; | |
1172 | } | |
1173 | ||
1174 | Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries(); | |
1175 | if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" " | |
1176 | <<nEntries<<endl; | |
1177 | TBranch * br= fTreeRecTracks->GetBranch("tracks"); | |
1178 | br->SetAddress(&fTPCTrack); | |
1179 | ||
1180 | for (Int_t iEntry=0; iEntry<nEntries;iEntry++){ | |
1181 | br->GetEntry(iEntry); | |
1182 | Int_t label = fTPCTrack->GetLabel(); | |
1183 | fIndexRecTracks[label] = iEntry; | |
1184 | } | |
1185 | ||
1186 | if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl; | |
1187 | ||
1188 | return 0; | |
1189 | } | |
1190 | //////////////////////////////////////////////////////////////////////// | |
1191 | Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr) | |
1192 | { | |
a7a1dd76 | 1193 | /// loop over all entries for a given event, find corresponding |
1194 | /// rec. track and store in the fTreeCmp | |
1195 | ||
e8558587 | 1196 | fFileGenTracks->cd(); |
1197 | Int_t entry = fNextTreeGenEntryToRead; | |
1198 | Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast(); | |
1199 | cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" " | |
1200 | <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl; | |
1201 | while (entry < nParticlesTR) { | |
1202 | fTreeGenTracks->GetEntry(entry); | |
1203 | entry++; | |
1204 | if (fEventNr < eventNr) continue; | |
1205 | if (fEventNr > eventNr) break; | |
1206 | fNextTreeGenEntryToRead = entry-1; | |
1207 | if (fDebug > 2 && fLabel < 10) { | |
1208 | cerr<<"Fill track with a label "<<fLabel<<endl; | |
1209 | } | |
1210 | ||
1211 | fReconstructed = 0; | |
1212 | fdEdx = 0.; | |
1213 | fRecPhi = fLambda = fRecPt_1 = 0.; | |
1214 | ||
1215 | if (fDebug > 2) { | |
1216 | cerr<<"fLabel, fIndexRecTracks[fLabel] "<<fLabel<<" "<<fIndexRecTracks[fLabel]<<endl; | |
1217 | } | |
1218 | if (fIndexRecTracks[fLabel] >= 0) { | |
1219 | Int_t nBytes = fTreeRecTracks->GetEvent(fIndexRecTracks[fLabel]); | |
1220 | if (nBytes > 0) { | |
1221 | fReconstructed = 1; | |
1222 | fdEdx = fTPCTrack->GetdEdx(); | |
1223 | Double_t Localx = TR2LocalX(fTrackRef,fParamTPC); | |
1224 | if (fDebug > 3) { | |
1225 | cerr<<"Track local X before prolongation: "<<fTPCTrack->GetX()<<endl; | |
1226 | } | |
1227 | fTPCTrack->PropagateTo(Localx); | |
1228 | Double_t par[5]; | |
1229 | if (fDebug > 3) { | |
1230 | cerr<<"Track local X after prolongation: "<<fTPCTrack->GetX()<<endl; | |
1231 | } | |
1232 | fTPCTrack->GetExternalParameters(Localx,par); | |
1233 | fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha(); | |
d8d3b5b8 | 1234 | if (fRecPhi<0) fRecPhi+=2*TMath::Pi(); |
1235 | if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi(); | |
e8558587 | 1236 | // fRecPhi = (fRecPhi)*kRaddeg; |
1237 | fLambda = TMath::ATan(par[3]); | |
1238 | fRecPt_1 = TMath::Abs(par[4]); | |
1239 | } | |
1240 | } | |
1241 | ||
1242 | ||
1243 | fTreeCmp->Fill(); | |
1244 | ||
1245 | } | |
1246 | fTreeCmp->AutoSave(); | |
1247 | // delete gAlice; gAlice = 0; | |
1248 | // fFileHits->Close(); | |
1249 | ||
1250 | if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl; | |
1251 | ||
1252 | return 0; | |
1253 | } | |
1254 | //////////////////////////////////////////////////////////////////////// |