ab21613989fd265dd3c9a3311c6fb4c4415020a6
[u/mrichter/AliRoot.git] / STEER / AliGenInfo.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
20 Origin: marian.ivanov@cern.ch
21
22 Macro to generate comples MC information - used for Comparison later on
23 How to use it?
24
25 .L $ALICE_ROOT/STEER/AliGenInfo.C+
26 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
27 t->Exec();
28
29 */
30
31 #if !defined(__CINT__) || defined(__MAKECINT__)
32 #include <stdio.h>
33 #include <string.h>
34 //ROOT includes
35 #include "Rtypes.h"
36 #include "TFile.h"
37 #include "TTree.h"
38 #include "TChain.h"
39 #include "TCut.h"
40 #include "TString.h"
41 #include "TBenchmark.h"
42 #include "TStopwatch.h"
43 #include "TParticle.h"
44 #include "TSystem.h"
45 #include "TTimer.h"
46 #include "TVector3.h"
47 #include "TH1F.h"
48 #include "TH2F.h"
49 #include "TCanvas.h"
50 #include "TPad.h"
51 #include "TF1.h"
52
53 //ALIROOT includes
54 #include "AliRun.h"
55 #include "AliStack.h"
56 #include "AliSimDigits.h"
57 #include "AliTPCParam.h"
58 #include "AliTPC.h"
59 #include "AliTPCLoader.h"
60 #include "AliDetector.h"
61 #include "AliTrackReference.h"
62 #include "AliTPCParamSR.h"
63 #include "AliTracker.h"
64 #include "AliMagF.h"
65 #endif
66 #include "AliGenInfo.h" 
67 //
68 //
69
70 AliTPCParam * GetTPCParam(){
71   AliTPCParamSR * par = new AliTPCParamSR;
72   par->Update();
73   return par;
74 }
75
76
77 //_____________________________________________________________________________
78 Float_t TPCBetheBloch(Float_t bg)
79 {
80   //
81   // Bethe-Bloch energy loss formula
82   //
83   const Double_t kp1=0.76176e-1;
84   const Double_t kp2=10.632;
85   const Double_t kp3=0.13279e-4;
86   const Double_t kp4=1.8631;
87   const Double_t kp5=1.9479;
88
89   Double_t dbg = (Double_t) bg;
90
91   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
92
93   Double_t aa = TMath::Power(beta,kp4);
94   Double_t bb = TMath::Power(1./dbg,kp5);
95
96   bb=TMath::Log(kp3+bb);
97   
98   return ((Float_t)((kp2-aa-bb)*kp1/aa));
99 }
100
101
102
103
104
105 ////////////////////////////////////////////////////////////////////////
106 AliMCInfo::AliMCInfo()
107 {
108   fTPCReferences  = new TClonesArray("AliTrackReference",10);
109   fITSReferences  = new TClonesArray("AliTrackReference",10);
110   fTRDReferences  = new TClonesArray("AliTrackReference",10);
111   fTRdecay.SetTrack(-1);
112 }
113
114 AliMCInfo::~AliMCInfo()
115 {
116   if (fTPCReferences) {
117     delete fTPCReferences;
118   }
119   if (fITSReferences){
120     delete fITSReferences;
121   }
122   if (fTRDReferences){    
123     delete fTRDReferences;  
124   }
125 }
126
127
128
129 void AliMCInfo::Update()
130 {
131   //
132   //
133   fMCtracks =1;
134   if (!fTPCReferences) {
135     fNTPCRef =0;
136     return;
137   }
138   Float_t direction=1;
139   //Float_t rlast=0;
140   fNTPCRef = fTPCReferences->GetEntriesFast();
141   fNITSRef = fITSReferences->GetEntriesFast();
142
143   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
144     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
145     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
146     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
147     if (iref==0) direction = newdirection;
148     if ( newdirection*direction<0){
149       //changed direction
150       direction = newdirection;
151       fMCtracks+=1;
152     }
153     //rlast=r;                      
154   }
155   //
156   // decay info
157   fTPCdecay=kFALSE;
158   if (fTRdecay.GetTrack()>0){
159     fDecayCoord[0] = fTRdecay.X();
160     fDecayCoord[1] = fTRdecay.Y();
161     fDecayCoord[2] = fTRdecay.Z();
162     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
163       fTPCdecay=kTRUE;     
164     }
165     else{
166       fDecayCoord[0] = 0;
167       fDecayCoord[1] = 0;
168       fDecayCoord[2] = 0;
169     }
170   }
171   //
172   //
173   //digits information update
174   fRowsWithDigits    = fTPCRow.RowsOn();    
175   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
176   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
177   //
178   //
179   // calculate primary ionization per cm  
180   if (fParticle.GetPDG()){
181     fMass = fParticle.GetMass();  
182     fCharge = fParticle.GetPDG()->Charge();
183     if (fTPCReferences->GetEntriesFast()>0){
184       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
185     }
186     if (fMass>0){
187       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
188                               fTrackRef.Py()*fTrackRef.Py()+
189                               fTrackRef.Pz()*fTrackRef.Pz());    
190       if (p>0.001){
191         Float_t betagama = p /fMass;
192         fPrim = TPCBetheBloch(betagama);
193       }else fPrim=0;
194     }
195   }else{
196     fMass =0;
197     fPrim =0;
198   }  
199 }
200
201 /////////////////////////////////////////////////////////////////////////////////
202 void AliGenV0Info::Update()
203 {
204   fMCPd[0] = fMCd.fParticle.Px();
205   fMCPd[1] = fMCd.fParticle.Py();
206   fMCPd[2] = fMCd.fParticle.Pz();
207   fMCPd[3] = fMCd.fParticle.P();
208   fMCX[0]  = fMCd.fParticle.Vx();
209   fMCX[1]  = fMCd.fParticle.Vy();
210   fMCX[2]  = fMCd.fParticle.Vz();
211   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
212   fPdg[0]    = fMCd.fParticle.GetPdgCode();
213   fPdg[1]    = fMCm.fParticle.GetPdgCode();
214   //
215   fLab[0]    = fMCd.fParticle.GetUniqueID();
216   fLab[1]    = fMCm.fParticle.GetUniqueID();
217
218 }
219
220
221
222
223   
224 ////////////////////////////////////////////////////////////////////////
225
226 ////////////////////////////////////////////////////////////////////////
227 //
228 // End of implementation of the class AliMCInfo
229 //
230 ////////////////////////////////////////////////////////////////////////
231
232
233
234 ////////////////////////////////////////////////////////////////////////
235 digitRow::digitRow()
236 {
237   Reset();
238 }
239 ////////////////////////////////////////////////////////////////////////
240 digitRow & digitRow::operator=(const digitRow &digOld)
241 {
242   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
243   return (*this);
244 }
245 ////////////////////////////////////////////////////////////////////////
246 void digitRow::SetRow(Int_t row) 
247 {
248   if (row >= 8*kgRowBytes) {
249     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
250     return;
251   }
252   Int_t iC = row/8;
253   Int_t iB = row%8;
254   SETBIT(fDig[iC],iB);
255 }
256
257 ////////////////////////////////////////////////////////////////////////
258 Bool_t digitRow::TestRow(Int_t row)
259 {
260 //
261 // return kTRUE if row is on
262 //
263   Int_t iC = row/8;
264   Int_t iB = row%8;
265   return TESTBIT(fDig[iC],iB);
266 }
267 ////////////////////////////////////////////////////////////////////////
268 Int_t digitRow::RowsOn(Int_t upto)
269 {
270 //
271 // returns number of rows with a digit  
272 // count only rows less equal row number upto
273 //
274   Int_t total = 0;
275   for (Int_t i = 0; i<kgRowBytes; i++) {
276     for (Int_t j = 0; j < 8; j++) {
277       if (i*8+j > upto) return total;
278       if (TESTBIT(fDig[i],j))  total++;
279     }
280   }
281   return total;
282 }
283 ////////////////////////////////////////////////////////////////////////
284 void digitRow::Reset()
285 {
286 //
287 // resets all rows to zero
288 //
289   for (Int_t i = 0; i<kgRowBytes; i++) {
290     fDig[i] <<= 8;
291   }
292 }
293 ////////////////////////////////////////////////////////////////////////
294 Int_t digitRow::Last()
295 {
296 //
297 // returns the last row number with a digit
298 // returns -1 if now digits 
299 //
300   for (Int_t i = kgRowBytes-1; i>=0; i--) {
301     for (Int_t j = 7; j >= 0; j--) {
302       if TESTBIT(fDig[i],j) return i*8+j;
303     }
304   }
305   return -1;
306 }
307 ////////////////////////////////////////////////////////////////////////
308 Int_t digitRow::First()
309 {
310 //
311 // returns the first row number with a digit
312 // returns -1 if now digits 
313 //
314   for (Int_t i = 0; i<kgRowBytes; i++) {
315     for (Int_t j = 0; j < 8; j++) {
316       if (TESTBIT(fDig[i],j)) return i*8+j;
317     }
318   }
319   return -1;
320 }
321
322 ////////////////////////////////////////////////////////////////////////
323 //
324 // end of implementation of a class digitRow
325 //
326 ////////////////////////////////////////////////////////////////////////
327   
328 ////////////////////////////////////////////////////////////////////////
329 AliGenInfoMaker::AliGenInfoMaker()
330 {
331   Reset();
332 }
333
334 ////////////////////////////////////////////////////////////////////////
335 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
336                                    Int_t nEvents, Int_t firstEvent)
337 {
338   Reset();
339   fFirstEventNr = firstEvent;
340   fEventNr = firstEvent;
341   fNEvents = nEvents;
342   //   fFnRes = fnRes;
343   sprintf(fFnRes,"%s",fnRes);
344   //
345   fLoader = AliRunLoader::Open(fnGalice);
346   if (gAlice){
347     delete gAlice->GetRunLoader();
348     delete gAlice;
349     gAlice = 0x0;
350   }
351   if (fLoader->LoadgAlice()){
352     cerr<<"Error occured while l"<<endl;
353   }
354   Int_t nall = fLoader->GetNumberOfEvents();
355   if (nEvents==0) {
356     nEvents =nall;
357     fNEvents=nall;
358     fFirstEventNr=0;
359   }    
360
361   if (nall<=0){
362     cerr<<"no events available"<<endl;
363     fEventNr = 0;
364     return;
365   }
366   if (firstEvent+nEvents>nall) {
367     fEventNr = nall-firstEvent;
368     cerr<<"restricted number of events availaible"<<endl;
369   }
370   AliMagF * magf = gAlice->Field();
371   AliTracker::SetFieldMap(magf);
372 }
373
374
375 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
376 {
377   // 
378   if (i<fNParticles) {
379     if (fGenInfo[i]) return  fGenInfo[i];
380     fGenInfo[i] = new AliMCInfo;  
381     fNInfos++;
382     return fGenInfo[i];
383   }
384   else 
385     return 0;  
386 }
387
388 ////////////////////////////////////////////////////////////////////////
389 void AliGenInfoMaker::Reset()
390 {
391   fEventNr = 0;
392   fNEvents = 0;
393   fTreeGenTracks = 0;
394   fFileGenTracks = 0;
395   fGenInfo = 0;
396   fNInfos  = 0;
397   //
398   //
399   fDebug = 0;
400   fVPrim[0] = -1000.;
401   fVPrim[1] = -1000.;
402   fVPrim[2] = -1000.;
403   fParamTPC = 0;
404 }
405 ////////////////////////////////////////////////////////////////////////
406 AliGenInfoMaker::~AliGenInfoMaker()
407 {
408   
409   if (fLoader){
410     fLoader->UnloadgAlice();
411     gAlice = 0;
412     delete fLoader;
413   }
414 }
415
416 Int_t  AliGenInfoMaker::SetIO()
417 {
418   //
419   // 
420   CreateTreeGenTracks();
421   if (!fTreeGenTracks) return 1;
422   //  AliTracker::SetFieldFactor(); 
423  
424   fParamTPC = GetTPCParam();
425   //
426   return 0;
427 }
428
429 ////////////////////////////////////////////////////////////////////////
430 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
431 {
432   //
433   // 
434   // SET INPUT
435   fLoader->SetEventNumber(eventNr);
436   //
437   fLoader->LoadHeader();
438   fLoader->LoadKinematics();  
439   fStack = fLoader->Stack();
440   //
441   fLoader->LoadTrackRefs();
442   fTreeTR = fLoader->TreeTR();
443   //
444   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
445   tpcl->LoadDigits();
446   fTreeD = tpcl->TreeD();  
447   return 0;
448 }
449
450 Int_t AliGenInfoMaker::CloseIOEvent()
451 {
452   fLoader->UnloadHeader();
453   fLoader->UnloadKinematics();
454   fLoader->UnloadTrackRefs();
455   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
456   tpcl->UnloadDigits();
457   return 0;
458 }
459
460 Int_t AliGenInfoMaker::CloseIO()
461 {
462   fLoader->UnloadgAlice();
463   return 0;
464 }
465
466
467
468 ////////////////////////////////////////////////////////////////////////
469 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
470 {
471   fNEvents = nEvents;
472   fFirstEventNr = firstEventNr;
473   return Exec();
474 }
475
476 ////////////////////////////////////////////////////////////////////////
477 Int_t AliGenInfoMaker::Exec()  
478 {
479   TStopwatch timer;
480   timer.Start();
481   Int_t status =SetIO();
482   if (status>0) return status;
483   //
484
485   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
486        fEventNr++) {
487     SetIO(fEventNr);
488     fNParticles = fStack->GetNtrack();
489     //
490     fGenInfo = new AliMCInfo*[fNParticles];
491     for (UInt_t i = 0; i<fNParticles; i++) {
492       fGenInfo[i]=0; 
493     }
494     //
495     cout<<"Start to process event "<<fEventNr<<endl;
496     cout<<"\tfNParticles = "<<fNParticles<<endl;
497     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
498     if (TreeTRLoop()>0) return 1;
499     //
500     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
501     if (TreeDLoop()>0) return 1;
502     //
503     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
504     if (TreeKLoop()>0) return 1;
505     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
506     for (UInt_t i = 0; i<fNParticles; i++) {
507       if (fGenInfo[i]) delete fGenInfo[i]; 
508     }
509     delete []fGenInfo;
510     CloseIOEvent();
511   }
512   //
513   CloseIO();
514   CloseOutputFile();
515
516   cerr<<"Exec finished"<<endl;
517
518   timer.Stop();
519   timer.Print();
520   return 0;
521 }
522 ////////////////////////////////////////////////////////////////////////
523 void AliGenInfoMaker::CreateTreeGenTracks() 
524 {
525   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
526   if (!fFileGenTracks) {
527     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
528     return;
529   }
530   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
531   AliMCInfo * info = new AliMCInfo;
532   //
533   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
534   delete info; 
535   fTreeGenTracks->AutoSave();
536 }
537 ////////////////////////////////////////////////////////////////////////
538 void AliGenInfoMaker::CloseOutputFile() 
539 {
540   if (!fFileGenTracks) {
541     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
542     return;
543   }
544   fFileGenTracks->cd();
545   fTreeGenTracks->Write();  
546   delete fTreeGenTracks;
547   fFileGenTracks->Close();
548   delete fFileGenTracks;
549   return;
550 }
551
552 ////////////////////////////////////////////////////////////////////////
553 Int_t AliGenInfoMaker::TreeKLoop()
554 {
555 //
556 // open the file with treeK
557 // loop over all entries there and save information about some tracks
558 //
559
560   AliStack * stack = fStack;
561   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
562   
563   if (fDebug > 0) {
564     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
565         <<fEventNr<<endl;
566   }  
567   Int_t  ipdg = 0;
568   TParticlePDG *ppdg = 0;
569   // not all generators give primary vertex position. Take the vertex
570   // of the particle 0 as primary vertex.
571   TDatabasePDG  pdg; //get pdg table  
572   //thank you very much root for this
573   TBranch * br = fTreeGenTracks->GetBranch("MC");
574   TParticle *particle = stack->ParticleFromTreeK(0);
575   fVPrim[0] = particle->Vx();
576   fVPrim[1] = particle->Vy();
577   fVPrim[2] = particle->Vz();
578   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
579     // load only particles with TR
580     AliMCInfo * info = GetInfo(iParticle);
581     if (!info) continue;
582     //////////////////////////////////////////////////////////////////////
583     info->fLabel = iParticle;
584     //
585     info->fParticle = *(stack->Particle(iParticle));
586     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
587     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
588     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; 
589     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
590                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
591     //
592     //
593     ipdg = info->fParticle.GetPdgCode();
594     info->fPdg = ipdg;
595     ppdg = pdg.GetParticle(ipdg);          
596     info->fEventNr = fEventNr;
597     info->Update();
598     //////////////////////////////////////////////////////////////////////    
599     br->SetAddress(&info);    
600     fTreeGenTracks->Fill();
601   }
602   fTreeGenTracks->AutoSave();
603   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
604   return 0;
605 }
606
607
608
609
610 ////////////////////////////////////////////////////////////////////////
611 Int_t AliGenInfoMaker::TreeDLoop()
612 {
613   //
614   // open the file with treeD
615   // loop over all entries there and save information about some tracks
616   //
617   
618   Int_t nInnerSector = fParamTPC->GetNInnerSector();
619   Int_t rowShift = 0;
620   Int_t zero=fParamTPC->GetZeroSup()+6;  
621   //
622   //
623   AliSimDigits digitsAddress, *digits=&digitsAddress;
624   fTreeD->GetBranch("Segment")->SetAddress(&digits);
625   
626   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
627   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
628   for (Int_t i=0; i<sectorsByRows; i++) {
629     if (!fTreeD->GetEvent(i)) continue;
630     Int_t sec,row;
631     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
632     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
633     // here I expect that upper sectors follow lower sectors
634     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
635     //
636     digits->ExpandTrackBuffer();
637     digits->First();        
638     do {
639       Int_t iRow=digits->CurrentRow();
640       Int_t iColumn=digits->CurrentColumn();
641       Short_t digitValue = digits->CurrentDigit();
642       if (digitValue >= zero) {
643         Int_t label;
644         for (Int_t j = 0; j<3; j++) {
645           //      label = digits->GetTrackID(iRow,iColumn,j);
646           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;       
647           if (label >= (Int_t)fNParticles) { //don't label from bakground event
648             continue;
649           }
650           if (label >= 0 && label <= (Int_t)fNParticles) {
651             if (fDebug > 6 ) {
652               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
653                   <<sec<<" "
654                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
655                   <<" "<<row<<endl;
656             }   
657             AliMCInfo * info = GetInfo(label);
658             if (info){
659               info->fTPCRow.SetRow(row+rowShift);
660             }
661           }
662         }
663       }
664     } while (digits->Next());
665   }
666   
667   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;  
668   return 0;
669 }
670
671
672 ////////////////////////////////////////////////////////////////////////
673 Int_t AliGenInfoMaker::TreeTRLoop()
674 {
675   //
676   // loop over TrackReferences and store the first one for each track
677   //  
678   TTree * treeTR = fTreeTR;
679   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
680   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
681   //
682   //
683   //track references for TPC
684   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
685   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
686   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
687   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
688   //
689   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
690   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
691   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
692   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
693   //
694   //
695   //
696   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
697     treeTR->GetEntry(iPrimPart);
698     //
699     // Loop over TPC references
700     //
701     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
702       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
703       //
704       if (trackRef->Pt()<fgTPCPtCut) continue;
705       Int_t label = trackRef->GetTrack();      
706       AliMCInfo * info = GetInfo(label);
707       if (!info) info = MakeInfo(label);
708       if (!info) continue;
709       TClonesArray & arr = *(info->fTPCReferences);
710       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
711     }
712     //
713     // Loop over ITS references
714     //
715     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
716       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);            
717       //
718       if (trackRef->Pt()<fgTPCPtCut) continue;
719       Int_t label = trackRef->GetTrack();      
720       AliMCInfo * info = GetInfo(label);
721       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
722       if (!info) info = MakeInfo(label);
723       if (!info) continue;
724       TClonesArray & arr = *(info->fITSReferences);
725       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
726     }
727     //
728     // Loop over TRD references
729     //
730     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
731       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);            
732       //
733       if (trackRef->Pt()<fgTPCPtCut) continue;
734       Int_t label = trackRef->GetTrack();      
735       AliMCInfo * info = GetInfo(label);
736       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
737       if (!info) info = MakeInfo(label);
738       if (!info) continue;
739       TClonesArray & arr = *(info->fTRDReferences);
740       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
741     }
742     //
743     // get dacay position
744     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
745       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
746       //
747       Int_t label = trackRef->GetTrack();
748       AliMCInfo * info = GetInfo(label);
749       if (!info) continue;
750       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
751       //      if (TMath::Abs(trackRef.X());
752       info->fTRdecay = *trackRef;      
753     }
754   }
755   //
756   TPCArrayTR->Delete();
757   delete  TPCArrayTR;
758   TRDArrayTR->Delete();
759   delete  TRDArrayTR;
760   ITSArrayTR->Delete();
761   delete  ITSArrayTR;
762   RunArrayTR->Delete();
763   delete  RunArrayTR;
764   //
765   return 0;
766 }
767
768 ////////////////////////////////////////////////////////////////////////
769 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
770                                     AliTPCParam *paramTPC) {
771
772   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
773   Int_t index[4];
774   paramTPC->Transform0to1(x,index);
775   paramTPC->Transform1to2(x,index);
776   return x[0];
777 }
778 ////////////////////////////////////////////////////////////////////////
779
780
781
782 TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection, 
783                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
784 {
785   //
786   Double_t* bins = CreateLogBins(nbins, minx, maxx);
787   Int_t nBinsRes = 30;
788   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
789   char cut[1000];
790   sprintf(cut,"%s&&%s",selection,quality);
791   char expression[1000];
792   sprintf(expression,"%s:%s>>hRes2",chy,chx);
793   fTree->Draw(expression, cut, "groff");
794   TH1F* hMean=0;
795   TH1F* hRes = CreateResHisto(hRes2, &hMean);
796   AliLabelAxes(hRes, chx, chy);
797   //
798   delete hRes2;
799   delete[] bins;
800   fRes  = hRes;
801   fMean = hMean;
802   return hRes;
803 }
804
805 TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, 
806                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
807 {
808   //
809   Double_t* bins = CreateLogBins(nbins, minx, maxx);
810   Int_t nBinsRes = 30;
811   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
812   char cut[1000];
813   sprintf(cut,"%s&&%s",selection,quality);
814   char expression[1000];
815   sprintf(expression,"%s:%s>>hRes2",chy,chx);
816   fTree->Draw(expression, cut, "groff");
817   TH1F* hMean=0;  
818   TH1F* hRes = CreateResHisto(hRes2, &hMean);
819   AliLabelAxes(hRes, chx, chy);
820   //
821   delete hRes2;
822   delete[] bins;
823   fRes  = hRes;
824   fMean = hMean;
825   return hRes;
826 }
827
828 ///////////////////////////////////////////////////////////////////////////////////
829 ///////////////////////////////////////////////////////////////////////////////////
830 TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality, 
831                               Int_t nbins, Float_t min, Float_t max)
832 {
833   //
834   //
835   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
836   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
837   char inputGen[1000];  
838   sprintf(inputGen,"%s>>hGen", variable);
839   fTree->Draw(inputGen, selection, "groff");
840   char selectionRec[256];
841   sprintf(selectionRec, "%s && %s", selection, quality);
842   char inputRec[1000];  
843   sprintf(inputRec,"%s>>hRec", variable);
844   fTree->Draw(inputRec, selectionRec, "groff");
845   //
846   TH1F* hEff = CreateEffHisto(hGen, hRec);
847   AliLabelAxes(hEff, variable, "#epsilon [%]");
848   fRes = hEff;
849   delete hRec;
850   delete hGen;
851   return hEff;
852 }
853
854
855
856 ///////////////////////////////////////////////////////////////////////////////////
857 ///////////////////////////////////////////////////////////////////////////////////
858 TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality, 
859                               Int_t nbins, Float_t min, Float_t max)
860 {
861   //
862   //
863   Double_t* bins = CreateLogBins(nbins, min, max);
864   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
865   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
866   char inputGen[1000];  
867   sprintf(inputGen,"%s>>hGen", variable);
868   fTree->Draw(inputGen, selection, "groff");
869   char selectionRec[256];
870   sprintf(selectionRec, "%s && %s", selection, quality);
871   char inputRec[1000];  
872   sprintf(inputRec,"%s>>hRec", variable);
873   fTree->Draw(inputRec, selectionRec, "groff");
874   //
875   TH1F* hEff = CreateEffHisto(hGen, hRec);
876   AliLabelAxes(hEff, variable, "#epsilon [%]");
877   fRes = hEff;
878   delete hRec;
879   delete hGen;
880   delete[] bins;
881   return hEff;
882 }
883
884
885 ///////////////////////////////////////////////////////////////////////////////////
886 ///////////////////////////////////////////////////////////////////////////////////
887
888 Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
889 {
890   Double_t* bins = new Double_t[nBins+1];
891   bins[0] = xMin;
892   Double_t factor = pow(xMax/xMin, 1./nBins);
893   for (Int_t i = 1; i <= nBins; i++)
894     bins[i] = factor * bins[i-1];
895   return bins;
896 }
897
898
899
900
901 void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
902 {
903   //
904   histo->GetXaxis()->SetTitle(xAxisTitle);
905   histo->GetYaxis()->SetTitle(yAxisTitle);
906 }
907
908
909 TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
910 {
911   //
912   Int_t nBins = hGen->GetNbinsX();
913   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
914   hEff->SetTitle("");
915   hEff->SetStats(kFALSE);
916   hEff->SetMinimum(0.);
917   hEff->SetMaximum(110.);
918   //
919   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
920     Double_t nGen = hGen->GetBinContent(iBin);
921     Double_t nRec = hRec->GetBinContent(iBin);
922     if (nGen > 0) {
923       Double_t eff = nRec/nGen;
924       hEff->SetBinContent(iBin, 100. * eff);
925       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
926       //      if (error == 0) error = sqrt(0.1/nGen);
927       //
928       if (error == 0) error = 0.0001;
929       hEff->SetBinError(iBin, 100. * error);
930     } else {
931       hEff->SetBinContent(iBin, 100. * 0.5);
932       hEff->SetBinError(iBin, 100. * 0.5);
933     }
934   }
935   return hEff;
936 }
937
938
939
940 TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
941                      Bool_t overflowBinFits)
942 {
943   TVirtualPad* currentPad = gPad;
944   TAxis* axis = hRes2->GetXaxis();
945   Int_t nBins = axis->GetNbins();
946   TH1F* hRes, *hMean;
947   if (axis->GetXbins()->GetSize()){
948     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
949     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
950   }
951   else{
952     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
953     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
954
955   }
956   hRes->SetStats(false);
957   hRes->SetOption("E");
958   hRes->SetMinimum(0.);
959   //
960   hMean->SetStats(false);
961   hMean->SetOption("E");
962  
963   // create the fit function
964   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
965   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
966   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
967   
968   fitFunc->SetLineWidth(2);
969   fitFunc->SetFillStyle(0);
970   // create canvas for fits
971   TCanvas* canBinFits = NULL;
972   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
973   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
974   Int_t ny = (nPads-1) / nx + 1;
975   if (drawBinFits) {
976     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
977     if (canBinFits) delete canBinFits;
978     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
979     canBinFits->Divide(nx, ny);
980   }
981
982   // loop over x bins and fit projection
983   Int_t dBin = ((overflowBinFits) ? 1 : 0);
984   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
985     if (drawBinFits) canBinFits->cd(bin + dBin);
986     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
987     //    
988     if (hBin->GetEntries() > 5) {
989       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
990       hBin->Fit(fitFunc,"s");
991       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
992
993       if (sigma > 0.){
994         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
995         hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
996       }
997       else{
998         hRes->SetBinContent(bin, 0.);
999         hMean->SetBinContent(bin,0);
1000       }
1001       hRes->SetBinError(bin, fitFunc->GetParError(2));
1002       hMean->SetBinError(bin, fitFunc->GetParError(1));
1003       
1004       //
1005       //
1006
1007     } else {
1008       hRes->SetBinContent(bin, 0.);
1009       hRes->SetBinError(bin, 0.);
1010       hMean->SetBinContent(bin, 0.);
1011       hMean->SetBinError(bin, 0.);
1012     }
1013     
1014
1015     if (drawBinFits) {
1016       char name[256];
1017       if (bin == 0) {
1018         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1019       } else if (bin == nBins+1) {
1020         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1021       } else {
1022         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1023                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1024       }
1025       canBinFits->cd(bin + dBin);
1026       hBin->SetTitle(name);
1027       hBin->SetStats(kTRUE);
1028       hBin->DrawCopy("E");
1029       canBinFits->Update();
1030       canBinFits->Modified();
1031       canBinFits->Update();
1032     }
1033     
1034     delete hBin;
1035   }
1036
1037   delete fitFunc;
1038   currentPad->cd();
1039   *phMean = hMean;
1040   return hRes;
1041 }
1042
1043
1044