Message commented out
[u/mrichter/AliRoot.git] / TPC / AliTPCComparisonMI.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Time Projection Chamber                                                  //
20 //  Comparison macro for TPC                                                 //
21 //  responsible: 
22 //  marian.ivanov@cern.ch                                                    //
23 //
24 //
25 // usage:
26
27 //
28 //  TO CREATE OF TREE WITH MC INFO
29
30 //  .L AliTPCComparisonMI.C+
31 //  TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",1,0)
32 //  t->Exec()
33 //  
34 //  TO CREATE COMPARISON TREE
35 //  TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",1,0);
36 //  t2->Exec()
37 //
38 //  EXAMPLE OF COMPARISON VISUALIZATION SESSION 
39 //
40 // .L AliTPCComparisonMI.C+
41 // TCut cprim("cprim","MC.fVDist[3]<1");
42 // TCut cprim("cprim","MC.fVDist[3]<1");
43 // TCut crec("crec","fReconstructed==1");
44 // TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
45 // TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
46 // TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
47 // TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
48 //
49 //
50 // AliTPCComparisonDraw comp;
51 // comp.SetIO();
52 // (comp.EffVsPt("MC.fRowsWithDigits>120"+cteta1+cprim,"1"))->Draw()
53 // (comp.EffVsRM("MC.fRowsWithDigits>20"+cteta1+cprim,"1"))->Draw()
54 // (comp.EffVsRS("MC.fRowsWithDigits>20"+cteta1+cpos1,"1"))->Draw()
55 //
56 // (comp.ResPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
57 // (comp.MeanPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
58 // (comp.ResdEdxvsN("RC.fReconstructed==1&&MC.fPrim<1.5&&abs(MC.fParticle.fPdgCode)==211&&MC.fParticle.P()>0.35"+cteta1+cpos1+cprim,"1",100,160,4))->Draw() 
59
60
61 ///////////////////////////////////////////////////////////////////////////////
62
63
64
65
66 #if !defined(__CINT__) || defined(__MAKECINT__)
67 #include <stdio.h>
68 #include <string.h>
69 //ROOT includes
70 #include "Rtypes.h"
71 #include "TFile.h"
72 #include "TTree.h"
73 #include "TChain.h"
74 #include "TCut.h"
75 #include "TString.h"
76 #include "TBenchmark.h"
77 #include "TStopwatch.h"
78 #include "TParticle.h"
79 #include "TSystem.h"
80 #include "TTimer.h"
81 #include "TVector3.h"
82 #include "TPad.h"
83 #include "TCanvas.h"
84 #include "TH1F.h"
85 #include "TH2F.h"
86 #include "TF1.h"
87 #include "TText.h"
88 #include "Getline.h"
89 #include "TStyle.h"
90
91 //ALIROOT includes
92 #include "AliRun.h"
93 #include "AliStack.h"
94 #include "AliTPCtrack.h"
95 #include "AliSimDigits.h"
96 #include "AliTPCParam.h"
97 #include "AliTPC.h"
98 #include "AliTPCLoader.h"
99 #include "AliDetector.h"
100 #include "AliTrackReference.h"
101 #include "AliRun.h"
102 #include "AliTPCParamSR.h"
103 #include "AliTracker.h"
104 #include "AliComplexCluster.h"
105 #include "AliMagF.h"
106 #endif
107 #include "AliTPCComparisonMI.h"
108
109
110
111
112
113
114
115  
116 //
117 //
118
119 AliTPCParam * GetTPCParam(){
120   AliTPCParamSR * par = new AliTPCParamSR;
121   par->Update();
122   return par;
123 }
124
125
126 //_____________________________________________________________________________
127 Float_t TPCBetheBloch(Float_t bg)
128 {
129   //
130   // Bethe-Bloch energy loss formula
131   //
132   const Double_t kp1=0.76176e-1;
133   const Double_t kp2=10.632;
134   const Double_t kp3=0.13279e-4;
135   const Double_t kp4=1.8631;
136   const Double_t kp5=1.9479;
137
138   Double_t dbg = (Double_t) bg;
139
140   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
141
142   Double_t aa = TMath::Power(beta,kp4);
143   Double_t bb = TMath::Power(1./dbg,kp5);
144
145   bb=TMath::Log(kp3+bb);
146   
147   return ((Float_t)((kp2-aa-bb)*kp1/aa));
148 }
149
150
151
152
153
154 ////////////////////////////////////////////////////////////////////////
155 AliTPCGenInfo::AliTPCGenInfo()
156 {
157   fReferences  = 0;
158   fTRdecay.SetTrack(-1);
159 }
160
161 AliTPCGenInfo::~AliTPCGenInfo()
162 {
163   if (fReferences) delete fReferences;
164   fReferences =0;
165   
166 }
167
168 void AliTPCGenInfo::Update()
169 {
170   //
171   //
172   fMCtracks =1;
173   if (!fReferences) return;
174   Float_t direction=1;
175   //Float_t rlast=0;
176   for (Int_t iref =0;iref<fReferences->GetEntriesFast();iref++){
177     AliTrackReference * ref = (AliTrackReference *) fReferences->At(iref);
178     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
179     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
180     if (iref==0) direction = newdirection;
181     if ( newdirection*direction<0){
182       //changed direction
183       direction = newdirection;
184       fMCtracks+=1;
185     }
186     //rlast=r;                      
187   }
188 }
189
190 ////////////////////////////////////////////////////////////////////////
191
192 ////////////////////////////////////////////////////////////////////////
193 //
194 // End of implementation of the class AliTPCGenInfo
195 //
196 ////////////////////////////////////////////////////////////////////////
197
198
199
200 ////////////////////////////////////////////////////////////////////////
201 digitRow::digitRow()
202 {
203   Reset();
204 }
205 ////////////////////////////////////////////////////////////////////////
206 digitRow & digitRow::operator=(const digitRow &digOld)
207 {
208   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
209   return (*this);
210 }
211 ////////////////////////////////////////////////////////////////////////
212 void digitRow::SetRow(Int_t row) 
213 {
214   if (row >= 8*kgRowBytes) {
215     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
216     return;
217   }
218   Int_t iC = row/8;
219   Int_t iB = row%8;
220   SETBIT(fDig[iC],iB);
221 }
222
223 ////////////////////////////////////////////////////////////////////////
224 Bool_t digitRow::TestRow(Int_t row)
225 {
226 //
227 // return kTRUE if row is on
228 //
229   Int_t iC = row/8;
230   Int_t iB = row%8;
231   return TESTBIT(fDig[iC],iB);
232 }
233 ////////////////////////////////////////////////////////////////////////
234 Int_t digitRow::RowsOn(Int_t upto)
235 {
236 //
237 // returns number of rows with a digit  
238 // count only rows less equal row number upto
239 //
240   Int_t total = 0;
241   for (Int_t i = 0; i<kgRowBytes; i++) {
242     for (Int_t j = 0; j < 8; j++) {
243       if (i*8+j > upto) return total;
244       if (TESTBIT(fDig[i],j))  total++;
245     }
246   }
247   return total;
248 }
249 ////////////////////////////////////////////////////////////////////////
250 void digitRow::Reset()
251 {
252 //
253 // resets all rows to zero
254 //
255   for (Int_t i = 0; i<kgRowBytes; i++) {
256     fDig[i] <<= 8;
257   }
258 }
259 ////////////////////////////////////////////////////////////////////////
260 Int_t digitRow::Last()
261 {
262 //
263 // returns the last row number with a digit
264 // returns -1 if now digits 
265 //
266   for (Int_t i = kgRowBytes-1; i>=0; i--) {
267     for (Int_t j = 7; j >= 0; j--) {
268       if TESTBIT(fDig[i],j) return i*8+j;
269     }
270   }
271   return -1;
272 }
273 ////////////////////////////////////////////////////////////////////////
274 Int_t digitRow::First()
275 {
276 //
277 // returns the first row number with a digit
278 // returns -1 if now digits 
279 //
280   for (Int_t i = 0; i<kgRowBytes; i++) {
281     for (Int_t j = 0; j < 8; j++) {
282       if (TESTBIT(fDig[i],j)) return i*8+j;
283     }
284   }
285   return -1;
286 }
287
288 ////////////////////////////////////////////////////////////////////////
289 //
290 // end of implementation of a class digitRow
291 //
292 ////////////////////////////////////////////////////////////////////////
293   
294 ////////////////////////////////////////////////////////////////////////
295 TPCFindGenTracks::TPCFindGenTracks()
296 {
297   fMCInfo = new AliTPCGenInfo;
298   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
299
300   Reset();
301 }
302
303 ////////////////////////////////////////////////////////////////////////
304 TPCFindGenTracks::TPCFindGenTracks(const char * fnGalice, const char* fnRes,
305                                    Int_t nEvents, Int_t firstEvent)
306 {
307   Reset();
308   fFirstEventNr = firstEvent;
309   fEventNr = firstEvent;
310   fNEvents = nEvents;
311   //   fFnRes = fnRes;
312   sprintf(fFnRes,"%s",fnRes);
313   //
314   fLoader = AliRunLoader::Open(fnGalice);
315   if (gAlice){
316     delete gAlice->GetRunLoader();
317     delete gAlice;
318     gAlice = 0x0;
319   }
320   if (fLoader->LoadgAlice()){
321     cerr<<"Error occured while l"<<endl;
322   }
323   Int_t nall = fLoader->GetNumberOfEvents();
324   if (nEvents==0) {
325     nEvents =nall;
326     fNEvents=nall;
327     fFirstEventNr=0;
328   }    
329
330   if (nall<=0){
331     cerr<<"no events available"<<endl;
332     fEventNr = 0;
333     return;
334   }
335   if (firstEvent+nEvents>nall) {
336     fEventNr = nall-firstEvent;
337     cerr<<"restricted number of events availaible"<<endl;
338   }
339   AliMagF * magf = gAlice->Field();
340   AliTracker::SetFieldMap(magf);
341 }
342
343 ////////////////////////////////////////////////////////////////////////
344 void TPCFindGenTracks::Reset()
345 {
346   fEventNr = 0;
347   fNEvents = 0;
348   fTreeGenTracks = 0;
349   fFileGenTracks = 0;
350   fContainerDigitRow = 0;
351   //
352   fReferences = 0;
353   fReferenceIndex0 = 0;
354   fReferenceIndex1 = 0;
355   fDecayRef   = 0;
356   //
357   fDebug = 0;
358   fVPrim[0] = -1000.;
359   fVPrim[1] = -1000.;
360   fVPrim[2] = -1000.;
361   fParamTPC = 0;
362 }
363 ////////////////////////////////////////////////////////////////////////
364 TPCFindGenTracks::~TPCFindGenTracks()
365 {
366   
367   if (fLoader){
368     fLoader->UnloadgAlice();
369     gAlice = 0;
370     delete fLoader;
371   }
372 }
373
374
375 Int_t  TPCFindGenTracks::SetIO()
376 {
377   //
378   // 
379   CreateTreeGenTracks();
380   if (!fTreeGenTracks) return 1;
381   //  AliTracker::SetFieldFactor(); 
382  
383   fParamTPC = GetTPCParam();
384   //
385   return 0;
386 }
387
388 ////////////////////////////////////////////////////////////////////////
389 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
390 {
391   //
392   // 
393   // SET INPUT
394   fLoader->SetEventNumber(eventNr);
395   //
396   fLoader->LoadHeader();
397   fLoader->LoadKinematics();  
398   fStack = fLoader->Stack();
399   //
400   fLoader->LoadTrackRefs();
401   fTreeTR = fLoader->TreeTR();
402   //
403   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
404   tpcl->LoadDigits();
405   fTreeD = tpcl->TreeD();  
406   return 0;
407 }
408
409 Int_t TPCFindGenTracks::CloseIOEvent()
410 {
411   fLoader->UnloadHeader();
412   fLoader->UnloadKinematics();
413   fLoader->UnloadTrackRefs();
414   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
415   tpcl->UnloadDigits();
416   return 0;
417 }
418
419 Int_t TPCFindGenTracks::CloseIO()
420 {
421   fLoader->UnloadgAlice();
422   return 0;
423 }
424
425
426
427 ////////////////////////////////////////////////////////////////////////
428 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
429 {
430   fNEvents = nEvents;
431   fFirstEventNr = firstEventNr;
432   return Exec();
433 }
434
435 ////////////////////////////////////////////////////////////////////////
436 Int_t TPCFindGenTracks::Exec()  
437 {
438   TStopwatch timer;
439   timer.Start();
440   Int_t status =SetIO();
441   if (status>0) return status;
442   //
443
444   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
445        fEventNr++) {
446     SetIO(fEventNr);
447     fNParticles = fStack->GetNtrack();
448     fContainerDigitRow = new digitRow[fNParticles];
449     //
450     fReferences      = new AliTrackReference[fgMaxTR];
451     fReferenceIndex0 = new Int_t[fNParticles];
452     fReferenceIndex1 = new Int_t[fNParticles];
453     fDecayRef        = new AliTrackReference[fNParticles];
454
455     for (Int_t i = 0; i<fNParticles; i++) {
456       fReferenceIndex0[i] = -1;
457       fReferenceIndex1[i] = -1;
458     }
459     
460     cout<<"Start to process event "<<fEventNr<<endl;
461     cout<<"\tfNParticles = "<<fNParticles<<endl;
462     if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
463     if (TreeDLoop()>0) return 1;
464     if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
465     if (TreeTRLoop()>0) return 1;
466     if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
467     if (TreeKLoop()>0) return 1;
468     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
469
470     delete [] fContainerDigitRow;
471     delete [] fReferences;
472     delete [] fReferenceIndex0;
473     delete [] fReferenceIndex1;    
474     delete [] fDecayRef;
475     CloseIOEvent();
476   }
477   //
478   CloseIO();
479   CloseOutputFile();
480
481   cerr<<"Exec finished"<<endl;
482
483   timer.Stop();
484   timer.Print();
485   return 0;
486 }
487 ////////////////////////////////////////////////////////////////////////
488 void TPCFindGenTracks::CreateTreeGenTracks() 
489 {
490   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
491   if (!fFileGenTracks) {
492     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
493     return;
494   }
495   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
496   
497
498
499   fMCInfo = new AliTPCGenInfo;
500   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
501
502   fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
503
504   fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
505
506   fTreeGenTracks->AutoSave();
507 }
508 ////////////////////////////////////////////////////////////////////////
509 void TPCFindGenTracks::CloseOutputFile() 
510 {
511   if (!fFileGenTracks) {
512     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
513     return;
514   }
515   fFileGenTracks->cd();
516   fTreeGenTracks->Write();  
517   delete fTreeGenTracks;
518   fFileGenTracks->Close();
519   delete fFileGenTracks;
520   return;
521 }
522
523 ////////////////////////////////////////////////////////////////////////
524 Int_t TPCFindGenTracks::TreeKLoop()
525 {
526 //
527 // open the file with treeK
528 // loop over all entries there and save information about some tracks
529 //
530
531   AliStack * stack = fStack;
532   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
533   
534   if (fDebug > 0) {
535     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
536         <<fEventNr<<endl;
537   }
538   
539   Int_t  ipdg = 0;
540   TParticlePDG *ppdg = 0;
541
542  
543   // not all generators give primary vertex position. Take the vertex
544   // of the particle 0 as primary vertex.
545   TDatabasePDG  pdg; //get pdg table
546   
547   //thank you very much root for this
548   TBranch * br = fTreeGenTracks->GetBranch("MC");
549
550   TParticle *particle = stack->ParticleFromTreeK(0);
551   fVPrim[0] = particle->Vx();
552   fVPrim[1] = particle->Vy();
553   fVPrim[2] = particle->Vz();
554   for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
555
556 // load only particles with TR
557     if (fReferenceIndex0[iParticle]<0) continue;
558     //////////////////////////////////////////////////////////////////////
559     fMCInfo->fLabel = iParticle;
560
561     fMCInfo->fRow = (fContainerDigitRow[iParticle]);
562     fMCInfo->fRowsWithDigits = fMCInfo->fRow.RowsOn();    
563     if (fMCInfo->fRowsWithDigits<10) continue;
564     fMCInfo->fRowsWithDigitsInn = fMCInfo->fRow.RowsOn(63); // 63 = number of inner rows
565     fMCInfo->fRowsTrackLength = fMCInfo->fRow.Last() - fMCInfo->fRow.First();
566     fMCInfo->fDigitsInSeed = 0;
567     if (fMCInfo->fRow.TestRow(seedRow11) && fMCInfo->fRow.TestRow(seedRow12)) 
568       fMCInfo->fDigitsInSeed = 1;
569     if (fMCInfo->fRow.TestRow(seedRow21) && fMCInfo->fRow.TestRow(seedRow22)) 
570       fMCInfo->fDigitsInSeed += 10;
571     //
572     //
573     //
574     fMCInfo->fParticle = *(stack->Particle(iParticle));
575     //
576     //
577     fMCInfo->fLabel = iParticle;
578     fMCInfo->fVDist[0] = fMCInfo->fParticle.Vx()-fVPrim[0];
579     fMCInfo->fVDist[1] = fMCInfo->fParticle.Vy()-fVPrim[1];
580     fMCInfo->fVDist[2] = fMCInfo->fParticle.Vz()-fVPrim[2]; 
581     fMCInfo->fVDist[3] = TMath::Sqrt(fMCInfo->fVDist[0]*fMCInfo->fVDist[0]+
582                                       fMCInfo->fVDist[1]*fMCInfo->fVDist[1]+fMCInfo->fVDist[2]*fMCInfo->fVDist[2]);
583
584     //
585     Int_t index = fReferenceIndex0[iParticle];
586     AliTrackReference  ref  = fReferences[index];
587     // if (ref.GetTrack()!=iParticle)
588     //  printf("problem2\n");
589     fMCInfo->fTrackRef = ref;
590   
591     Int_t rfindex =0;
592     if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
593     fMCInfo->fReferences = 
594       new TClonesArray("AliTrackReference");
595     fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
596     //
597     for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
598       AliTrackReference  ref  = fReferences[i];
599       AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
600       if (ref.GetTrack()!=iParticle){
601         //printf("problem5\n"); 
602         continue;
603       }
604       *ref2 = ref;
605       rfindex++;
606     }   
607     //
608     //
609     ipdg = fMCInfo->fParticle.GetPdgCode();
610     fMCInfo->fPdg = ipdg;
611     ppdg = pdg.GetParticle(ipdg);
612     // calculate primary ionization per cm
613     if (ppdg){
614       Float_t mass = ppdg->Mass();
615       Float_t p = TMath::Sqrt(fMCInfo->fTrackRef.Px()*fMCInfo->fTrackRef.Px()+
616                               fMCInfo->fTrackRef.Py()*fMCInfo->fTrackRef.Py()+
617                               fMCInfo->fTrackRef.Pz()*fMCInfo->fTrackRef.Pz());
618       
619       //      Float_t betagama = fMCInfo->fParticle.P()/mass;
620       Float_t betagama = p /mass;
621       fMCInfo->fPrim = TPCBetheBloch(betagama);
622     }   
623     fMCInfo->fTPCdecay=kFALSE;
624     if (fDecayRef[iParticle].GetTrack()>0){
625       fMCInfo->fTRdecay  = fDecayRef[iParticle];
626       fMCInfo->fDecayCoord[0] = fMCInfo->fTRdecay.X();
627       fMCInfo->fDecayCoord[1] = fMCInfo->fTRdecay.Y();
628       fMCInfo->fDecayCoord[2] = fMCInfo->fTRdecay.Z();
629       if ( (fMCInfo->fTRdecay.R()<250)&&(fMCInfo->fTRdecay.R()>85) && (TMath::Abs(fMCInfo->fTRdecay.Z())<250) ){
630         fMCInfo->fTPCdecay=kTRUE;
631       }
632     }
633     else{
634       fMCInfo->fTRdecay.SetTrack(-1);
635       fMCInfo->fDecayCoord[0] = 0;
636       fMCInfo->fDecayCoord[1] = 0;
637       fMCInfo->fDecayCoord[2] = 0;
638     }
639     fMCInfo->Update();
640     //////////////////////////////////////////////////////////////////////
641     
642     br->SetAddress(&fMCInfo);
643     fTreeGenTracks->Fill();
644
645   }
646   fTreeGenTracks->AutoSave();
647
648   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
649
650   return 0;
651 }
652
653
654
655
656 ////////////////////////////////////////////////////////////////////////
657 Int_t TPCFindGenTracks::TreeDLoop()
658 {
659 //
660 // open the file with treeD
661 // loop over all entries there and save information about some tracks
662 //
663
664   Int_t nInnerSector = fParamTPC->GetNInnerSector();
665   Int_t rowShift = 0;
666   //  Int_t zero=fParamTPC->GetZeroSup();
667   Int_t zero=fParamTPC->GetZeroSup()+6;
668
669   
670   //char treeDName[100]; 
671   //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
672   //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
673   //
674   //
675   AliSimDigits digitsAddress, *digits=&digitsAddress;
676   fTreeD->GetBranch("Segment")->SetAddress(&digits);
677
678   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
679   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
680   for (Int_t i=0; i<sectorsByRows; i++) {
681 //  for (Int_t i=5720; i<sectorsByRows; i++) {
682     if (!fTreeD->GetEvent(i)) continue;
683     Int_t sec,row;
684     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
685     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
686 //    cerr<<sec<<' '<<row<<endl;
687
688 // here I expect that upper sectors follow lower sectors
689     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
690     //
691     //digits->ExpandTrackBuffer();     
692     //Int_t *tracks = digits->GetTracks();
693     digits->First();    
694     do {
695       Int_t iRow=digits->CurrentRow();
696       Int_t iColumn=digits->CurrentColumn();
697       Short_t digitValue = digits->CurrentDigit();
698 //      cout<<"Inner loop: sector, iRow, iColumn "
699 //        <<sec<<" "<<iRow<<" "<<iColumn<<endl;
700       if (digitValue >= zero) {
701         Int_t label;
702         for (Int_t j = 0; j<3; j++) {
703           label = digits->GetTrackID(iRow,iColumn,j); 
704           //label = digits->GetTrackIDFast(iRow,iColumn,j); 
705           
706           if (label >= fNParticles) { //don't label from bakground event
707             continue;
708           }
709           if (label >= 0 && label <= fNParticles) {
710 //        if (label >= 0 && label <= fDebug) {
711             if (fDebug > 6 ) {
712               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
713                   <<sec<<" "
714                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
715                   <<" "<<row<<endl;
716             }   
717             fContainerDigitRow[label].SetRow(row+rowShift);
718           }
719         }
720       }
721     } while (digits->Next());
722   }
723
724   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
725
726   return 0;
727 }
728
729
730 ////////////////////////////////////////////////////////////////////////
731 Int_t TPCFindGenTracks::TreeTRLoop()
732 {
733   //
734   // loop over TrackReferences and store the first one for each track
735   //
736   
737   TTree * treeTR = fTreeTR;
738   if (!treeTR) {
739     cerr<<"TreeTR not found"<<endl;
740     return 1;
741   }
742   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
743   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
744   //
745   //
746   //track references for TPC
747   TBranch *TPCBranchTR  = treeTR->GetBranch("TPC");
748   if (!TPCBranchTR) {
749     cerr<<"TPC branch in TR not found"<<endl;
750     return 1;
751   }
752   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
753   TPCBranchTR->SetAddress(&TPCArrayTR);
754   //get decay point if exist
755   TBranch *runbranch  = treeTR->GetBranch("AliRun");
756   if (!runbranch) {
757     cerr<<"Run branch in TR not found"<<endl;
758     return 1;
759   }
760   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
761   runbranch->SetAddress(&RunArrayTR);
762   //
763   //
764   //
765   Int_t index     =  0;
766   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
767     TPCBranchTR->GetEntry(iPrimPart);
768     Float_t ptstart = 0;
769     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
770       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
771       //
772       Int_t label = trackRef->GetTrack();
773       if (label<0  || label > fNParticles) {
774         continue;
775       }
776       if (fReferenceIndex0[label]<0) ptstart = trackRef->Pt();  //store pt at the TPC entrance
777       if (ptstart<fgPtCut) continue;
778
779       if (index>=fgMaxTR) continue;     //restricted size of buffer for TR
780       fReferences[index] = *trackRef;
781       fReferenceIndex1[label] = index;  // the last ref with given label
782       if (fReferenceIndex0[label]==-1) fReferenceIndex0[label] = index;   //the first index with label
783       index++;           
784     }
785     // get dacay position
786     runbranch->GetEntry(iPrimPart);    
787     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
788       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
789       //
790       if (trackRef->Pt() < fgPtCut) continue;      
791       Int_t label = trackRef->GetTrack();
792       if (label<0  || label > fNParticles) {
793         continue;
794       }
795       if (trackRef->R()>450) continue;   //not decay  in TPC
796       if (trackRef->Z()>450) continue;   //not decay  in TPC
797       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
798       
799       if (label>=fgMaxTR) continue;     //restricted size of buffer for TR
800       fDecayRef[label] = *trackRef;      
801     }
802
803   }
804   TPCArrayTR->Delete();
805   delete  TPCArrayTR;
806   RunArrayTR->Delete();
807   delete  RunArrayTR;
808
809   return 0;
810 }
811
812 ////////////////////////////////////////////////////////////////////////
813 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
814                                     AliTPCParam *paramTPC) {
815
816   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
817   Int_t index[4];
818   paramTPC->Transform0to1(x,index);
819   paramTPC->Transform1to2(x,index);
820   return x[0];
821 }
822 ////////////////////////////////////////////////////////////////////////
823
824
825
826   
827 ////////////////////////////////////////////////////////////////////////
828 TPCCmpTr::TPCCmpTr()
829 {
830   Reset();
831 }
832
833 ////////////////////////////////////////////////////////////////////////
834 TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
835                    const char* fnCmp,
836                    const char* fnGalice,
837                    Int_t nEvents, Int_t firstEvent)
838 {
839   Reset();
840   //  fFnGenTracks = fnGenTracks;
841   //  fFnCmp = fnCmp;
842   sprintf(fFnGenTracks,"%s",fnGenTracks);
843   sprintf(fFnCmp,"%s",fnCmp);
844
845   fFirstEventNr = firstEvent;
846   fEventNr = firstEvent;
847   fNEvents = nEvents;
848   
849   //
850   fLoader = AliRunLoader::Open(fnGalice);
851   if (gAlice){
852     //delete gAlice->GetRunLoader();
853     delete gAlice;
854     gAlice = 0x0;
855   }
856   if (fLoader->LoadgAlice()){
857     cerr<<"Error occured while l"<<endl;
858   }
859   Int_t nall = fLoader->GetNumberOfEvents();
860   if (nEvents==0) {
861     nEvents =nall;
862     fNEvents=nall;
863     fFirstEventNr=0;
864   }    
865
866   if (nall<=0){
867     cerr<<"no events available"<<endl;
868     fEventNr = 0;
869     return;
870   }
871   if (firstEvent+nEvents>nall) {
872     fEventNr = nall-firstEvent;
873     cerr<<"restricted number of events availaible"<<endl;
874   }
875   AliMagF * magf = gAlice->Field();
876   AliTracker::SetFieldMap(magf);
877
878 }
879
880
881 ////////////////////////////////////////////////////////////////////////
882 TPCCmpTr::~TPCCmpTr()
883 {
884   if (fLoader) {
885     delete fLoader;
886   }
887 }
888
889 //////////////////////////////////////////////////////////////
890 Int_t TPCCmpTr::SetIO()
891 {
892   //
893   // 
894   CreateTreeCmp();
895   if (!fTreeCmp) return 1;
896   fParamTPC = GetTPCParam();
897   //
898   if (!ConnectGenTree()) {
899     cerr<<"Cannot connect tree with generated tracks"<<endl;
900     return 1;
901   }
902   return 0;
903 }
904
905 //////////////////////////////////////////////////////////////
906
907 Int_t TPCCmpTr::SetIO(Int_t eventNr)
908 {
909   //
910   // 
911   // SET INPUT
912   //gAlice->GetEvent(eventNr);
913   fLoader->SetEventNumber(eventNr);  
914   //
915   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
916   tpcl->LoadTracks();
917   fTreeRecTracks = tpcl->TreeT();
918   //
919   return 0;
920 }
921
922
923
924 ////////////////////////////////////////////////////////////////////////
925 void TPCCmpTr::Reset()
926 {
927   fEventNr = 0;
928   fNEvents = 0;
929   fTreeCmp = 0;
930   //  fFnCmp = "cmpTracks.root";
931   fFileGenTracks = 0;
932   fDebug = 0;
933   //
934   fParamTPC = 0;
935   fTreeRecTracks = 0;
936   fTreePoints =0;
937
938   fTPCTrack = 0; 
939   fTracks   = 0;
940   fTrackPoints =0;
941 }
942
943 ////////////////////////////////////////////////////////////////////////
944 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
945 {
946   fNEvents = nEvents;
947   fFirstEventNr = firstEventNr;
948   return Exec();
949 }
950
951 ////////////////////////////////////////////////////////////////////////
952 Int_t TPCCmpTr::Exec()
953 {
954   TStopwatch timer;
955   timer.Start();
956
957   if (SetIO()==1) 
958     return 1;
959    
960   fNextTreeGenEntryToRead = 0;
961   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
962   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
963        eventNr++) {
964     SetIO(fEventNr);
965     fNParticles = gAlice->GetEvent(fEventNr);    
966
967     fIndexRecTracks = new Int_t[fNParticles*4];  //write at maximum 4 tracks corresponding to particle
968     fFakeRecTracks = new Int_t[fNParticles];
969     fMultiRecTracks = new Int_t[fNParticles];
970     for (Int_t i = 0; i<fNParticles; i++) {
971       fIndexRecTracks[4*i] = -1;
972       fIndexRecTracks[4*i+1] = -1;
973       fIndexRecTracks[4*i+2] = -1;
974       fIndexRecTracks[4*i+3] = -1;
975
976       fFakeRecTracks[i] = 0;
977       fMultiRecTracks[i] = 0;
978     }
979   
980     cout<<"Start to process event "<<fEventNr<<endl;
981     cout<<"\tfNParticles = "<<fNParticles<<endl;
982     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
983     if (TreeTLoop(eventNr)>0) return 1;
984
985     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
986     if (TreeGenLoop(eventNr)>0) return 1;
987     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
988
989     delete fTreeRecTracks;
990     delete [] fIndexRecTracks;
991     delete [] fFakeRecTracks;
992     delete [] fMultiRecTracks;
993   }
994
995   CloseOutputFile();
996
997   cerr<<"Exec finished"<<endl;
998   timer.Stop();
999   timer.Print();
1000   return 0;
1001
1002 }
1003 ////////////////////////////////////////////////////////////////////////
1004 Bool_t TPCCmpTr::ConnectGenTree()
1005 {
1006 //
1007 // connect all branches from the genTracksTree
1008 // use the same variables as for the new cmp tree, it may work
1009 //
1010   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1011   if (!fFileGenTracks) {
1012     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1013     return kFALSE;
1014   }
1015   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1016   if (!fTreeGenTracks) {
1017     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1018         <<fFnGenTracks<<endl;
1019     return kFALSE;
1020   }
1021   //
1022   fMCInfo = new AliTPCGenInfo;
1023   fMCInfo->fReferences = new TClonesArray("AliTrackReference");  
1024   fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1025
1026   //
1027   //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1028
1029
1030   if (fDebug > 1) {
1031     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1032   }
1033   return kTRUE;
1034 }
1035
1036
1037 ////////////////////////////////////////////////////////////////////////
1038 void TPCCmpTr::CreateTreeCmp() 
1039 {
1040   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1041   if (!fFileCmp) {
1042     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1043     return;
1044   }
1045
1046
1047   fTreeCmp    = new TTree("TPCcmpTracks","TPCcmpTracks");
1048   //
1049   fMCInfo = new AliTPCGenInfo;
1050   fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1051   fRecInfo = new AliTPCRecInfo;
1052   //
1053   fTPCTrack = new AliTPCtrack;
1054    //
1055   fTreeCmp->Branch("MC","AliTPCGenInfo",&fMCInfo);
1056   fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
1057   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
1058   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1059   //
1060   fTreeCmp->AutoSave(); 
1061
1062 }
1063 ////////////////////////////////////////////////////////////////////////
1064 void TPCCmpTr::CloseOutputFile()  
1065 {
1066   if (!fFileCmp) {
1067     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1068     return;
1069   }
1070   fFileCmp->cd();
1071   fTreeCmp->Write();    
1072   delete fTreeCmp;
1073   
1074   fFileCmp->Close();
1075   delete fFileCmp;
1076   return;
1077 }
1078 ////////////////////////////////////////////////////////////////////////
1079
1080 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
1081                             AliTPCParam *paramTPC) {
1082
1083   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1084   Int_t index[4];
1085   paramTPC->Transform0to1(x,index);
1086   paramTPC->Transform1to2(x,index);
1087   return TVector3(x);
1088 }
1089 ////////////////////////////////////////////////////////////////////////
1090
1091 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1092 {
1093   //
1094   // loop over all TPC reconstructed tracks and store info in memory
1095   //
1096   TStopwatch  timer;
1097   timer.Start();
1098   
1099   if (!fTreeRecTracks) {
1100     cerr<<"Can't get a tree with TPC rec. tracks  "<<endl;
1101     return 1;
1102   }
1103   //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
1104   
1105   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1106   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1107                       <<nEntries<<endl;
1108   TBranch * br= fTreeRecTracks->GetBranch("tracks");
1109   br->SetAddress(&fTPCTrack);
1110   TBranch *brp = 0;
1111   if (fTreePoints) brp = fTreePoints->GetBranch("debug");
1112
1113   if (fTracks){
1114     fTracks->Delete();    
1115     delete fTracks;
1116   }
1117   if (fTrackPoints){
1118     fTrackPoints->Delete();
1119     delete fTrackPoints;
1120     fTrackPoints =0;
1121   }
1122   fTracks      = new TObjArray(nEntries);
1123   if (brp){
1124     fTrackPoints = new TObjArray(nEntries);
1125   }
1126   else fTrackPoints = 0;
1127
1128   //
1129   //load tracks to the memory
1130   for (Int_t i=0; i<nEntries;i++){
1131     AliTPCtrack * track = new AliTPCtrack;
1132     br->SetAddress(&track);
1133     br->GetEntry(i);
1134     fTracks->AddAt(track,i);
1135   }
1136   //
1137   //load track points to the memory
1138   if (brp) for (Int_t i=0; i<nEntries;i++){
1139     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
1140     brp->SetAddress(&arr);
1141     brp->GetEntry(i);
1142     if (arr!=0)
1143       for (Int_t j=0;j<arr->GetEntriesFast();j++){
1144         AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
1145         if (point && point->fID>=0){
1146           fTrackPoints->AddAt(arr,point->fID);
1147           break;
1148         }
1149       }    
1150   }
1151   //
1152
1153   //
1154   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1155     //br->GetEntry(iEntry);
1156     fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
1157     //
1158     Int_t label = fTPCTrack->GetLabel();
1159     Int_t absLabel = abs(label);
1160     if (absLabel < fNParticles) {
1161       //      fIndexRecTracks[absLabel] =  iEntry;
1162       if (label < 0) fFakeRecTracks[absLabel]++;
1163       
1164       if (fMultiRecTracks[absLabel]>0){
1165         if (fMultiRecTracks[absLabel]<4)
1166           fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] =  iEntry;      
1167       }
1168       else      
1169         fIndexRecTracks[absLabel*4] =  iEntry;
1170       fMultiRecTracks[absLabel]++;
1171     }
1172   }  
1173   printf("Time spended in TreeTLoop\n");
1174   timer.Print();
1175   
1176   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1177
1178   return 0;
1179 }
1180 ////////////////////////////////////////////////////////////////////////
1181 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1182 {
1183 //
1184 // loop over all entries for a given event, find corresponding 
1185 // rec. track and store in the fTreeCmp
1186 //
1187   TStopwatch timer;
1188   timer.Start();
1189   Int_t entry = fNextTreeGenEntryToRead;
1190   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1191   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1192       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1193   TBranch * branch = fTreeCmp->GetBranch("RC");
1194   branch->SetAddress(&fRecInfo); // set all pointers
1195
1196   while (entry < nParticlesTR) {
1197     fTreeGenTracks->GetEntry(entry);
1198     entry++;
1199     if (fEventNr < eventNr) continue;
1200     if (fEventNr > eventNr) break;
1201     fNextTreeGenEntryToRead = entry-1;
1202     if (fDebug > 2 && fMCInfo->fLabel < 10) {
1203       cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1204     }
1205
1206     fRecInfo->Reset();
1207
1208
1209     if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1210       fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1211       
1212       //      if (nBytes > 0) {
1213       if (fTPCTrack) {
1214         //
1215         TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1216         local.GetXYZ(fRecInfo->fTRLocalCoord);
1217         //
1218         // find nearest track if multifound
1219         if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
1220           Float_t dz = TMath::Abs(local.Z()-fTPCTrack->GetZ());
1221           for (Int_t i=1;i<4;i++){
1222             if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
1223               AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
1224               if  (TMath::Abs(local.Z()-track->GetZ())<dz)
1225                 fTPCTrack = track;                 
1226             }
1227           }
1228         }
1229         fRecInfo->fTP=0;
1230         if (fTrackPoints){
1231           Int_t id = fTPCTrack->GetUniqueID();
1232           if (fTrackPoints->UncheckedAt(id)){
1233             fRecInfo->fTP = (TClonesArray*)fTrackPoints->UncheckedAt(id);
1234             //      fTrackPoints->AddAt(0,id);   //not owner anymore
1235           }
1236         }
1237         fRecInfo->fTPCTrack =*fTPCTrack; 
1238         fRecInfo->fReconstructed = 1;
1239         fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1240         fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1241         fRecInfo->fdEdx = fTPCTrack->GetdEdx();
1242         //
1243         fRecInfo->fTPCTrack.PropagateTo(local.X());
1244         Double_t par[5];
1245         //
1246         Double_t localX = local.X();
1247         fTPCTrack->GetExternalParameters(localX,par);
1248         fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1249         if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*TMath::Pi();
1250         if (fRecInfo->fRecPhi>=2*TMath::Pi()) fRecInfo->fRecPhi-=2*TMath::Pi();
1251 //        fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1252         fRecInfo->fLambda = TMath::ATan(par[3]);
1253         fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1254       }
1255
1256     } 
1257
1258     fTreeCmp->Fill();
1259   }
1260   fTreeCmp->AutoSave();
1261   fTracks->Delete();
1262   fTPCTrack =0;
1263   if (fTrackPoints){
1264     fTrackPoints->Delete();
1265     delete fTrackPoints;
1266     fTrackPoints =0;
1267   } 
1268   printf("Time spended in TreeKLoop\n");
1269   timer.Print();
1270   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1271
1272   return 0;
1273 }
1274 ////////////////////////////////////////////////////////////////////////
1275 ////////////////////////////////////////////////////////////////////////
1276
1277 void AliTPCComparisonDraw::SetIO(const char *fname)
1278 {
1279   //
1280    TFile* file = TFile::Open(fname);
1281    if (!file) {
1282      printf("Could not open file  - generated new one\n"); 
1283      TFile* filegen = TFile::Open("genTracks.root");
1284      if (!filegen){
1285        printf("FILE with MC information is generated\n"); 
1286        TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",0);
1287        t->Exec();
1288        delete t;
1289      }
1290      filegen = TFile::Open("genTracks.root");
1291      if (!filegen){
1292        printf("COMPARISON  FILE COULDNT BE GENERATED \n");
1293        return;
1294      }
1295      printf("COMPARISON  FILE IS GENERATED \n");
1296      TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",0);
1297      t2->Exec();
1298    }
1299    file = TFile::Open(fname);
1300    if (!file){
1301      printf("Comparison file couldn't be generated\n"); 
1302      return;
1303    }
1304    //
1305    fTree = (TTree*) file->Get("TPCcmpTracks");
1306    if (!fTree) {
1307     printf("no track comparison tree found\n");
1308     file->Close();
1309     delete file;
1310   }
1311 }
1312
1313 void AliTPCComparisonDraw::ResPt()
1314 {
1315  //
1316   //
1317   gStyle->SetOptFit();
1318   TCanvas *c1=new TCanvas("TPC pt resolution","TPC pt resolution",0,0,700,850);
1319   c1->Draw(); c1->cd();
1320   TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99); 
1321   p1->Draw();p1->cd(); 
1322   p1->SetGridx(); p1->SetGridy();
1323   //
1324   c1->cd();
1325   TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49); 
1326   p2->Draw();p2->cd();
1327   p2->SetGridx(); p2->SetGridy();
1328   // 
1329   //Default cuts
1330   TCut cprim("cprim","MC.fVDist[3]<1");
1331   TCut cnprim("cnprim","MC.fVDist[3]>1");
1332   TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1333   TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1334   //
1335   c1->cd();  p1->cd();
1336   TH1F* hisp = ResPtvsPt("MC.fRowsWithDigits>100"+cteta1+cpos1,"1",0.15,2.,6);
1337   c1->cd(); 
1338   c1->Draw();
1339   p1->cd();
1340   p1->Draw();
1341   hisp->Draw(); 
1342   //
1343   c1->cd();  
1344   c1->Draw();
1345   p2->cd();
1346   p2->Draw();
1347   TH1F* his2 =  new TH1F("Ptresolution","Ptresolution",40,-5.,5.);
1348   fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt()>>Ptresolution","MC.fRowsWithDigits>100&&RC.fTPCTrack.fN>50&&RC.fMultiple==1"+cteta1+cpos1+cprim);
1349   AliLabelAxes(his2, "#Delta p_{t} / p_{t} [%]", "entries");
1350   his2->Fit("gaus");
1351   his2->Draw();
1352  
1353 }
1354
1355 void AliTPCComparisonDraw::Eff()
1356 {
1357   //
1358   //
1359   TCanvas *c1=new TCanvas("TPC efficiency","TPC efficiency",0,0,700,850);
1360   c1->Draw(); c1->cd();
1361   TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99); 
1362   p1->Draw();p1->cd(); 
1363   p1->SetGridx(); p1->SetGridy();
1364   //
1365   c1->cd();
1366   TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49); 
1367   p2->Draw();p2->cd();
1368   p2->SetGridx(); p2->SetGridy();
1369   // 
1370   //Default cuts
1371   TCut cprim("cprim","MC.fVDist[3]<1");
1372   TCut cnprim("cnprim","MC.fVDist[3]>1");
1373   TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1374   TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1375   //
1376   c1->cd();  p1->cd();
1377   TH1F* hisp = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cprim,"RC.fTPCTrack.fN>50");
1378   hisp->Draw();
1379   //hisp->DrawClone(); 
1380   TText * text = new TText(0.25,102.,"Primary particles");
1381   text->SetTextSize(0.05);
1382   text->Draw();
1383   //
1384   c1->cd();  p2->cd();
1385   TH1F* hiss = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cnprim,"RC.fTPCTrack.fN>50");
1386   hiss->Draw();
1387   //hiss->DrawClone();
1388   text = new TText(0.25,102.,"Secondary particles");
1389   text->SetTextSize(0.05);
1390   text->Draw();
1391
1392  
1393 }
1394
1395
1396 TH1F * AliTPCComparisonDraw::EffVsPt(const char* selection, const char * quality, Float_t min, Float_t max)
1397 {
1398   //
1399   //
1400   Int_t nBins = 10;
1401   Double_t* bins = CreateLogBins(nBins, min, max);
1402   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, bins);
1403   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, bins);
1404   
1405   fTree->Draw("MC.fParticle.Pt()>>hGen", selection, "groff");
1406   char selectionRec[256];
1407   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1408   fTree->Draw("MC.fParticle.Pt()>>hRec", selectionRec, "groff");
1409
1410   TH1F* hEff = CreateEffHisto(hGen, hRec);
1411   AliLabelAxes(hEff, "p_{t} [GeV/c]", "#epsilon [%]");
1412
1413   delete hRec;
1414   delete hGen;
1415   delete[] bins;
1416   return hEff;
1417 }
1418
1419
1420 TH1F * AliTPCComparisonDraw::EffVsRM(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1421 {
1422   //
1423   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1424   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1425   //  
1426   fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hGen", selection, "groff");
1427   char selectionRec[256];
1428   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1429   fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hRec", selectionRec, "groff");
1430   //
1431   TH1F* hEff = CreateEffHisto(hGen, hRec);
1432   AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1433   //
1434   delete hRec;
1435   delete hGen;
1436   return hEff;
1437 }
1438
1439 TH1F * AliTPCComparisonDraw::EffVsRS(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1440 {
1441   //
1442   TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1443   TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1444   //  
1445   fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hGen", selection, "groff");
1446   char selectionRec[256];
1447   sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1448   fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hRec", selectionRec, "groff");
1449   //
1450   TH1F* hEff = CreateEffHisto(hGen, hRec);
1451   AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1452   //
1453   delete hRec;
1454   delete hGen;
1455   return hEff;
1456
1457 }
1458
1459 TH1F * AliTPCComparisonDraw::ResPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1460 {
1461   //
1462   Double_t* bins = CreateLogBins(nBins, min, max);
1463   Int_t nBinsRes = 30;
1464   Double_t maxRes = 10.;
1465   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1466   
1467   fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1468
1469   TH1F* hMean=0;
1470   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1471   AliLabelAxes(hRes, "p_{t} [GeV/c]", "#Delta p_{t} / p_{t} [%]");
1472   //
1473   delete hRes2;
1474   delete[] bins;
1475   return hRes;
1476
1477 }
1478
1479 TH1F * AliTPCComparisonDraw::MeanPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1480 {
1481   //
1482   Double_t* bins = CreateLogBins(nBins, min, max);
1483   Int_t nBinsRes = 30;
1484   Double_t maxRes = 10.;
1485   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1486   
1487   fTree->Draw("100.*(1./fTPCTrack.Get1Pt()-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1488
1489   TH1F* hMean=0;
1490   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1491   AliLabelAxes(hRes, "p_{t} [GeV/c]", "mean p_{t} / p_{t} [%]");
1492   //
1493   delete hRes2;
1494   delete[] bins;
1495   if (!hMean) return 0;
1496   return hMean;
1497
1498 }
1499
1500
1501 TH1F * AliTPCComparisonDraw::ResdEdxvsN(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1502 {
1503   //
1504   Int_t nBinsRes = 15;
1505
1506   TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, min, max, nBinsRes, 34, 60);
1507   
1508   fTree->Draw("RC.fTPCTrack.fdEdx/MC.fPrim:RC.fTPCTrack.fN>>hRes2", selection, "groff");
1509
1510   TH1F* hMean=0;
1511   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1512   AliLabelAxes(hRes, "N points", "sigma dEdx/Nprim [%]");
1513   //
1514   delete hRes2;
1515   return hRes;
1516
1517 }
1518
1519
1520
1521 Double_t* AliTPCComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1522 {
1523   Double_t* bins = new Double_t[nBins+1];
1524   bins[0] = xMin;
1525   Double_t factor = pow(xMax/xMin, 1./nBins);
1526   for (Int_t i = 1; i <= nBins; i++)
1527     bins[i] = factor * bins[i-1];
1528   return bins;
1529 }
1530
1531
1532
1533
1534 void AliTPCComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1535 {
1536   //
1537   histo->GetXaxis()->SetTitle(xAxisTitle);
1538   histo->GetYaxis()->SetTitle(yAxisTitle);
1539 }
1540
1541
1542 TH1F* AliTPCComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1543 {
1544   //
1545   Int_t nBins = hGen->GetNbinsX();
1546   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1547   hEff->SetTitle("");
1548   hEff->SetStats(kFALSE);
1549   hEff->SetMinimum(0.);
1550   hEff->SetMaximum(110.);
1551   //
1552   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1553     Double_t nGen = hGen->GetBinContent(iBin);
1554     Double_t nRec = hRec->GetBinContent(iBin);
1555     if (nGen > 0) {
1556       Double_t eff = nRec/nGen;
1557       hEff->SetBinContent(iBin, 100. * eff);
1558       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
1559       //      if (error == 0) error = sqrt(0.1/nGen);
1560       //
1561       if (error == 0) error = 0.0001;
1562       hEff->SetBinError(iBin, 100. * error);
1563     } else {
1564       hEff->SetBinContent(iBin, 100. * 0.5);
1565       hEff->SetBinError(iBin, 100. * 0.5);
1566     }
1567   }
1568   return hEff;
1569 }
1570
1571
1572
1573 TH1F* AliTPCComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
1574                      Bool_t overflowBinFits)
1575 {
1576   TVirtualPad* currentPad = gPad;
1577   TAxis* axis = hRes2->GetXaxis();
1578   Int_t nBins = axis->GetNbins();
1579   TH1F* hRes, *hMean;
1580   if (axis->GetXbins()->GetSize()){
1581     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1582     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1583   }
1584   else{
1585     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1586     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1587
1588   }
1589   hRes->SetStats(false);
1590   hRes->SetOption("E");
1591   hRes->SetMinimum(0.);
1592   //
1593   hMean->SetStats(false);
1594   hMean->SetOption("E");
1595  
1596   // create the fit function
1597   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1598   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1599   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1600   
1601   fitFunc->SetLineWidth(2);
1602   fitFunc->SetFillStyle(0);
1603   // create canvas for fits
1604   TCanvas* canBinFits = NULL;
1605   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1606   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1607   Int_t ny = (nPads-1) / nx + 1;
1608   if (drawBinFits) {
1609     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1610     if (canBinFits) delete canBinFits;
1611     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1612     canBinFits->Divide(nx, ny);
1613   }
1614
1615   // loop over x bins and fit projection
1616   Int_t dBin = ((overflowBinFits) ? 1 : 0);
1617   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1618     if (drawBinFits) canBinFits->cd(bin + dBin);
1619     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1620     //    
1621     if (hBin->GetEntries() > 5) {
1622       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1623       hBin->Fit(fitFunc,"s");
1624       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1625
1626       if (sigma > 0.){
1627         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1628         hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
1629       }
1630       else{
1631         hRes->SetBinContent(bin, 0.);
1632         hMean->SetBinContent(bin,0);
1633       }
1634       hRes->SetBinError(bin, fitFunc->GetParError(2));
1635       hMean->SetBinError(bin, fitFunc->GetParError(1));
1636       
1637       //
1638       //
1639
1640     } else {
1641       hRes->SetBinContent(bin, 0.);
1642       hRes->SetBinError(bin, 0.);
1643       hMean->SetBinContent(bin, 0.);
1644       hMean->SetBinError(bin, 0.);
1645     }
1646     
1647
1648     if (drawBinFits) {
1649       char name[256];
1650       if (bin == 0) {
1651         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1652       } else if (bin == nBins+1) {
1653         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1654       } else {
1655         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1656                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1657       }
1658       canBinFits->cd(bin + dBin);
1659       hBin->SetTitle(name);
1660       hBin->SetStats(kTRUE);
1661       hBin->DrawCopy("E");
1662       canBinFits->Update();
1663       canBinFits->Modified();
1664       canBinFits->Update();
1665     }
1666     
1667     delete hBin;
1668   }
1669
1670   delete fitFunc;
1671   currentPad->cd();
1672   *phMean = hMean;
1673   return hRes;
1674 }
1675
1676
1677