Adding new code for kink and V0 reconstruction (AliESDkink and AliESDV0MI classes...
[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 #include "TView.h"
53 #include "TGeometry.h"
54 #include "TPolyLine3D.h"
55
56 //ALIROOT includes
57 #include "AliRun.h"
58 #include "AliStack.h"
59 #include "AliSimDigits.h"
60 #include "AliTPCParam.h"
61 #include "AliTPC.h"
62 #include "AliTPCLoader.h"
63 #include "AliDetector.h"
64 #include "AliTrackReference.h"
65 #include "AliTPCParamSR.h"
66 #include "AliTracker.h"
67 #include "AliMagF.h"
68 #include "AliHelix.h"
69 #include "AliPoints.h"
70
71 #endif
72 #include "AliGenInfo.h" 
73 //
74 // 
75
76 AliTPCParam * GetTPCParam(){
77   AliTPCParamSR * par = new AliTPCParamSR;
78   par->Update();
79   return par;
80 }
81
82
83 //_____________________________________________________________________________
84 Float_t TPCBetheBloch(Float_t bg)
85 {
86   //
87   // Bethe-Bloch energy loss formula
88   //
89   const Double_t kp1=0.76176e-1;
90   const Double_t kp2=10.632;
91   const Double_t kp3=0.13279e-4;
92   const Double_t kp4=1.8631;
93   const Double_t kp5=1.9479;
94
95   Double_t dbg = (Double_t) bg;
96
97   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
98
99   Double_t aa = TMath::Power(beta,kp4);
100   Double_t bb = TMath::Power(1./dbg,kp5);
101
102   bb=TMath::Log(kp3+bb);
103   
104   return ((Float_t)((kp2-aa-bb)*kp1/aa));
105 }
106
107 AliPointsMI::AliPointsMI(){
108   fN=0;
109   fX=0;
110   fY=0;
111   fZ=0;  
112   fCapacity = 0;
113   fLabel0=0;
114   fLabel1=0;
115 }
116
117
118 AliPointsMI::AliPointsMI(Int_t n, Float_t *x,Float_t *y, Float_t *z){
119   //
120   //
121   fN=n;
122   fCapacity = 1000+2*n;
123   fX= new Float_t[n];
124   fY= new Float_t[n];
125   fZ= new Float_t[n];
126   memcpy(fX,x,n*sizeof(Float_t));
127   memcpy(fY,y,n*sizeof(Float_t));
128   memcpy(fZ,z,n*sizeof(Float_t));
129   fLabel0=0;
130   fLabel1=0;
131 }
132
133 void AliPointsMI::Reset()
134 {
135   fN=0;
136 }
137
138 void AliPointsMI::Reset(AliDetector * det, Int_t particle){
139   //
140   // write points from detector points
141   //
142   Reset();
143   if (!det) return;
144   TObjArray *array = det->Points();
145   if (!array) return;
146   for (Int_t i=0;i<array->GetEntriesFast();i++){
147     AliPoints * points = (AliPoints*) array->At(i);
148     if (!points) continue;
149     if (points->GetIndex()!=particle) continue;
150     Int_t npoints = points->GetN();
151     if (npoints<2) continue;
152     Int_t delta = npoints/100;
153     if (delta<1) delta=1;
154     if (delta>10) delta=10;
155     Int_t mypoints = npoints/delta;
156     //
157     fN = mypoints;
158     if (fN>fCapacity){
159       fCapacity = 1000+2*fN;
160       delete []fX;
161       delete []fY;
162       delete []fZ;
163       fX = new Float_t[fCapacity];
164       fY = new Float_t[fCapacity];
165       fZ = new Float_t[fCapacity]; 
166     }   
167     Float_t *xyz = points->GetP();
168     for (Int_t ipoint=0;ipoint<mypoints;ipoint++){
169       Int_t index = 3*ipoint*delta;
170       fX[ipoint]=0;
171       fY[ipoint]=0;
172       fZ[ipoint]=0;
173       if (index+2<npoints*3){
174         fX[ipoint] = xyz[index];
175         fY[ipoint] = xyz[index+1];
176         fZ[ipoint] = xyz[index+2];
177       }
178     }
179   }
180   fLabel0 = particle;
181 }
182
183
184 AliPointsMI::~AliPointsMI(){
185   fN=0;
186   fCapacity =0;
187   delete[] fX;
188   delete[]fY;
189   delete []fZ;
190   fX=0;fY=0;fZ=0;  
191 }
192
193
194
195
196
197 ////////////////////////////////////////////////////////////////////////
198 AliMCInfo::AliMCInfo()
199 {
200   fTPCReferences  = new TClonesArray("AliTrackReference",10);
201   fITSReferences  = new TClonesArray("AliTrackReference",10);
202   fTRDReferences  = new TClonesArray("AliTrackReference",10);
203   fTOFReferences  = new TClonesArray("AliTrackReference",10);
204   fTRdecay.SetTrack(-1);
205 }
206
207 AliMCInfo::~AliMCInfo()
208 {
209   if (fTPCReferences) {
210     delete fTPCReferences;
211   }
212   if (fITSReferences){
213     delete fITSReferences;
214   }
215   if (fTRDReferences){    
216     delete fTRDReferences;  
217   }
218   if (fTOFReferences){    
219     delete fTOFReferences;  
220   }
221
222 }
223
224
225
226 void AliMCInfo::Update()
227 {
228   //
229   //
230   fMCtracks =1;
231   if (!fTPCReferences) {
232     fNTPCRef =0;
233     return;
234   }
235   Float_t direction=1;
236   //Float_t rlast=0;
237   fNTPCRef = fTPCReferences->GetEntriesFast();
238   fNITSRef = fITSReferences->GetEntriesFast();
239   fNTRDRef = fTRDReferences->GetEntriesFast();
240   fNTOFRef = fTOFReferences->GetEntriesFast();
241   
242   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
243     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
244     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
245     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
246     if (iref==0) direction = newdirection;
247     if ( newdirection*direction<0){
248       //changed direction
249       direction = newdirection;
250       fMCtracks+=1;
251     }
252     //rlast=r;                      
253   }
254   //
255   // decay info
256   fTPCdecay=kFALSE;
257   if (fTRdecay.GetTrack()>0){
258     fDecayCoord[0] = fTRdecay.X();
259     fDecayCoord[1] = fTRdecay.Y();
260     fDecayCoord[2] = fTRdecay.Z();
261     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
262       fTPCdecay=kTRUE;     
263     }
264     else{
265       fDecayCoord[0] = 0;
266       fDecayCoord[1] = 0;
267       fDecayCoord[2] = 0;
268     }
269   }
270   //
271   //
272   //digits information update
273   fRowsWithDigits    = fTPCRow.RowsOn();    
274   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
275   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
276   //
277   //
278   // calculate primary ionization per cm  
279   if (fParticle.GetPDG()){
280     fMass = fParticle.GetMass();  
281     fCharge = fParticle.GetPDG()->Charge();
282     if (fTPCReferences->GetEntriesFast()>0){
283       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
284     }
285     if (fMass>0){
286       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
287                               fTrackRef.Py()*fTrackRef.Py()+
288                               fTrackRef.Pz()*fTrackRef.Pz());    
289       if (p>0.001){
290         Float_t betagama = p /fMass;
291         fPrim = TPCBetheBloch(betagama);
292       }else fPrim=0;
293     }
294   }else{
295     fMass =0;
296     fPrim =0;
297   }  
298 }
299
300 /////////////////////////////////////////////////////////////////////////////////
301 /*
302 void AliGenV0Info::Update()
303 {
304   fMCPd[0] = fMCd.fParticle.Px();
305   fMCPd[1] = fMCd.fParticle.Py();
306   fMCPd[2] = fMCd.fParticle.Pz();
307   fMCPd[3] = fMCd.fParticle.P();
308   //
309   fMCX[0]  = fMCd.fParticle.Vx();
310   fMCX[1]  = fMCd.fParticle.Vy();
311   fMCX[2]  = fMCd.fParticle.Vz();
312   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
313   //
314   fPdg[0]    = fMCd.fParticle.GetPdgCode();
315   fPdg[1]    = fMCm.fParticle.GetPdgCode();
316   //
317   fLab[0]    = fMCd.fParticle.GetUniqueID();
318   fLab[1]    = fMCm.fParticle.GetUniqueID();
319   //  
320 }
321 */
322
323 void AliGenV0Info::Update(Float_t vertex[3])
324 {
325   fMCPd[0] = fMCd.fParticle.Px();
326   fMCPd[1] = fMCd.fParticle.Py();
327   fMCPd[2] = fMCd.fParticle.Pz();
328   fMCPd[3] = fMCd.fParticle.P();
329   //
330   fMCX[0]  = fMCd.fParticle.Vx();
331   fMCX[1]  = fMCd.fParticle.Vy();
332   fMCX[2]  = fMCd.fParticle.Vz();
333   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
334   //
335   fPdg[0]    = fMCd.fParticle.GetPdgCode();
336   fPdg[1]    = fMCm.fParticle.GetPdgCode();
337   //
338   fLab[0]    = fMCd.fParticle.GetUniqueID();
339   fLab[1]    = fMCm.fParticle.GetUniqueID();
340   //
341   //
342   //
343   Double_t x1[3],p1[3];
344   TParticle & pdaughter = fMCd.fParticle;
345   x1[0] = pdaughter.Vx();      
346   x1[1] = pdaughter.Vy();
347   x1[2] = pdaughter.Vz();
348   p1[0] = pdaughter.Px();
349   p1[1] = pdaughter.Py();
350   p1[2] = pdaughter.Pz();
351   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
352   AliHelix dhelix1(x1,p1,sign);
353   //
354   //
355   Double_t x2[3],p2[3];            
356   //
357   TParticle & pmother = fMCm.fParticle;
358   x2[0] = pmother.Vx();      
359   x2[1] = pmother.Vy();      
360   x2[2] = pmother.Vz();      
361   p2[0] = pmother.Px();
362   p2[1] = pmother.Py();
363   p2[2] = pmother.Pz();
364   //
365   //
366   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
367   AliHelix mhelix(x2,p2,sign);
368   
369   //
370   //
371   //
372   //find intersection linear
373   //
374   Double_t distance1, distance2;
375   Double_t phase[2][2],radius[2];
376   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
377   Double_t delta1=10000,delta2=10000;    
378   if (points>0){
379     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
380     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
381     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
382   }
383   else{
384     fInvMass=-1;
385     return;
386   }
387   if (points==2){    
388     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
389     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
390     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
391   }
392   distance1 = TMath::Min(delta1,delta2);
393   //
394   //find intersection parabolic
395   //
396   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
397   delta1=10000,delta2=10000;  
398   
399   if (points>0){
400     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
401   }
402   if (points==2){    
403     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
404   }
405   
406   distance2 = TMath::Min(delta1,delta2);
407   //
408   if (delta1<delta2){
409     //get V0 info
410     dhelix1.Evaluate(phase[0][0],fMCXr);
411     dhelix1.GetMomentum(phase[0][0],fMCPdr);
412     mhelix.GetMomentum(phase[0][1],fMCPm);
413     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
414     fMCRr = TMath::Sqrt(radius[0]);
415   }
416   else{
417     dhelix1.Evaluate(phase[1][0],fMCXr);
418     dhelix1.GetMomentum(phase[1][0],fMCPdr);
419     mhelix.GetMomentum(phase[1][1],fMCPm);
420     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
421     fMCRr = TMath::Sqrt(radius[1]);
422   }     
423   //            
424   //
425   fMCDist1 = TMath::Sqrt(distance1);
426   fMCDist2 = TMath::Sqrt(distance2);      
427   //
428   //
429   // 
430   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
431   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
432   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
433   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
434   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
435   vnorm2 = TMath::Sqrt(vnorm2);
436   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
437   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
438   pnorm2 = TMath::Sqrt(pnorm2);
439   //
440   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
441   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
442   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
443   Float_t mass1 = fMCd.fMass;
444   Float_t mass2 = fMCm.fMass;
445
446   
447   Float_t e1    = TMath::Sqrt(mass1*mass1+
448                               fMCPdr[0]*fMCPdr[0]+
449                               fMCPdr[1]*fMCPdr[1]+
450                               fMCPdr[2]*fMCPdr[2]);
451   Float_t e2    = TMath::Sqrt(mass2*mass2+
452                               fMCPm[0]*fMCPm[0]+
453                               fMCPm[1]*fMCPm[1]+
454                               fMCPm[2]*fMCPm[2]);
455   
456   fInvMass =  
457     (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
458     (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
459     (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
460   
461   fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
462
463     
464 }
465
466
467
468
469
470 /////////////////////////////////////////////////////////////////////////////////
471 void AliGenKinkInfo::Update()
472 {
473   fMCPd[0] = fMCd.fParticle.Px();
474   fMCPd[1] = fMCd.fParticle.Py();
475   fMCPd[2] = fMCd.fParticle.Pz();
476   fMCPd[3] = fMCd.fParticle.P();
477   //
478   fMCX[0]  = fMCd.fParticle.Vx();
479   fMCX[1]  = fMCd.fParticle.Vy();
480   fMCX[2]  = fMCd.fParticle.Vz();
481   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
482   //
483   fPdg[0]    = fMCd.fParticle.GetPdgCode();
484   fPdg[1]    = fMCm.fParticle.GetPdgCode();
485   //
486   fLab[0]    = fMCd.fParticle.GetUniqueID();
487   fLab[1]    = fMCm.fParticle.GetUniqueID();
488   //
489   //
490   //
491   Double_t x1[3],p1[3];
492   TParticle & pdaughter = fMCd.fParticle;
493   x1[0] = pdaughter.Vx();      
494   x1[1] = pdaughter.Vy();
495   x1[2] = pdaughter.Vz();
496   p1[0] = pdaughter.Px();
497   p1[1] = pdaughter.Py();
498   p1[2] = pdaughter.Pz();
499   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
500   AliHelix dhelix1(x1,p1,sign);
501   //
502   //
503   Double_t x2[3],p2[3];            
504   //
505   TParticle & pmother = fMCm.fParticle;
506   x2[0] = pmother.Vx();      
507   x2[1] = pmother.Vy();      
508   x2[2] = pmother.Vz();      
509   p2[0] = pmother.Px();
510   p2[1] = pmother.Py();
511   p2[2] = pmother.Pz();
512   //
513   AliTrackReference & pdecay = fMCm.fTRdecay;
514   x2[0] = pdecay.X();      
515   x2[1] = pdecay.Y();      
516   x2[2] = pdecay.Z();      
517   p2[0] = pdecay.Px();
518   p2[1] = pdecay.Py();
519   p2[2] = pdecay.Pz();
520   //
521   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
522   AliHelix mhelix(x2,p2,sign);
523   
524   //
525   //
526   //
527   //find intersection linear
528   //
529   Double_t distance1, distance2;
530   Double_t phase[2][2],radius[2];
531   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
532   Double_t delta1=10000,delta2=10000;    
533   if (points>0){
534     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
535     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
536     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
537   }
538   if (points==2){    
539     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
540     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
541     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
542   }
543   distance1 = TMath::Min(delta1,delta2);
544   //
545   //find intersection parabolic
546   //
547   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
548   delta1=10000,delta2=10000;  
549   
550   if (points>0){
551     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
552   }
553   if (points==2){    
554     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
555   }
556   
557   distance2 = TMath::Min(delta1,delta2);
558   //
559   if (delta1<delta2){
560     //get V0 info
561     dhelix1.Evaluate(phase[0][0],fMCXr);
562     dhelix1.GetMomentum(phase[0][0],fMCPdr);
563     mhelix.GetMomentum(phase[0][1],fMCPm);
564     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
565     fMCRr = TMath::Sqrt(radius[0]);
566   }
567   else{
568     dhelix1.Evaluate(phase[1][0],fMCXr);
569     dhelix1.GetMomentum(phase[1][0],fMCPdr);
570     mhelix.GetMomentum(phase[1][1],fMCPm);
571     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
572     fMCRr = TMath::Sqrt(radius[1]);
573   }     
574   //            
575   //
576   fMCDist1 = TMath::Sqrt(distance1);
577   fMCDist2 = TMath::Sqrt(distance2);      
578     
579 }
580
581
582 Float_t AliGenKinkInfo::GetQt(){
583   //
584   //
585   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
586   return TMath::Sin(fMCAngle[2])*momentumd;
587 }
588
589
590
591   
592 ////////////////////////////////////////////////////////////////////////
593
594 ////////////////////////////////////////////////////////////////////////
595 //
596 // End of implementation of the class AliMCInfo
597 //
598 ////////////////////////////////////////////////////////////////////////
599
600
601
602 ////////////////////////////////////////////////////////////////////////
603 digitRow::digitRow()
604 {
605   Reset();
606 }
607 ////////////////////////////////////////////////////////////////////////
608 digitRow & digitRow::operator=(const digitRow &digOld)
609 {
610   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
611   return (*this);
612 }
613 ////////////////////////////////////////////////////////////////////////
614 void digitRow::SetRow(Int_t row) 
615 {
616   if (row >= 8*kgRowBytes) {
617     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
618     return;
619   }
620   Int_t iC = row/8;
621   Int_t iB = row%8;
622   SETBIT(fDig[iC],iB);
623 }
624
625 ////////////////////////////////////////////////////////////////////////
626 Bool_t digitRow::TestRow(Int_t row)
627 {
628 //
629 // return kTRUE if row is on
630 //
631   Int_t iC = row/8;
632   Int_t iB = row%8;
633   return TESTBIT(fDig[iC],iB);
634 }
635 ////////////////////////////////////////////////////////////////////////
636 Int_t digitRow::RowsOn(Int_t upto)
637 {
638 //
639 // returns number of rows with a digit  
640 // count only rows less equal row number upto
641 //
642   Int_t total = 0;
643   for (Int_t i = 0; i<kgRowBytes; i++) {
644     for (Int_t j = 0; j < 8; j++) {
645       if (i*8+j > upto) return total;
646       if (TESTBIT(fDig[i],j))  total++;
647     }
648   }
649   return total;
650 }
651 ////////////////////////////////////////////////////////////////////////
652 void digitRow::Reset()
653 {
654 //
655 // resets all rows to zero
656 //
657   for (Int_t i = 0; i<kgRowBytes; i++) {
658     fDig[i] <<= 8;
659   }
660 }
661 ////////////////////////////////////////////////////////////////////////
662 Int_t digitRow::Last()
663 {
664 //
665 // returns the last row number with a digit
666 // returns -1 if now digits 
667 //
668   for (Int_t i = kgRowBytes-1; i>=0; i--) {
669     for (Int_t j = 7; j >= 0; j--) {
670       if TESTBIT(fDig[i],j) return i*8+j;
671     }
672   }
673   return -1;
674 }
675 ////////////////////////////////////////////////////////////////////////
676 Int_t digitRow::First()
677 {
678 //
679 // returns the first row number with a digit
680 // returns -1 if now digits 
681 //
682   for (Int_t i = 0; i<kgRowBytes; i++) {
683     for (Int_t j = 0; j < 8; j++) {
684       if (TESTBIT(fDig[i],j)) return i*8+j;
685     }
686   }
687   return -1;
688 }
689
690 ////////////////////////////////////////////////////////////////////////
691 //
692 // end of implementation of a class digitRow
693 //
694 ////////////////////////////////////////////////////////////////////////
695   
696 ////////////////////////////////////////////////////////////////////////
697 AliGenInfoMaker::AliGenInfoMaker()
698 {
699   Reset();
700 }
701
702 ////////////////////////////////////////////////////////////////////////
703 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
704                                    Int_t nEvents, Int_t firstEvent)
705 {
706   Reset();
707   fFirstEventNr = firstEvent;
708   fEventNr = firstEvent;
709   fNEvents = nEvents;
710   //   fFnRes = fnRes;
711   sprintf(fFnRes,"%s",fnRes);
712   //
713   fLoader = AliRunLoader::Open(fnGalice);
714   if (gAlice){
715     delete gAlice->GetRunLoader();
716     delete gAlice;
717     gAlice = 0x0;
718   }
719   if (fLoader->LoadgAlice()){
720     cerr<<"Error occured while l"<<endl;
721   }
722   Int_t nall = fLoader->GetNumberOfEvents();
723   if (nEvents==0) {
724     nEvents =nall;
725     fNEvents=nall;
726     fFirstEventNr=0;
727   }    
728
729   if (nall<=0){
730     cerr<<"no events available"<<endl;
731     fEventNr = 0;
732     return;
733   }
734   if (firstEvent+nEvents>nall) {
735     fEventNr = nall-firstEvent;
736     cerr<<"restricted number of events availaible"<<endl;
737   }
738   AliMagF * magf = gAlice->Field();
739   AliTracker::SetFieldMap(magf);
740 }
741
742
743 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
744 {
745   // 
746   if (i<fNParticles) {
747     if (fGenInfo[i]) return  fGenInfo[i];
748     fGenInfo[i] = new AliMCInfo;  
749     fNInfos++;
750     return fGenInfo[i];
751   }
752   else 
753     return 0;  
754 }
755
756 ////////////////////////////////////////////////////////////////////////
757 void AliGenInfoMaker::Reset()
758 {
759   fEventNr = 0;
760   fNEvents = 0;
761   fTreeGenTracks = 0;
762   fFileGenTracks = 0;
763   fGenInfo = 0;
764   fNInfos  = 0;
765   //
766   //
767   fDebug = 0;
768   fVPrim[0] = -1000.;
769   fVPrim[1] = -1000.;
770   fVPrim[2] = -1000.;
771   fParamTPC = 0;
772 }
773 ////////////////////////////////////////////////////////////////////////
774 AliGenInfoMaker::~AliGenInfoMaker()
775 {
776   
777   if (fLoader){
778     fLoader->UnloadgAlice();
779     gAlice = 0;
780     delete fLoader;
781   }
782 }
783
784 Int_t  AliGenInfoMaker::SetIO()
785 {
786   //
787   // 
788   CreateTreeGenTracks();
789   if (!fTreeGenTracks) return 1;
790   //  AliTracker::SetFieldFactor(); 
791  
792   fParamTPC = GetTPCParam();
793   //
794   return 0;
795 }
796
797 ////////////////////////////////////////////////////////////////////////
798 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
799 {
800   //
801   // 
802   // SET INPUT
803   fLoader->SetEventNumber(eventNr);
804   //
805   fLoader->LoadHeader();
806   fLoader->LoadKinematics();  
807   fStack = fLoader->Stack();
808   //
809   fLoader->LoadTrackRefs();
810   fLoader->LoadHits();
811   fTreeTR = fLoader->TreeTR();
812   //
813   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
814   tpcl->LoadDigits();
815   fTreeD = tpcl->TreeD();  
816   return 0;
817 }
818
819 Int_t AliGenInfoMaker::CloseIOEvent()
820 {
821   fLoader->UnloadHeader();
822   fLoader->UnloadKinematics();
823   fLoader->UnloadTrackRefs();
824   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
825   tpcl->UnloadDigits();
826   return 0;
827 }
828
829 Int_t AliGenInfoMaker::CloseIO()
830 {
831   fLoader->UnloadgAlice();
832   return 0;
833 }
834
835
836
837 ////////////////////////////////////////////////////////////////////////
838 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
839 {
840   fNEvents = nEvents;
841   fFirstEventNr = firstEventNr;
842   return Exec();
843 }
844
845 ////////////////////////////////////////////////////////////////////////
846 Int_t AliGenInfoMaker::Exec()  
847 {
848   TStopwatch timer;
849   timer.Start();
850   Int_t status =SetIO();
851   if (status>0) return status;
852   //
853
854   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
855        fEventNr++) {
856     SetIO(fEventNr);
857     fNParticles = fStack->GetNtrack();
858     //
859     fGenInfo = new AliMCInfo*[fNParticles];
860     for (UInt_t i = 0; i<fNParticles; i++) {
861       fGenInfo[i]=0; 
862     }
863     //
864     cout<<"Start to process event "<<fEventNr<<endl;
865     cout<<"\tfNParticles = "<<fNParticles<<endl;
866     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
867     if (TreeTRLoop()>0) return 1;
868     //
869     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
870     if (TreeDLoop()>0) return 1;
871     //
872     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
873     if (TreeKLoop()>0) return 1;
874     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
875     //
876     if (BuildKinkInfo()>0) return 1;
877     if (BuildV0Info()>0) return 1;
878     //if (BuildHitLines()>0) return 1;
879     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
880     //
881     for (UInt_t i = 0; i<fNParticles; i++) {
882       if (fGenInfo[i]) delete fGenInfo[i]; 
883     }
884     delete []fGenInfo;
885     CloseIOEvent();
886   }
887   //
888   CloseIO();
889   CloseOutputFile();
890
891   cerr<<"Exec finished"<<endl;
892
893   timer.Stop();
894   timer.Print();
895   return 0;
896 }
897 ////////////////////////////////////////////////////////////////////////
898 void AliGenInfoMaker::CreateTreeGenTracks() 
899 {
900   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
901   if (!fFileGenTracks) {
902     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
903     return;
904   }
905   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
906   AliMCInfo * info = new AliMCInfo;
907   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
908   delete info; 
909   //
910   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
911   fTreeKinks = new TTree("genKinksTree","genKinksTree");  
912   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
913   delete kinkinfo;
914   //
915   AliGenV0Info *v0info = new AliGenV0Info;
916   fTreeV0 = new TTree("genV0Tree","genV0Tree");  
917   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
918   delete v0info;
919   //
920   //
921   AliPointsMI * points0 = new AliPointsMI;  
922   AliPointsMI * points1 = new AliPointsMI;  
923   AliPointsMI * points2 = new AliPointsMI;  
924   fTreeHitLines = new TTree("HitLines","HitLines");
925   fTreeHitLines->Branch("TPC.","AliPointsMI",&points0,32000,10);
926   fTreeHitLines->Branch("TRD.","AliPointsMI",&points1,32000,10);
927   fTreeHitLines->Branch("ITS.","AliPointsMI",&points2,32000,10);
928   Int_t label=0;
929   fTreeHitLines->Branch("Label",&label,"label/I");
930   //
931   fTreeGenTracks->AutoSave();
932   fTreeKinks->AutoSave();
933   fTreeV0->AutoSave();
934   fTreeHitLines->AutoSave();
935 }
936 ////////////////////////////////////////////////////////////////////////
937 void AliGenInfoMaker::CloseOutputFile() 
938 {
939   if (!fFileGenTracks) {
940     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
941     return;
942   }
943   fFileGenTracks->cd();
944   fTreeGenTracks->Write();  
945   delete fTreeGenTracks;
946   fTreeKinks->Write();
947   delete fTreeKinks;
948   fTreeV0->Write();
949   delete fTreeV0;
950
951   fFileGenTracks->Close();
952   delete fFileGenTracks;
953   return;
954 }
955
956 ////////////////////////////////////////////////////////////////////////
957 Int_t AliGenInfoMaker::TreeKLoop()
958 {
959 //
960 // open the file with treeK
961 // loop over all entries there and save information about some tracks
962 //
963
964   AliStack * stack = fStack;
965   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
966   
967   if (fDebug > 0) {
968     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
969         <<fEventNr<<endl;
970   }  
971   Int_t  ipdg = 0;
972   TParticlePDG *ppdg = 0;
973   // not all generators give primary vertex position. Take the vertex
974   // of the particle 0 as primary vertex.
975   TDatabasePDG  pdg; //get pdg table  
976   //thank you very much root for this
977   TBranch * br = fTreeGenTracks->GetBranch("MC");
978   TParticle *particle = stack->ParticleFromTreeK(0);
979   fVPrim[0] = particle->Vx();
980   fVPrim[1] = particle->Vy();
981   fVPrim[2] = particle->Vz();
982   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
983     // load only particles with TR
984     AliMCInfo * info = GetInfo(iParticle);
985     if (!info) continue;
986     //////////////////////////////////////////////////////////////////////
987     info->fLabel = iParticle;
988     //
989     info->fParticle = *(stack->Particle(iParticle));
990     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
991     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
992     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; 
993     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
994                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
995     //
996     //
997     ipdg = info->fParticle.GetPdgCode();
998     info->fPdg = ipdg;
999     ppdg = pdg.GetParticle(ipdg);          
1000     info->fEventNr = fEventNr;
1001     info->Update();
1002     //////////////////////////////////////////////////////////////////////    
1003     br->SetAddress(&info);    
1004     fTreeGenTracks->Fill();    
1005   }
1006   fTreeGenTracks->AutoSave();
1007   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1008   return 0;
1009 }
1010
1011
1012 ////////////////////////////////////////////////////////////////////////////////////
1013 Int_t  AliGenInfoMaker::BuildKinkInfo()
1014 {
1015   //
1016   // Build tree with Kink Information
1017   //
1018   AliStack * stack = fStack;
1019   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1020   
1021   if (fDebug > 0) {
1022     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1023         <<fEventNr<<endl;
1024   }  
1025   //  Int_t  ipdg = 0;
1026   //TParticlePDG *ppdg = 0;
1027   // not all generators give primary vertex position. Take the vertex
1028   // of the particle 0 as primary vertex.
1029   TDatabasePDG  pdg; //get pdg table  
1030   //thank you very much root for this
1031   TBranch * br = fTreeKinks->GetBranch("MC");
1032   //
1033   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
1034   //
1035   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1036     // load only particles with TR
1037     AliMCInfo * info = GetInfo(iParticle);
1038     if (!info) continue;
1039     if (info->fCharge==0) continue;  
1040     if (info->fTRdecay.P()<0.13) continue;  //momenta cut 
1041     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
1042     TParticle & particle = info->fParticle;
1043     Int_t first = particle.GetDaughter(0);
1044     if (first ==0) continue;
1045
1046     Int_t last = particle.GetDaughter(1);
1047     if (last ==0) last=first;
1048     AliMCInfo * info2 =  0;
1049     AliMCInfo * dinfo =  0;
1050     
1051     for (Int_t id2=first;id2<=last;id2++){
1052       info2 = GetInfo(id2);
1053       if (!info2) continue;
1054       TParticle &p2 = info2->fParticle;
1055       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
1056       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
1057       if (!(p2.GetPDG())) continue;
1058       dinfo =info2;
1059      
1060       kinkinfo->fMCm = (*info);
1061       kinkinfo->fMCd = (*dinfo);
1062       if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
1063       kinkinfo->Update();
1064       br->SetAddress(&kinkinfo);    
1065       fTreeKinks->Fill();
1066     }
1067     /*
1068     if (dinfo){
1069       kinkinfo->fMCm = (*info);
1070       kinkinfo->fMCd = (*dinfo);
1071       kinkinfo->Update();
1072       br->SetAddress(&kinkinfo);    
1073       fTreeKinks->Fill();
1074     }
1075     */
1076   }
1077   fTreeGenTracks->AutoSave();
1078   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1079   return 0;
1080 }
1081
1082
1083 ////////////////////////////////////////////////////////////////////////////////////
1084 Int_t  AliGenInfoMaker::BuildV0Info()
1085 {
1086   //
1087   // Build tree with V0 Information
1088   //
1089   AliStack * stack = fStack;
1090   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1091   
1092   if (fDebug > 0) {
1093     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1094         <<fEventNr<<endl;
1095   }  
1096   //  Int_t  ipdg = 0;
1097   //TParticlePDG *ppdg = 0;
1098   // not all generators give primary vertex position. Take the vertex
1099   // of the particle 0 as primary vertex.
1100   TDatabasePDG  pdg; //get pdg table  
1101   //thank you very much root for this
1102   TBranch * br = fTreeV0->GetBranch("MC");
1103   //
1104   AliGenV0Info * v0info = new AliGenV0Info;
1105   //
1106   //
1107   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1108     // load only particles with TR
1109     AliMCInfo * info = GetInfo(iParticle);
1110     if (!info) continue;
1111     if (info->fCharge==0) continue;  
1112     //
1113     //
1114     TParticle & particle = info->fParticle;
1115     Int_t mother = particle.GetMother(0);
1116     if (mother <=0) continue;
1117     //
1118     TParticle * motherparticle = stack->Particle(mother);
1119     if (!motherparticle) continue;
1120     //
1121     Int_t last = motherparticle->GetDaughter(1);
1122     if (last==(int)iParticle) continue;
1123     AliMCInfo * info2 =  info;
1124     AliMCInfo * dinfo =  GetInfo(last);
1125     if (!dinfo) continue;
1126     if (!dinfo->fParticle.GetPDG()) continue;
1127     if (!info2->fParticle.GetPDG()) continue;
1128     if (dinfo){
1129       v0info->fMCm = (*info);
1130       v0info->fMCd = (*dinfo);
1131       v0info->Update(fVPrim);
1132       br->SetAddress(&v0info);    
1133       fTreeV0->Fill();
1134     }
1135   }
1136   fTreeV0->AutoSave();
1137   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1138   return 0;
1139 }
1140
1141
1142
1143
1144 ////////////////////////////////////////////////////////////////////////
1145 Int_t AliGenInfoMaker::BuildHitLines()
1146 {
1147
1148 //
1149 // open the file with treeK
1150 // loop over all entries there and save information about some tracks
1151 //
1152
1153   AliStack * stack = fStack;
1154   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1155   
1156   if (fDebug > 0) {
1157     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1158         <<fEventNr<<endl;
1159   }  
1160   Int_t  ipdg = 0;
1161   // TParticlePDG *ppdg = 0;
1162   // not all generators give primary vertex position. Take the vertex
1163   // of the particle 0 as primary vertex.
1164   TDatabasePDG  pdg; //get pdg table  
1165   //thank you very much root for this
1166   AliPointsMI *tpcp = new AliPointsMI;
1167   AliPointsMI *trdp = new AliPointsMI;
1168   AliPointsMI *itsp = new AliPointsMI;
1169   Int_t label =0;
1170   //
1171   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1172   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");  
1173   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1174   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1175   brlabel->SetAddress(&label);
1176   brtpc->SetAddress(&tpcp);
1177   brtrd->SetAddress(&trdp);
1178   brits->SetAddress(&itsp);
1179   //
1180   AliDetector *dtpc = gAlice->GetDetector("TPC");
1181   AliDetector *dtrd = gAlice->GetDetector("TRD");
1182   AliDetector *dits = gAlice->GetDetector("ITS");
1183  
1184   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1185     // load only particles with TR
1186     AliMCInfo * info = GetInfo(iParticle);
1187     if (!info) continue;
1188     Int_t primpart = info->fPrimPart;
1189     ipdg = info->fParticle.GetPdgCode();
1190     label = iParticle;
1191     //
1192     gAlice->ResetHits();
1193     tpcp->Reset();
1194     itsp->Reset();
1195     trdp->Reset();
1196     tpcp->fLabel1 = ipdg;
1197     trdp->fLabel1 = ipdg;
1198     itsp->fLabel1 = ipdg;
1199     if (dtpc->TreeH()->GetEvent(primpart)){
1200       dtpc->LoadPoints(primpart);
1201       tpcp->Reset(dtpc,iParticle);
1202     }
1203     if (dtrd->TreeH()->GetEvent(primpart)){
1204       dtrd->LoadPoints(primpart);
1205       trdp->Reset(dtrd,iParticle);
1206     }
1207     if (dits->TreeH()->GetEvent(primpart)){
1208       dits->LoadPoints(primpart);
1209       itsp->Reset(dits,iParticle);
1210     }    
1211     //    
1212     fTreeHitLines->Fill();
1213     dtpc->ResetPoints();
1214     dtrd->ResetPoints();
1215     dits->ResetPoints();
1216   }
1217   delete tpcp;
1218   delete trdp;
1219   delete itsp;
1220   fTreeHitLines->AutoSave();
1221   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1222   return 0;
1223 }
1224
1225
1226 ////////////////////////////////////////////////////////////////////////
1227 Int_t AliGenInfoMaker::TreeDLoop()
1228 {
1229   //
1230   // open the file with treeD
1231   // loop over all entries there and save information about some tracks
1232   //
1233   
1234   Int_t nInnerSector = fParamTPC->GetNInnerSector();
1235   Int_t rowShift = 0;
1236   Int_t zero=fParamTPC->GetZeroSup()+6;  
1237   //
1238   //
1239   AliSimDigits digitsAddress, *digits=&digitsAddress;
1240   fTreeD->GetBranch("Segment")->SetAddress(&digits);
1241   
1242   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1243   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1244   for (Int_t i=0; i<sectorsByRows; i++) {
1245     if (!fTreeD->GetEvent(i)) continue;
1246     Int_t sec,row;
1247     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1248     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
1249     // here I expect that upper sectors follow lower sectors
1250     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1251     //
1252     digits->ExpandTrackBuffer();
1253     digits->First();        
1254     do {
1255       Int_t iRow=digits->CurrentRow();
1256       Int_t iColumn=digits->CurrentColumn();
1257       Short_t digitValue = digits->CurrentDigit();
1258       if (digitValue >= zero) {
1259         Int_t label;
1260         for (Int_t j = 0; j<3; j++) {
1261           //      label = digits->GetTrackID(iRow,iColumn,j);
1262           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;       
1263           if (label >= (Int_t)fNParticles) { //don't label from bakground event
1264             continue;
1265           }
1266           if (label >= 0 && label <= (Int_t)fNParticles) {
1267             if (fDebug > 6 ) {
1268               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1269                   <<sec<<" "
1270                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1271                   <<" "<<row<<endl;
1272             }   
1273             AliMCInfo * info = GetInfo(label);
1274             if (info){
1275               info->fTPCRow.SetRow(row+rowShift);
1276             }
1277           }
1278         }
1279       }
1280     } while (digits->Next());
1281   }
1282   
1283   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;  
1284   return 0;
1285 }
1286
1287
1288 ////////////////////////////////////////////////////////////////////////
1289 Int_t AliGenInfoMaker::TreeTRLoop()
1290 {
1291   //
1292   // loop over TrackReferences and store the first one for each track
1293   //  
1294   TTree * treeTR = fTreeTR;
1295   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1296   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1297   //
1298   //
1299   //track references for TPC
1300   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
1301   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
1302   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
1303   TClonesArray* TOFArrayTR = new TClonesArray("AliTrackReference");
1304   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
1305   //
1306   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
1307   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
1308   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
1309   if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
1310   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
1311   //
1312   //
1313   //
1314   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1315     treeTR->GetEntry(iPrimPart);
1316     //
1317     // Loop over TPC references
1318     //
1319     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
1320       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
1321       //
1322       if (trackRef->P()<fgTPCPtCut) continue;
1323       Int_t label = trackRef->GetTrack();      
1324       AliMCInfo * info = GetInfo(label);
1325       if (!info) info = MakeInfo(label);
1326       if (!info) continue;
1327       info->fPrimPart =  iPrimPart;
1328       TClonesArray & arr = *(info->fTPCReferences);
1329       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1330     }
1331     //
1332     // Loop over ITS references
1333     //
1334     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
1335       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);            
1336       //
1337       if (trackRef->P()<fgTPCPtCut) continue;
1338       Int_t label = trackRef->GetTrack();      
1339       AliMCInfo * info = GetInfo(label);
1340       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
1341       if (!info) info = MakeInfo(label);
1342       if (!info) continue;
1343       info->fPrimPart =  iPrimPart;
1344       TClonesArray & arr = *(info->fITSReferences);
1345       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1346     }
1347     //
1348     // Loop over TRD references
1349     //
1350     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
1351       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);            
1352       //
1353       if (trackRef->P()<fgTPCPtCut) continue;
1354       Int_t label = trackRef->GetTrack();      
1355       AliMCInfo * info = GetInfo(label);
1356       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
1357       if (!info) info = MakeInfo(label);
1358       if (!info) continue;
1359       info->fPrimPart =  iPrimPart;
1360       TClonesArray & arr = *(info->fTRDReferences);
1361       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1362     }
1363     //
1364     // Loop over TOF references
1365     //
1366     for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
1367       AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);            
1368       Int_t label = trackRef->GetTrack();      
1369       AliMCInfo * info = GetInfo(label);
1370       if (!info){
1371         if (trackRef->Pt()<fgTPCPtCut) continue;
1372         if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
1373       }
1374       if (!info) info = MakeInfo(label);
1375       if (!info) continue;
1376       info->fPrimPart =  iPrimPart;
1377       TClonesArray & arr = *(info->fTOFReferences);
1378       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1379     }
1380     //
1381     // get dacay position
1382     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
1383       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
1384       //
1385       Int_t label = trackRef->GetTrack();
1386       AliMCInfo * info = GetInfo(label);
1387       if (!info) continue;
1388       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
1389       //      if (TMath::Abs(trackRef.X());
1390       info->fTRdecay = *trackRef;      
1391     }
1392   }
1393   //
1394   TPCArrayTR->Delete();
1395   delete  TPCArrayTR;
1396   TRDArrayTR->Delete();
1397   delete  TRDArrayTR;
1398   TOFArrayTR->Delete();
1399   delete  TOFArrayTR;
1400
1401   ITSArrayTR->Delete();
1402   delete  ITSArrayTR;
1403   RunArrayTR->Delete();
1404   delete  RunArrayTR;
1405   //
1406   return 0;
1407 }
1408
1409 ////////////////////////////////////////////////////////////////////////
1410 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
1411                                     AliTPCParam *paramTPC) {
1412
1413   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1414   Int_t index[4];
1415   paramTPC->Transform0to1(x,index);
1416   paramTPC->Transform1to2(x,index);
1417   return x[0];
1418 }
1419 ////////////////////////////////////////////////////////////////////////
1420
1421
1422
1423 TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection, 
1424                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
1425 {
1426   //
1427   Double_t* bins = CreateLogBins(nbins, minx, maxx);
1428   Int_t nBinsRes = 30;
1429   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
1430   char cut[1000];
1431   sprintf(cut,"%s&&%s",selection,quality);
1432   char expression[1000];
1433   sprintf(expression,"%s:%s>>hRes2",chy,chx);
1434   fTree->Draw(expression, cut, "groff");
1435   TH1F* hMean=0;
1436   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1437   AliLabelAxes(hRes, chx, chy);
1438   //
1439   delete hRes2;
1440   delete[] bins;
1441   fRes  = hRes;
1442   fMean = hMean;
1443   return hRes;
1444 }
1445
1446 TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, 
1447                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
1448 {
1449   //
1450   Double_t* bins = CreateLogBins(nbins, minx, maxx);
1451   Int_t nBinsRes = 30;
1452   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
1453   char cut[1000];
1454   sprintf(cut,"%s&&%s",selection,quality);
1455   char expression[1000];
1456   sprintf(expression,"%s:%s>>hRes2",chy,chx);
1457   fTree->Draw(expression, cut, "groff");
1458   TH1F* hMean=0;  
1459   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1460   AliLabelAxes(hRes, chx, chy);
1461   //
1462   delete hRes2;
1463   delete[] bins;
1464   fRes  = hRes;
1465   fMean = hMean;
1466   return hRes;
1467 }
1468
1469 ///////////////////////////////////////////////////////////////////////////////////
1470 ///////////////////////////////////////////////////////////////////////////////////
1471 TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality, 
1472                               Int_t nbins, Float_t min, Float_t max)
1473 {
1474   //
1475   //
1476   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
1477   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
1478   char inputGen[1000];  
1479   sprintf(inputGen,"%s>>hGen", variable);
1480   fTree->Draw(inputGen, selection, "groff");
1481   char selectionRec[256];
1482   sprintf(selectionRec, "%s && %s", selection, quality);
1483   char inputRec[1000];  
1484   sprintf(inputRec,"%s>>hRec", variable);
1485   fTree->Draw(inputRec, selectionRec, "groff");
1486   //
1487   TH1F* hEff = CreateEffHisto(hGen, hRec);
1488   AliLabelAxes(hEff, variable, "#epsilon [%]");
1489   fRes = hEff;
1490   delete hRec;
1491   delete hGen;
1492   return hEff;
1493 }
1494
1495
1496
1497 ///////////////////////////////////////////////////////////////////////////////////
1498 ///////////////////////////////////////////////////////////////////////////////////
1499 TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality, 
1500                               Int_t nbins, Float_t min, Float_t max)
1501 {
1502   //
1503   //
1504   Double_t* bins = CreateLogBins(nbins, min, max);
1505   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
1506   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
1507   char inputGen[1000];  
1508   sprintf(inputGen,"%s>>hGen", variable);
1509   fTree->Draw(inputGen, selection, "groff");
1510   char selectionRec[256];
1511   sprintf(selectionRec, "%s && %s", selection, quality);
1512   char inputRec[1000];  
1513   sprintf(inputRec,"%s>>hRec", variable);
1514   fTree->Draw(inputRec, selectionRec, "groff");
1515   //
1516   TH1F* hEff = CreateEffHisto(hGen, hRec);
1517   AliLabelAxes(hEff, variable, "#epsilon [%]");
1518   fRes = hEff;
1519   delete hRec;
1520   delete hGen;
1521   delete[] bins;
1522   return hEff;
1523 }
1524
1525
1526 ///////////////////////////////////////////////////////////////////////////////////
1527 ///////////////////////////////////////////////////////////////////////////////////
1528
1529 Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1530 {
1531   Double_t* bins = new Double_t[nBins+1];
1532   bins[0] = xMin;
1533   Double_t factor = pow(xMax/xMin, 1./nBins);
1534   for (Int_t i = 1; i <= nBins; i++)
1535     bins[i] = factor * bins[i-1];
1536   return bins;
1537 }
1538
1539
1540
1541
1542 void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1543 {
1544   //
1545   histo->GetXaxis()->SetTitle(xAxisTitle);
1546   histo->GetYaxis()->SetTitle(yAxisTitle);
1547 }
1548
1549
1550 TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1551 {
1552   //
1553   Int_t nBins = hGen->GetNbinsX();
1554   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1555   hEff->SetTitle("");
1556   hEff->SetStats(kFALSE);
1557   hEff->SetMinimum(0.);
1558   hEff->SetMaximum(110.);
1559   //
1560   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1561     Double_t nGen = hGen->GetBinContent(iBin);
1562     Double_t nRec = hRec->GetBinContent(iBin);
1563     if (nGen > 0) {
1564       Double_t eff = nRec/nGen;
1565       hEff->SetBinContent(iBin, 100. * eff);
1566       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
1567       //      if (error == 0) error = sqrt(0.1/nGen);
1568       //
1569       if (error == 0) error = 0.0001;
1570       hEff->SetBinError(iBin, 100. * error);
1571     } else {
1572       hEff->SetBinContent(iBin, 100. * 0.5);
1573       hEff->SetBinError(iBin, 100. * 0.5);
1574     }
1575   }
1576   return hEff;
1577 }
1578
1579
1580
1581 TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
1582                      Bool_t overflowBinFits)
1583 {
1584   TVirtualPad* currentPad = gPad;
1585   TAxis* axis = hRes2->GetXaxis();
1586   Int_t nBins = axis->GetNbins();
1587   TH1F* hRes, *hMean;
1588   if (axis->GetXbins()->GetSize()){
1589     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1590     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1591   }
1592   else{
1593     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1594     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1595
1596   }
1597   hRes->SetStats(false);
1598   hRes->SetOption("E");
1599   hRes->SetMinimum(0.);
1600   //
1601   hMean->SetStats(false);
1602   hMean->SetOption("E");
1603  
1604   // create the fit function
1605   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1606   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1607   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1608   
1609   fitFunc->SetLineWidth(2);
1610   fitFunc->SetFillStyle(0);
1611   // create canvas for fits
1612   TCanvas* canBinFits = NULL;
1613   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1614   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1615   Int_t ny = (nPads-1) / nx + 1;
1616   if (drawBinFits) {
1617     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1618     if (canBinFits) delete canBinFits;
1619     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1620     canBinFits->Divide(nx, ny);
1621   }
1622
1623   // loop over x bins and fit projection
1624   Int_t dBin = ((overflowBinFits) ? 1 : 0);
1625   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1626     if (drawBinFits) canBinFits->cd(bin + dBin);
1627     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1628     //    
1629     if (hBin->GetEntries() > 5) {
1630       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1631       hBin->Fit(fitFunc,"s");
1632       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1633
1634       if (sigma > 0.){
1635         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1636         hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
1637       }
1638       else{
1639         hRes->SetBinContent(bin, 0.);
1640         hMean->SetBinContent(bin,0);
1641       }
1642       hRes->SetBinError(bin, fitFunc->GetParError(2));
1643       hMean->SetBinError(bin, fitFunc->GetParError(1));
1644       
1645       //
1646       //
1647
1648     } else {
1649       hRes->SetBinContent(bin, 0.);
1650       hRes->SetBinError(bin, 0.);
1651       hMean->SetBinContent(bin, 0.);
1652       hMean->SetBinError(bin, 0.);
1653     }
1654     
1655
1656     if (drawBinFits) {
1657       char name[256];
1658       if (bin == 0) {
1659         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1660       } else if (bin == nBins+1) {
1661         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1662       } else {
1663         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1664                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1665       }
1666       canBinFits->cd(bin + dBin);
1667       hBin->SetTitle(name);
1668       hBin->SetStats(kTRUE);
1669       hBin->DrawCopy("E");
1670       canBinFits->Update();
1671       canBinFits->Modified();
1672       canBinFits->Update();
1673     }
1674     
1675     delete hBin;
1676   }
1677
1678   delete fitFunc;
1679   currentPad->cd();
1680   *phMean = hMean;
1681   return hRes;
1682 }
1683
1684
1685 void   AliComparisonDraw::DrawFriend2D(const char * chx, const char *chy, const char* selection, TTree * tfriend)
1686 {
1687   
1688 }
1689
1690
1691 void   AliComparisonDraw::GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
1692   //
1693   //
1694   if (!fPoints) fPoints = new TObjArray;
1695   if (tpoints->GetIndex()==0) tpoints->BuildIndex("Label","Label");
1696   AliPointsMI * points = new AliPointsMI;
1697   TBranch * br = tpoints->GetBranch(chpoints);
1698   br->SetAddress(&points);
1699   Int_t npoints = fTree->Draw(label,selection);
1700   Float_t xyz[30000];
1701   rmin*=rmin;
1702   for (Int_t i=0;i<npoints;i++){
1703     Int_t index = (Int_t)fTree->GetV1()[i];
1704     tpoints->GetEntryWithIndex(index,index);
1705     if (points->fN<2) continue;
1706     Int_t goodpoints=0;
1707     for (Int_t i=0;i<points->fN; i++){
1708       xyz[goodpoints*3]   = points->fX[i];
1709       xyz[goodpoints*3+1] = points->fY[i];
1710       xyz[goodpoints*3+2] = points->fZ[i];
1711       if ( points->fX[i]*points->fX[i]+points->fY[i]*points->fY[i]>rmin) goodpoints++;
1712     } 
1713     TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
1714     //    marker->SeWidth(1); 
1715     marker->SetMarkerColor(color);
1716     marker->SetMarkerStyle(1);
1717     fPoints->AddLast(marker);
1718   }
1719 }
1720
1721 void AliComparisonDraw::Draw3D(Int_t min, Int_t max){
1722   if (!fPoints) return;
1723   Int_t npoints = fPoints->GetEntriesFast();
1724   max = TMath::Min(npoints,max);
1725   min = TMath::Min(npoints,min);
1726
1727   for (Int_t i =min; i<max;i++){
1728     TObject * obj = fPoints->At(i);
1729     if (obj) obj->Draw();
1730   }
1731 }
1732
1733 void AliComparisonDraw:: SavePoints(const char* name)
1734 {
1735   char fname[25];
1736   sprintf(fname,"TH_%s.root",name);
1737   TFile f(fname,"new");
1738   fPoints->Write("Tracks",TObject::kSingleKey);
1739   f.Close();
1740   sprintf(fname,"TH%s.gif",name);
1741   fCanvas->SaveAs(fname);
1742 }
1743
1744 void   AliComparisonDraw::InitView(){
1745   //
1746   //
1747   AliRunLoader* rl = AliRunLoader::Open();
1748   rl->GetEvent(0); 
1749   rl->CdGAFile(); 
1750   //
1751   fCanvas = new TCanvas("cdisplay", "Cluster display",0,0,700,730);
1752   fView   = new TView(1);
1753   fView->SetRange(-330,-360,-330, 360,360, 330);
1754   fCanvas->Clear();
1755   fCanvas->SetFillColor(1);
1756   fCanvas->SetTheta(90.);
1757   fCanvas->SetPhi(0.);
1758   //
1759   TGeometry *geom =(TGeometry*)gDirectory->Get("AliceGeom");
1760   geom->Draw("same");
1761
1762 }