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