New TPC comparison added
[u/mrichter/AliRoot.git] / TPC / AliTPCComparisonMI.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "iostream.h"
3 #include <stdio.h>
4 #include <string.h>
5 #include "Rtypes.h"
6 #include "TFile.h"
7 #include "TTree.h"
8 #include "TString.h"
9 #include "TBenchmark.h"
10 #include "TStopwatch.h"
11 #include "TParticle.h"
12 #include "AliRun.h"
13 #include "AliStack.h"
14 #include "AliTPCtrack.h"
15 #include "AliSimDigits.h"
16 #include "AliTPCParam.h"
17 #include "TParticle.h"
18 #include "AliTPC.h"
19 #include "AliDetector.h"
20 #include "AliTrackReference.h"
21 #include "TSystem.h"
22 #include "TTimer.h"
23 #include "TVector3.h"
24 //
25 #include "iostream.h"  
26 #include "Rtypes.h"
27 #include "TSystem.h"
28 #include "TTimer.h"
29 #include "Getline.h"
30 #include "TChain.h"
31 #include "TString.h"
32 #include "TFile.h"
33 #include "AliRun.h"
34 #include "AliTPCComparisonMI.h"
35 #include "AliTPCParamSR.h"
36 #include "AliTracker.h"
37 #include "AliComplexCluster.h"
38
39
40 #endif
41  
42 //
43 //
44 Bool_t ImportgAlice(TFile *file);
45 TFile* OpenAliceFile(char *fn, Bool_t importgAlice = kFALSE, char *mode = "read");
46
47 AliTPCParam * GetTPCParam(){
48   AliTPCParamSR * par = new AliTPCParamSR;
49   par->Update();
50   return par;
51 }
52
53
54
55 ////////////////////////////////////////////////////////////////////////
56 Bool_t ImportgAlice(TFile *file) {
57 // read in gAlice object from the file
58   gAlice = (AliRun*)file->Get("gAlice");
59   if (!gAlice)  return kFALSE;
60   return kTRUE;
61 }
62 ////////////////////////////////////////////////////////////////////////
63 TFile* OpenAliceFile(char *fn, Bool_t importgAlice, char *mode) {
64   TFile *file = TFile::Open(fn,mode);
65   if (!file->IsOpen()) {
66     cerr<<"OpenAliceFile: cannot open file "<<fn<<" in mode "
67         <<mode<<endl;
68     return 0;
69   }
70   if (!importgAlice) return file;
71   if (ImportgAlice(file)) return file;
72   return 0;
73 }
74 ////////////////////////////////////////////////////////////////////////
75
76 //_____________________________________________________________________________
77 Float_t TPCBetheBloch(Float_t bg)
78 {
79   //
80   // Bethe-Bloch energy loss formula
81   //
82   const Double_t kp1=0.76176e-1;
83   const Double_t kp2=10.632;
84   const Double_t kp3=0.13279e-4;
85   const Double_t kp4=1.8631;
86   const Double_t kp5=1.9479;
87
88   Double_t dbg = (Double_t) bg;
89
90   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
91
92   Double_t aa = TMath::Power(beta,kp4);
93   Double_t bb = TMath::Power(1./dbg,kp5);
94
95   bb=TMath::Log(kp3+bb);
96   
97   return ((Float_t)((kp2-aa-bb)*kp1/aa));
98 }
99
100
101
102
103
104 ////////////////////////////////////////////////////////////////////////
105 AliTPCGenInfo::AliTPCGenInfo()
106 {
107   fReferences  = 0;
108 }
109
110 AliTPCGenInfo::~AliTPCGenInfo()
111 {
112   if (fReferences) delete fReferences;
113   fReferences =0;
114 }
115 ////////////////////////////////////////////////////////////////////////
116
117 ////////////////////////////////////////////////////////////////////////
118 //
119 // End of implementation of the class AliTPCGenInfo
120 //
121 ////////////////////////////////////////////////////////////////////////
122
123
124
125 ////////////////////////////////////////////////////////////////////////
126 digitRow::digitRow()
127 {
128   Reset();
129 }
130 ////////////////////////////////////////////////////////////////////////
131 digitRow & digitRow::operator=(const digitRow &digOld)
132 {
133   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
134   return (*this);
135 }
136 ////////////////////////////////////////////////////////////////////////
137 void digitRow::SetRow(Int_t row) 
138 {
139   if (row >= 8*kgRowBytes) {
140     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
141     return;
142   }
143   Int_t iC = row/8;
144   Int_t iB = row%8;
145   SETBIT(fDig[iC],iB);
146 }
147
148 ////////////////////////////////////////////////////////////////////////
149 Bool_t digitRow::TestRow(Int_t row)
150 {
151 //
152 // return kTRUE if row is on
153 //
154   Int_t iC = row/8;
155   Int_t iB = row%8;
156   return TESTBIT(fDig[iC],iB);
157 }
158 ////////////////////////////////////////////////////////////////////////
159 Int_t digitRow::RowsOn(Int_t upto)
160 {
161 //
162 // returns number of rows with a digit  
163 // count only rows less equal row number upto
164 //
165   Int_t total = 0;
166   for (Int_t i = 0; i<kgRowBytes; i++) {
167     for (Int_t j = 0; j < 8; j++) {
168       if (i*8+j > upto) return total;
169       if (TESTBIT(fDig[i],j))  total++;
170     }
171   }
172   return total;
173 }
174 ////////////////////////////////////////////////////////////////////////
175 void digitRow::Reset()
176 {
177 //
178 // resets all rows to zero
179 //
180   for (Int_t i = 0; i<kgRowBytes; i++) {
181     fDig[i] <<= 8;
182   }
183 }
184 ////////////////////////////////////////////////////////////////////////
185 Int_t digitRow::Last()
186 {
187 //
188 // returns the last row number with a digit
189 // returns -1 if now digits 
190 //
191   for (Int_t i = kgRowBytes-1; i>=0; i--) {
192     for (Int_t j = 7; j >= 0; j--) {
193       if TESTBIT(fDig[i],j) return i*8+j;
194     }
195   }
196   return -1;
197 }
198 ////////////////////////////////////////////////////////////////////////
199 Int_t digitRow::First()
200 {
201 //
202 // returns the first row number with a digit
203 // returns -1 if now digits 
204 //
205   for (Int_t i = 0; i<kgRowBytes; i++) {
206     for (Int_t j = 0; j < 8; j++) {
207       if (TESTBIT(fDig[i],j)) return i*8+j;
208     }
209   }
210   return -1;
211 }
212
213 ////////////////////////////////////////////////////////////////////////
214 //
215 // end of implementation of a class digitRow
216 //
217 ////////////////////////////////////////////////////////////////////////
218   
219 ////////////////////////////////////////////////////////////////////////
220 TPCFindGenTracks::TPCFindGenTracks()
221 {
222   fMCInfo = new AliTPCGenInfo;
223   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
224   Reset();
225 }
226
227 ////////////////////////////////////////////////////////////////////////
228 TPCFindGenTracks::TPCFindGenTracks(char * fnGalice, char* fnRes,
229                                    Int_t nEvents, Int_t firstEvent)
230 {
231   Reset();
232   fFirstEventNr = firstEvent;
233   fEventNr = firstEvent;
234   fNEvents = nEvents;
235   fFnRes = fnRes;
236   //
237   fLoader = AliRunLoader::Open(fnGalice);
238   if (gAlice){
239     delete gAlice->GetRunLoader();
240     delete gAlice;
241     gAlice = 0x0;
242   }
243   if (fLoader->LoadgAlice()){
244     cerr<<"Error occured while l"<<endl;
245   }
246   Int_t nall = fLoader->GetNumberOfEvents();
247   if (nall<=0){
248     cerr<<"no events available"<<endl;
249     fEventNr = 0;
250     return;
251   }
252   if (firstEvent+nEvents>nall) {
253     fEventNr = nall-firstEvent;
254     cerr<<"restricted number of events availaible"<<endl;
255   }
256 }
257
258 ////////////////////////////////////////////////////////////////////////
259 void TPCFindGenTracks::Reset()
260 {
261   fEventNr = 0;
262   fNEvents = 0;
263   fTreeGenTracks = 0;
264   fFnRes = "genTracks.root";
265   fFileGenTracks = 0;
266   fContainerDigitRow = 0;
267   //
268   fReferences = 0;
269   fReferenceIndex0 = 0;
270   fReferenceIndex1 = 0;
271   //
272   fDebug = 0;
273   fVPrim[0] = -1000.;
274   fVPrim[1] = -1000.;
275   fVPrim[2] = -1000.;
276   fParamTPC = 0;
277 }
278 ////////////////////////////////////////////////////////////////////////
279 TPCFindGenTracks::~TPCFindGenTracks()
280 {
281   ;
282 }
283
284
285 Int_t  TPCFindGenTracks::SetIO()
286 {
287   //
288   // 
289   CreateTreeGenTracks();
290   if (!fTreeGenTracks) return 1;
291   AliTracker::SetFieldFactor(); 
292   fParamTPC = GetTPCParam();
293   //
294   return 0;
295 }
296
297 ////////////////////////////////////////////////////////////////////////
298 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
299 {
300   //
301   // 
302   // SET INPUT
303   gAlice->GetEvent(eventNr);
304   fLoader->SetEventNumber(eventNr);
305   //
306   fLoader->LoadHeader();
307   fLoader->LoadKinematics();  
308   fStack = fLoader->Stack();
309   //
310   fLoader->LoadTrackRefs();
311   fTreeTR = fLoader->TreeTR();
312   //
313   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
314   tpcl->LoadDigits();
315   fTreeD = tpcl->TreeD();
316   
317   return 0;
318 }
319
320 ////////////////////////////////////////////////////////////////////////
321 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
322 {
323   fNEvents = nEvents;
324   fFirstEventNr = firstEventNr;
325   return Exec();
326 }
327
328 ////////////////////////////////////////////////////////////////////////
329 Int_t TPCFindGenTracks::Exec()  
330 {
331   TStopwatch timer;
332   timer.Start();
333   Int_t status =SetIO();
334   if (status>0) return status;
335   //
336
337   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
338        fEventNr++) {
339     SetIO(fEventNr);
340     fNParticles = fStack->GetNtrack();
341     fContainerDigitRow = new digitRow[fNParticles];
342     //
343     fReferences = new AliTrackReference[fgMaxTR];
344     fReferenceIndex0 = new Int_t[fNParticles];
345     fReferenceIndex1 = new Int_t[fNParticles];
346
347     for (Int_t i = 0; i<fNParticles; i++) {
348       fReferenceIndex0[i] = -1;
349       fReferenceIndex1[i] = -1;
350     }
351     
352     cout<<"Start to process event "<<fEventNr<<endl;
353     cout<<"\tfNParticles = "<<fNParticles<<endl;
354     if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
355     if (TreeDLoop()>0) return 1;
356     if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
357     if (TreeTRLoop()>0) return 1;
358     if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
359     if (TreeKLoop()>0) return 1;
360     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
361
362     delete [] fContainerDigitRow;
363     delete [] fReferences;
364     delete [] fReferenceIndex0;
365     delete [] fReferenceIndex1;    
366   }
367
368   CloseOutputFile();
369
370   cerr<<"Exec finished"<<endl;
371   delete gAlice;
372
373   timer.Stop();
374   timer.Print();
375   return 0;
376 }
377 ////////////////////////////////////////////////////////////////////////
378 void TPCFindGenTracks::CreateTreeGenTracks() 
379 {
380   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
381   if (!fFileGenTracks) {
382     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
383     return;
384   }
385   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
386   
387
388
389   fMCInfo = new AliTPCGenInfo;
390   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
391
392   fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
393
394   fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
395
396   fTreeGenTracks->AutoSave();
397 }
398 ////////////////////////////////////////////////////////////////////////
399 void TPCFindGenTracks::CloseOutputFile() 
400 {
401   if (!fFileGenTracks) {
402     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
403     return;
404   }
405   fFileGenTracks->cd();
406   fTreeGenTracks->Write();  
407   delete fTreeGenTracks;
408   fFileGenTracks->Close();
409   delete fFileGenTracks;
410   return;
411 }
412
413 ////////////////////////////////////////////////////////////////////////
414 Int_t TPCFindGenTracks::TreeKLoop()
415 {
416 //
417 // open the file with treeK
418 // loop over all entries there and save information about some tracks
419 //
420
421   AliStack * stack = fStack;
422   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
423   
424   if (fDebug > 0) {
425     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
426         <<fEventNr<<endl;
427   }
428   
429   Int_t  ipdg = 0;
430   TParticlePDG *ppdg = 0;
431
432  
433   // not all generators give primary vertex position. Take the vertex
434   // of the particle 0 as primary vertex.
435   TDatabasePDG  pdg; //get pdg table
436   
437   //thank you very much root for this
438   TBranch * br = fTreeGenTracks->GetBranch("MC");
439
440   TParticle *particle = stack->ParticleFromTreeK(0);
441   fVPrim[0] = particle->Vx();
442   fVPrim[1] = particle->Vy();
443   fVPrim[2] = particle->Vz();
444   for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
445
446 // load only particles with TR
447     if (fReferenceIndex0[iParticle]<0) continue;
448     //////////////////////////////////////////////////////////////////////
449     fMCInfo->fLabel = iParticle;
450
451     fMCInfo->fRow = (fContainerDigitRow[iParticle]);
452     fMCInfo->fRowsWithDigits = fMCInfo->fRow.RowsOn();    
453     if (fMCInfo->fRowsWithDigits<10) continue;
454     fMCInfo->fRowsWithDigitsInn = fMCInfo->fRow.RowsOn(63); // 63 = number of inner rows
455     fMCInfo->fRowsTrackLength = fMCInfo->fRow.Last() - fMCInfo->fRow.First();
456     fMCInfo->fDigitsInSeed = 0;
457     if (fMCInfo->fRow.TestRow(seedRow11) && fMCInfo->fRow.TestRow(seedRow12)) 
458       fMCInfo->fDigitsInSeed = 1;
459     if (fMCInfo->fRow.TestRow(seedRow21) && fMCInfo->fRow.TestRow(seedRow22)) 
460       fMCInfo->fDigitsInSeed += 10;
461     //
462     //
463     //
464     fMCInfo->fParticle = *(stack->Particle(iParticle));
465     if (fMCInfo->fParticle.GetLastDaughter()>0&&fMCInfo->fParticle.GetLastDaughter()<fNParticles){
466       TParticle * dparticle0 = 0;
467       if (fMCInfo->fParticle.GetDaughter(0)>0) dparticle0 = stack->Particle(fMCInfo->fParticle.GetDaughter(0));
468       TParticle * dparticle1 = 0;
469       if (fMCInfo->fParticle.GetDaughter(1)>0) dparticle1 = stack->Particle(fMCInfo->fParticle.GetDaughter(1));
470       if ((dparticle1)&&(dparticle1->P()>0.15)) {
471         fMCInfo->fDecayCoord[0] = dparticle1->Vx();
472         fMCInfo->fDecayCoord[1] = dparticle1->Vy();
473         fMCInfo->fDecayCoord[2] = dparticle1->Vz();
474       }
475       else 
476         fMCInfo->fDecayCoord[0]=1000000;
477     }else
478       fMCInfo->fDecayCoord[0]=1000000;
479     fMCInfo->fParticle = *(stack->ParticleFromTreeK(iParticle));
480     //
481     //
482     //
483     fMCInfo->fLabel = iParticle;
484     fMCInfo->fVDist[0] = fMCInfo->fParticle.Vx()-fVPrim[0];
485     fMCInfo->fVDist[1] = fMCInfo->fParticle.Vy()-fVPrim[1];
486     fMCInfo->fVDist[2] = fMCInfo->fParticle.Vz()-fVPrim[2]; 
487     fMCInfo->fVDist[3] = TMath::Sqrt(fMCInfo->fVDist[0]*fMCInfo->fVDist[0]+
488                                       fMCInfo->fVDist[1]*fMCInfo->fVDist[1]+fMCInfo->fVDist[2]*fMCInfo->fVDist[2]);
489
490     //
491     Int_t index = fReferenceIndex0[iParticle];
492     AliTrackReference  ref  = fReferences[index];
493     // if (ref.GetTrack()!=iParticle)
494     //  printf("problem2\n");
495     fMCInfo->fTrackRef = ref;
496   
497     Int_t rfindex =0;
498     if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
499     fMCInfo->fReferences = 
500       new TClonesArray("AliTrackReference");
501     fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
502     //
503     for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
504       AliTrackReference  ref  = fReferences[i];
505       AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
506       if (ref.GetTrack()!=iParticle)
507         printf("problem5\n");   
508       *ref2 = ref;
509       rfindex++;
510     }   
511     //
512     //
513     ipdg = fMCInfo->fParticle.GetPdgCode();
514     ppdg = pdg.GetParticle(ipdg);
515     // calculate primary ionization per cm
516     if (ppdg){
517       Float_t mass = ppdg->Mass();
518       Float_t p = TMath::Sqrt(fMCInfo->fTrackRef.Px()*fMCInfo->fTrackRef.Px()+
519                               fMCInfo->fTrackRef.Py()*fMCInfo->fTrackRef.Py()+
520                               fMCInfo->fTrackRef.Pz()*fMCInfo->fTrackRef.Pz());
521       
522       //      Float_t betagama = fMCInfo->fParticle.P()/mass;
523       Float_t betagama = p /mass;
524       fMCInfo->fPrim = TPCBetheBloch(betagama);
525     }
526     //////////////////////////////////////////////////////////////////////
527     br->SetAddress(&fMCInfo);
528     fTreeGenTracks->Fill();
529
530   }
531   fTreeGenTracks->AutoSave();
532
533   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
534
535   return 0;
536 }
537
538
539
540
541 ////////////////////////////////////////////////////////////////////////
542 Int_t TPCFindGenTracks::TreeDLoop()
543 {
544 //
545 // open the file with treeD
546 // loop over all entries there and save information about some tracks
547 //
548
549   Int_t nInnerSector = fParamTPC->GetNInnerSector();
550   Int_t rowShift = 0;
551   //  Int_t zero=fParamTPC->GetZeroSup();
552   Int_t zero=fParamTPC->GetZeroSup()+6;
553
554   
555   //char treeDName[100]; 
556   //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
557   //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
558   //
559   //
560   AliSimDigits digitsAddress, *digits=&digitsAddress;
561   fTreeD->GetBranch("Segment")->SetAddress(&digits);
562
563   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
564   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
565   for (Int_t i=0; i<sectorsByRows; i++) {
566 //  for (Int_t i=5720; i<sectorsByRows; i++) {
567     if (!fTreeD->GetEvent(i)) continue;
568     Int_t sec,row;
569     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
570     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
571 //    cerr<<sec<<' '<<row<<endl;
572
573 // here I expect that upper sectors follow lower sectors
574     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
575     //
576     //digits->ExpandTrackBuffer();     
577     //Int_t *tracks = digits->GetTracks();
578     digits->First();    
579     do {
580       Int_t iRow=digits->CurrentRow();
581       Int_t iColumn=digits->CurrentColumn();
582       Short_t digitValue = digits->CurrentDigit();
583 //      cout<<"Inner loop: sector, iRow, iColumn "
584 //        <<sec<<" "<<iRow<<" "<<iColumn<<endl;
585       if (digitValue >= zero) {
586         Int_t label;
587         for (Int_t j = 0; j<3; j++) {
588           label = digits->GetTrackID(iRow,iColumn,j); 
589           //label = digits->GetTrackIDFast(iRow,iColumn,j); 
590           
591           if (label >= fNParticles) { //don't label from bakground event
592             continue;
593           }
594           if (label >= 0 && label <= fNParticles) {
595 //        if (label >= 0 && label <= fDebug) {
596             if (fDebug > 6 ) {
597               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
598                   <<sec<<" "
599                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
600                   <<" "<<row<<endl;
601             }   
602             fContainerDigitRow[label].SetRow(row+rowShift);
603           }
604         }
605       }
606     } while (digits->Next());
607   }
608
609   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
610
611   return 0;
612 }
613
614
615 ////////////////////////////////////////////////////////////////////////
616 Int_t TPCFindGenTracks::TreeTRLoop()
617 {
618   //
619   // loop over TrackReferences and store the first one for each track
620   //
621   
622   TTree * treeTR = fTreeTR;
623   if (!treeTR) {
624     cerr<<"TreeTR not found"<<endl;
625     return 1;
626   }
627   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
628   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
629   //
630   //
631   TBranch *TPCBranchTR  = treeTR->GetBranch("TPC");
632   if (!TPCBranchTR) {
633     cerr<<"TPC branch in TR not found"<<endl;
634     return 1;
635   }
636   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
637   TPCBranchTR->SetAddress(&TPCArrayTR);
638   //
639   //
640   //
641   Int_t index     =  0;
642   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
643     TPCBranchTR->GetEntry(iPrimPart);
644     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
645       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);      
646       //
647       if (trackRef->Pt() < fgPtCut) continue;      
648       Int_t label = trackRef->GetTrack();
649       if (label<0  || label > fNParticles) {
650         continue;
651       }
652       if (index>=fgMaxTR) continue;     //restricted size of buffer for TR
653       fReferences[index] = *trackRef;
654       fReferenceIndex1[label] = index;  // the last ref with given label
655       if (fReferenceIndex0[label]==-1) fReferenceIndex0[label] = index;   //the first index with label
656       index++;
657            
658     }
659   }
660   return 0;
661 }
662
663 ////////////////////////////////////////////////////////////////////////
664 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
665                                     AliTPCParam *paramTPC) {
666
667   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
668   Int_t index[4];
669   paramTPC->Transform0to1(x,index);
670   paramTPC->Transform1to2(x,index);
671   return x[0];
672 }
673 ////////////////////////////////////////////////////////////////////////
674
675
676
677   
678 ////////////////////////////////////////////////////////////////////////
679 TPCCmpTr::TPCCmpTr()
680 {
681   Reset();
682 }
683
684 ////////////////////////////////////////////////////////////////////////
685 TPCCmpTr::TPCCmpTr(char* fnGenTracks,
686                    char* fnCmp,
687                    char* fnGalice,
688                    Int_t nEvents, Int_t firstEvent)
689 {
690   Reset();
691   fFnGenTracks = fnGenTracks;
692   fFnCmp = fnCmp;
693   fFirstEventNr = firstEvent;
694   fEventNr = firstEvent;
695   fNEvents = nEvents;
696   
697   //
698   fLoader = AliRunLoader::Open(fnGalice);
699   if (gAlice){
700     delete gAlice->GetRunLoader();
701     delete gAlice;
702     gAlice = 0x0;
703   }
704   if (fLoader->LoadgAlice()){
705     cerr<<"Error occured while l"<<endl;
706   }
707   Int_t nall = fLoader->GetNumberOfEvents();
708   if (nall<=0){
709     cerr<<"no events available"<<endl;
710     fEventNr = 0;
711     return;
712   }
713   if (firstEvent+nEvents>nall) {
714     fEventNr = nall-firstEvent;
715     cerr<<"restricted number of events availaible"<<endl;
716   }
717 }
718
719 //////////////////////////////////////////////////////////////
720 Int_t TPCCmpTr::SetIO()
721 {
722   //
723   // 
724   CreateTreeCmp();
725   if (!fTreeCmp) return 1;
726   AliTracker::SetFieldFactor(); 
727   fParamTPC = GetTPCParam();
728   //
729   if (!ConnectGenTree()) {
730     cerr<<"Cannot connect tree with generated tracks"<<endl;
731     return 1;
732   }
733   return 0;
734 }
735
736 //////////////////////////////////////////////////////////////
737
738 Int_t TPCCmpTr::SetIO(Int_t eventNr)
739 {
740   //
741   // 
742   // SET INPUT
743   gAlice->GetEvent(eventNr);
744   fLoader->SetEventNumber(eventNr);  
745   //
746   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
747   tpcl->LoadTracks();
748   fTreeRecTracks = tpcl->TreeT();
749   //
750   return 0;
751 }
752
753
754
755 ////////////////////////////////////////////////////////////////////////
756 void TPCCmpTr::Reset()
757 {
758   fEventNr = 0;
759   fNEvents = 0;
760   fTreeCmp = 0;
761   fFnCmp = "cmpTracks.root";
762   fFileGenTracks = 0;
763   fDebug = 0;
764   //
765   fParamTPC = 0;
766   fTreeRecTracks = 0;
767   fTreePoints =0;
768
769   fTPCTrack = 0; 
770   fTracks   = 0;
771   fTrackPoints =0;
772 }
773 ////////////////////////////////////////////////////////////////////////
774 TPCCmpTr::~TPCCmpTr()
775 {
776   ;
777 }
778 ////////////////////////////////////////////////////////////////////////
779 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
780 {
781   fNEvents = nEvents;
782   fFirstEventNr = firstEventNr;
783   return Exec();
784 }
785
786 ////////////////////////////////////////////////////////////////////////
787 Int_t TPCCmpTr::Exec()
788 {
789   TStopwatch timer;
790   timer.Start();
791
792   if (SetIO()==1) 
793     return 1;
794    
795   fNextTreeGenEntryToRead = 0;
796   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
797   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
798        eventNr++) {
799     SetIO(fEventNr);
800     fNParticles = gAlice->GetEvent(fEventNr);    
801
802     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
803     fFakeRecTracks = new Int_t[fNParticles];
804     fMultiRecTracks = new Int_t[fNParticles];
805     for (Int_t i = 0; i<fNParticles; i++) {
806       fIndexRecTracks[4*i] = -1;
807       fIndexRecTracks[4*i+1] = -1;
808       fIndexRecTracks[4*i+2] = -1;
809       fIndexRecTracks[4*i+3] = -1;
810
811       fFakeRecTracks[i] = 0;
812       fMultiRecTracks[i] = 0;
813     }
814   
815     cout<<"Start to process event "<<fEventNr<<endl;
816     cout<<"\tfNParticles = "<<fNParticles<<endl;
817     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
818     if (TreeTLoop(eventNr)>0) return 1;
819
820     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
821     if (TreeGenLoop(eventNr)>0) return 1;
822     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
823
824     delete fTreeRecTracks;
825     delete [] fIndexRecTracks;
826     delete [] fFakeRecTracks;
827     delete [] fMultiRecTracks;
828   }
829
830   CloseOutputFile();
831
832   cerr<<"Exec finished"<<endl;
833   delete gAlice;
834
835   timer.Stop();
836   timer.Print();
837   return 0;
838 }
839 ////////////////////////////////////////////////////////////////////////
840 Bool_t TPCCmpTr::ConnectGenTree()
841 {
842 //
843 // connect all branches from the genTracksTree
844 // use the same variables as for the new cmp tree, it may work
845 //
846   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
847   if (!fFileGenTracks) {
848     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
849     return kFALSE;
850   }
851   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
852   if (!fTreeGenTracks) {
853     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
854         <<fFnGenTracks<<endl;
855     return kFALSE;
856   }
857   //
858   fMCInfo = new AliTPCGenInfo;
859   fMCInfo->fReferences = new TClonesArray("AliTrackReference");  
860   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
861
862   //
863   //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
864
865
866   if (fDebug > 1) {
867     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
868   }
869   return kTRUE;
870 }
871
872
873 ////////////////////////////////////////////////////////////////////////
874 void TPCCmpTr::CreateTreeCmp() 
875 {
876   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
877   if (!fFileCmp) {
878     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
879     return;
880   }
881
882
883   fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
884   //
885   fMCInfo = new AliTPCGenInfo;
886   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
887   fRecInfo = new AliTPCRecInfo;
888   //
889   fTPCTrack = new AliTPCtrack;
890    //
891   fTreeCmp->Branch("MC","AliTPCGenInfo",&fMCInfo);
892   fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
893   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
894   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
895   //
896   fTreeCmp->AutoSave(); 
897
898 }
899 ////////////////////////////////////////////////////////////////////////
900 void TPCCmpTr::CloseOutputFile()  
901 {
902   if (!fFileCmp) {
903     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
904     return;
905   }
906   fFileCmp->cd();
907   fTreeCmp->Write();    
908   delete fTreeCmp;
909   
910   fFileCmp->Close();
911   delete fFileCmp;
912   return;
913 }
914 ////////////////////////////////////////////////////////////////////////
915
916 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
917                             AliTPCParam *paramTPC) {
918
919   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
920   Int_t index[4];
921   paramTPC->Transform0to1(x,index);
922   paramTPC->Transform1to2(x,index);
923   return TVector3(x);
924 }
925 ////////////////////////////////////////////////////////////////////////
926
927 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
928 {
929   //
930   // loop over all TPC reconstructed tracks and store info in memory
931   //
932   
933   
934   if (!fTreeRecTracks) {
935     cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
936     return 1;
937   }
938   //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
939   
940   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
941   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
942                       <<nEntries<<endl;
943   TBranch * br= fTreeRecTracks->GetBranch("tracks");
944   br->SetAddress(&fTPCTrack);
945   TBranch *brp = 0;
946   if (fTreePoints) brp = fTreePoints->GetBranch("debug");
947
948   if (fTracks){
949     fTracks->Delete();    
950     delete fTracks;
951   }
952   if (fTrackPoints){
953     fTrackPoints->Delete();
954     delete fTrackPoints;
955     fTrackPoints =0;
956   }
957   fTracks      = new TObjArray(nEntries);
958   if (brp){
959     fTrackPoints = new TObjArray(nEntries);
960   }
961   else fTrackPoints = 0;
962
963   //
964   //load tracks to the memory
965   for (Int_t i=0; i<nEntries;i++){
966     AliTPCtrack * track = new AliTPCtrack;
967     br->SetAddress(&track);
968     br->GetEntry(i);
969     fTracks->AddAt(track,i);
970   }
971   //
972   //load track points to the memory
973   if (brp) for (Int_t i=0; i<nEntries;i++){
974     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
975     brp->SetAddress(&arr);
976     brp->GetEntry(i);
977     if (arr!=0)
978       for (Int_t j=0;j<arr->GetEntriesFast();j++){
979         AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
980         if (point && point->fID>=0){
981           fTrackPoints->AddAt(arr,point->fID);
982           break;
983         }
984       }    
985   }
986   //
987
988   //
989   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
990     //br->GetEntry(iEntry);
991     fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
992     //
993     Int_t label = fTPCTrack->GetLabel();
994     Int_t absLabel = abs(label);
995     if (absLabel < fNParticles) {
996       //      fIndexRecTracks[absLabel] =  iEntry;
997       if (label < 0) fFakeRecTracks[absLabel]++;
998       
999       if (fMultiRecTracks[absLabel]>0){
1000         if (fMultiRecTracks[absLabel]<4)
1001           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
1002       }
1003       else      
1004         fIndexRecTracks[absLabel*4] =  iEntry;
1005       fMultiRecTracks[absLabel]++;
1006     }
1007   }  
1008
1009   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1010
1011   return 0;
1012 }
1013 ////////////////////////////////////////////////////////////////////////
1014 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1015 {
1016 //
1017 // loop over all entries for a given event, find corresponding 
1018 // rec. track and store in the fTreeCmp
1019 //
1020   
1021   Int_t entry = fNextTreeGenEntryToRead;
1022   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1023   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1024       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1025   TBranch * branch = fTreeCmp->GetBranch("RC");
1026   while (entry < nParticlesTR) {
1027     fTreeGenTracks->GetEntry(entry);
1028     entry++;
1029     if (fEventNr < eventNr) continue;
1030     if (fEventNr > eventNr) break;
1031     fNextTreeGenEntryToRead = entry-1;
1032     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1033       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1034     }
1035
1036     fRecInfo->Reset();
1037
1038
1039     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1040       fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1041       
1042       //      if (nBytes > 0) {
1043       if (fTPCTrack) {
1044         //
1045         TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1046         local.GetXYZ(fRecInfo->fTRLocalCoord);
1047         //
1048         // find nearest track if multifound
1049         if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
1050           Float_t dz = TMath::Abs(local.Z()-fTPCTrack->GetZ());
1051           for (Int_t i=1;i<4;i++){
1052             if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
1053               AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
1054               if  (TMath::Abs(local.Z()-track->GetZ())<dz)
1055                 fTPCTrack = track;                 
1056             }
1057           }
1058         }
1059         fRecInfo->fTP=0;
1060         if (fTrackPoints){
1061           Int_t id = fTPCTrack->GetUniqueID();
1062           if (fTrackPoints->UncheckedAt(id)){
1063             fRecInfo->fTP = (TClonesArray*)fTrackPoints->UncheckedAt(id);
1064             //      fTrackPoints->AddAt(0,id);   //not owner anymore
1065           }
1066         }
1067         fRecInfo->fTPCTrack =*fTPCTrack; 
1068         fRecInfo->fReconstructed = 1;
1069         fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1070         fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1071         fRecInfo->fdEdx = fTPCTrack->GetdEdx();
1072         //
1073         fRecInfo->fTPCTrack.PropagateTo(local.X());
1074         Double_t par[5];
1075         //
1076         Double_t localX = local.X();
1077         fTPCTrack->GetExternalParameters(localX,par);
1078         fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1079         if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*kPI;
1080         if (fRecInfo->fRecPhi>=2*kPI) fRecInfo->fRecPhi-=2*kPI;
1081 //        fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1082         fRecInfo->fLambda = TMath::ATan(par[3]);
1083         fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1084       }
1085
1086     } 
1087
1088     branch->SetAddress(&fRecInfo); // set all pointers
1089     fTreeCmp->Fill();
1090
1091   }
1092   fTreeCmp->AutoSave();
1093   fTracks->Delete();
1094   fTPCTrack =0;
1095   if (fTrackPoints){
1096     fTrackPoints->Delete();
1097     delete fTrackPoints;
1098     fTrackPoints =0;
1099   }
1100   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1101
1102   return 0;
1103 }
1104 ////////////////////////////////////////////////////////////////////////