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