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