0214cf07db9f22966e5fa8c4067545dbb87253b5
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfo.cxx
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 Generate complex MC information - used for Comparison later on
23 How to use it?
24
25 .L 
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 "TROOT.h"
36 #include "Rtypes.h"
37 #include "TFile.h"
38 #include "TTree.h"
39 #include "TChain.h"
40 #include "TCut.h"
41 #include "TString.h"
42 #include "TBenchmark.h"
43 #include "TStopwatch.h"
44 #include "TParticle.h"
45 #include "TSystem.h"
46 #include "TTimer.h"
47 #include "TVector3.h"
48 #include "TH1F.h"
49 #include "TH2F.h"
50 #include "TCanvas.h"
51 #include "TPad.h"
52 #include "TF1.h"
53 #include "TView.h"
54 #include "TView3D.h"
55 #include "TGeometry.h"
56 #include "TPolyLine3D.h"
57
58 //ALIROOT includes
59 #include "AliRun.h"
60 #include "AliStack.h"
61 #include "AliSimDigits.h"
62 #include "AliTPCParam.h"
63 #include "AliTPC.h"
64 #include "AliTPCLoader.h"
65 #include "AliDetector.h"
66 #include "AliTrackReference.h"
67 #include "AliTPCParamSR.h"
68 #include "AliTracker.h"
69 #include "AliMagF.h"
70 #include "AliHelix.h"
71 #include "AliPoints.h"
72 #include "AliTrackPointArray.h"
73
74 #endif
75 #include "AliGenInfo.h" 
76 //
77 // 
78
79 ClassImp(AliTPCdigitRow);
80 ClassImp(AliMCInfo);
81 ClassImp(AliGenV0Info)
82 ClassImp(AliGenKinkInfo)
83 ClassImp(AliGenInfoMaker)
84
85
86
87 AliTPCParam * GetTPCParam(){
88   AliTPCParamSR * par = new AliTPCParamSR;
89   par->Update();
90   return par;
91 }
92
93
94 //_____________________________________________________________________________
95 Float_t TPCBetheBloch(Float_t bg)
96 {
97   //
98   // Bethe-Bloch energy loss formula
99   //
100   const Double_t kp1=0.76176e-1;
101   const Double_t kp2=10.632;
102   const Double_t kp3=0.13279e-4;
103   const Double_t kp4=1.8631;
104   const Double_t kp5=1.9479;
105
106   Double_t dbg = (Double_t) bg;
107
108   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
109
110   Double_t aa = TMath::Power(beta,kp4);
111   Double_t bb = TMath::Power(1./dbg,kp5);
112
113   bb=TMath::Log(kp3+bb);
114   
115   return ((Float_t)((kp2-aa-bb)*kp1/aa));
116 }
117
118
119
120
121
122 ////////////////////////////////////////////////////////////////////////
123 AliMCInfo::AliMCInfo()
124 {
125   fTPCReferences  = new TClonesArray("AliTrackReference",10);
126   fITSReferences  = new TClonesArray("AliTrackReference",10);
127   fTRDReferences  = new TClonesArray("AliTrackReference",10);
128   fTOFReferences  = new TClonesArray("AliTrackReference",10);
129   fTRdecay.SetTrack(-1);
130   fCharge = 0;
131 }
132
133 AliMCInfo::~AliMCInfo()
134 {
135   if (fTPCReferences) {
136     delete fTPCReferences;
137   }
138   if (fITSReferences){
139     delete fITSReferences;
140   }
141   if (fTRDReferences){    
142     delete fTRDReferences;  
143   }
144   if (fTOFReferences){    
145     delete fTOFReferences;  
146   }
147
148 }
149
150
151
152 void AliMCInfo::Update()
153 {
154   //
155   //
156   fMCtracks =1;
157   if (!fTPCReferences) {
158     fNTPCRef =0;
159     return;
160   }
161   Float_t direction=1;
162   //Float_t rlast=0;
163   fNTPCRef = fTPCReferences->GetEntriesFast();
164   fNITSRef = fITSReferences->GetEntriesFast();
165   fNTRDRef = fTRDReferences->GetEntriesFast();
166   fNTOFRef = fTOFReferences->GetEntriesFast();
167   
168   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
169     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
170     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
171     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
172     if (iref==0) direction = newdirection;
173     if ( newdirection*direction<0){
174       //changed direction
175       direction = newdirection;
176       fMCtracks+=1;
177     }
178     //rlast=r;                      
179   }
180   //
181   // decay info
182   fTPCdecay=kFALSE;
183   if (fTRdecay.GetTrack()>0){
184     fDecayCoord[0] = fTRdecay.X();
185     fDecayCoord[1] = fTRdecay.Y();
186     fDecayCoord[2] = fTRdecay.Z();
187     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
188       fTPCdecay=kTRUE;     
189     }
190     else{
191       fDecayCoord[0] = 0;
192       fDecayCoord[1] = 0;
193       fDecayCoord[2] = 0;
194     }
195   }
196   //
197   //
198   //digits information update
199   fRowsWithDigits    = fTPCRow.RowsOn();    
200   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
201   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
202   //
203   //
204   // calculate primary ionization per cm  
205   if (fParticle.GetPDG()){
206     fMass = fParticle.GetMass();  
207     fCharge = fParticle.GetPDG()->Charge();
208     if (fTPCReferences->GetEntriesFast()>0){
209       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
210     }
211     if (fMass>0){
212       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
213                               fTrackRef.Py()*fTrackRef.Py()+
214                               fTrackRef.Pz()*fTrackRef.Pz());    
215       if (p>0.001){
216         Float_t betagama = p /fMass;
217         fPrim = TPCBetheBloch(betagama);
218       }else fPrim=0;
219     }
220   }else{
221     fMass =0;
222     fPrim =0;
223   }  
224 }
225
226 /////////////////////////////////////////////////////////////////////////////////
227 /*
228 void AliGenV0Info::Update()
229 {
230   fMCPd[0] = fMCd.fParticle.Px();
231   fMCPd[1] = fMCd.fParticle.Py();
232   fMCPd[2] = fMCd.fParticle.Pz();
233   fMCPd[3] = fMCd.fParticle.P();
234   //
235   fMCX[0]  = fMCd.fParticle.Vx();
236   fMCX[1]  = fMCd.fParticle.Vy();
237   fMCX[2]  = fMCd.fParticle.Vz();
238   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
239   //
240   fPdg[0]    = fMCd.fParticle.GetPdgCode();
241   fPdg[1]    = fMCm.fParticle.GetPdgCode();
242   //
243   fLab[0]    = fMCd.fParticle.GetUniqueID();
244   fLab[1]    = fMCm.fParticle.GetUniqueID();
245   //  
246 }
247 */
248
249 void AliGenV0Info::Update(Float_t vertex[3])
250 {
251   fMCPd[0] = fMCd.fParticle.Px();
252   fMCPd[1] = fMCd.fParticle.Py();
253   fMCPd[2] = fMCd.fParticle.Pz();
254   fMCPd[3] = fMCd.fParticle.P();
255   //
256   fMCX[0]  = fMCd.fParticle.Vx();
257   fMCX[1]  = fMCd.fParticle.Vy();
258   fMCX[2]  = fMCd.fParticle.Vz();
259   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
260   //
261   fPdg[0]    = fMCd.fParticle.GetPdgCode();
262   fPdg[1]    = fMCm.fParticle.GetPdgCode();
263   //
264   fLab[0]    = fMCd.fParticle.GetUniqueID();
265   fLab[1]    = fMCm.fParticle.GetUniqueID();
266   //
267   //
268   //
269   Double_t x1[3],p1[3];
270   TParticle & pdaughter = fMCd.fParticle;
271   x1[0] = pdaughter.Vx();      
272   x1[1] = pdaughter.Vy();
273   x1[2] = pdaughter.Vz();
274   p1[0] = pdaughter.Px();
275   p1[1] = pdaughter.Py();
276   p1[2] = pdaughter.Pz();
277   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
278   AliHelix dhelix1(x1,p1,sign);
279   //
280   //
281   Double_t x2[3],p2[3];            
282   //
283   TParticle & pmother = fMCm.fParticle;
284   x2[0] = pmother.Vx();      
285   x2[1] = pmother.Vy();      
286   x2[2] = pmother.Vz();      
287   p2[0] = pmother.Px();
288   p2[1] = pmother.Py();
289   p2[2] = pmother.Pz();
290   //
291   //
292   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
293   AliHelix mhelix(x2,p2,sign);
294   
295   //
296   //
297   //
298   //find intersection linear
299   //
300   Double_t distance1, distance2;
301   Double_t phase[2][2],radius[2];
302   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
303   Double_t delta1=10000,delta2=10000;    
304   if (points>0){
305     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
306     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
307     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
308   }
309   else{
310     fInvMass=-1;
311     return;
312   }
313   if (points==2){    
314     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
315     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
316     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
317   }
318   distance1 = TMath::Min(delta1,delta2);
319   //
320   //find intersection parabolic
321   //
322   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
323   delta1=10000,delta2=10000;  
324   
325   if (points>0){
326     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
327   }
328   if (points==2){    
329     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
330   }
331   
332   distance2 = TMath::Min(delta1,delta2);
333   //
334   if (delta1<delta2){
335     //get V0 info
336     dhelix1.Evaluate(phase[0][0],fMCXr);
337     dhelix1.GetMomentum(phase[0][0],fMCPdr);
338     mhelix.GetMomentum(phase[0][1],fMCPm);
339     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
340     fMCRr = TMath::Sqrt(radius[0]);
341   }
342   else{
343     dhelix1.Evaluate(phase[1][0],fMCXr);
344     dhelix1.GetMomentum(phase[1][0],fMCPdr);
345     mhelix.GetMomentum(phase[1][1],fMCPm);
346     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
347     fMCRr = TMath::Sqrt(radius[1]);
348   }     
349   //            
350   //
351   fMCDist1 = TMath::Sqrt(distance1);
352   fMCDist2 = TMath::Sqrt(distance2);      
353   //
354   //
355   // 
356   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
357   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
358   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
359   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
360   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
361   vnorm2 = TMath::Sqrt(vnorm2);
362   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
363   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
364   pnorm2 = TMath::Sqrt(pnorm2);
365   //
366   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
367   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
368   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
369   Double_t mass1 = fMCd.fMass;
370   Double_t mass2 = fMCm.fMass;
371
372   
373   Double_t e1    = TMath::Sqrt(mass1*mass1+
374                               fMCPd[0]*fMCPd[0]+
375                               fMCPd[1]*fMCPd[1]+
376                               fMCPd[2]*fMCPd[2]);
377   Double_t e2    = TMath::Sqrt(mass2*mass2+
378                               fMCPm[0]*fMCPm[0]+
379                               fMCPm[1]*fMCPm[1]+
380                               fMCPm[2]*fMCPm[2]);
381   
382   fInvMass =  
383     (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
384     (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
385     (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
386   
387   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
388   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
389   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
390
391     
392 }
393
394
395
396 /////////////////////////////////////////////////////////////////////////////////
397 void AliGenKinkInfo::Update()
398 {
399   fMCPd[0] = fMCd.fParticle.Px();
400   fMCPd[1] = fMCd.fParticle.Py();
401   fMCPd[2] = fMCd.fParticle.Pz();
402   fMCPd[3] = fMCd.fParticle.P();
403   //
404   fMCX[0]  = fMCd.fParticle.Vx();
405   fMCX[1]  = fMCd.fParticle.Vy();
406   fMCX[2]  = fMCd.fParticle.Vz();
407   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
408   //
409   fPdg[0]    = fMCd.fParticle.GetPdgCode();
410   fPdg[1]    = fMCm.fParticle.GetPdgCode();
411   //
412   fLab[0]    = fMCd.fParticle.GetUniqueID();
413   fLab[1]    = fMCm.fParticle.GetUniqueID();
414   //
415   //
416   //
417   Double_t x1[3],p1[3];
418   TParticle & pdaughter = fMCd.fParticle;
419   x1[0] = pdaughter.Vx();      
420   x1[1] = pdaughter.Vy();
421   x1[2] = pdaughter.Vz();
422   p1[0] = pdaughter.Px();
423   p1[1] = pdaughter.Py();
424   p1[2] = pdaughter.Pz();
425   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
426   AliHelix dhelix1(x1,p1,sign);
427   //
428   //
429   Double_t x2[3],p2[3];            
430   //
431   TParticle & pmother = fMCm.fParticle;
432   x2[0] = pmother.Vx();      
433   x2[1] = pmother.Vy();      
434   x2[2] = pmother.Vz();      
435   p2[0] = pmother.Px();
436   p2[1] = pmother.Py();
437   p2[2] = pmother.Pz();
438   //
439   AliTrackReference & pdecay = fMCm.fTRdecay;
440   x2[0] = pdecay.X();      
441   x2[1] = pdecay.Y();      
442   x2[2] = pdecay.Z();      
443   p2[0] = pdecay.Px();
444   p2[1] = pdecay.Py();
445   p2[2] = pdecay.Pz();
446   //
447   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
448   AliHelix mhelix(x2,p2,sign);
449   
450   //
451   //
452   //
453   //find intersection linear
454   //
455   Double_t distance1, distance2;
456   Double_t phase[2][2],radius[2];
457   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
458   Double_t delta1=10000,delta2=10000;    
459   if (points>0){
460     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
461     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
462     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
463   }
464   if (points==2){    
465     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
466     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
467     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
468   }
469   distance1 = TMath::Min(delta1,delta2);
470   //
471   //find intersection parabolic
472   //
473   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
474   delta1=10000,delta2=10000;  
475   
476   if (points>0){
477     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
478   }
479   if (points==2){    
480     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
481   }
482   
483   distance2 = TMath::Min(delta1,delta2);
484   //
485   if (delta1<delta2){
486     //get V0 info
487     dhelix1.Evaluate(phase[0][0],fMCXr);
488     dhelix1.GetMomentum(phase[0][0],fMCPdr);
489     mhelix.GetMomentum(phase[0][1],fMCPm);
490     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
491     fMCRr = TMath::Sqrt(radius[0]);
492   }
493   else{
494     dhelix1.Evaluate(phase[1][0],fMCXr);
495     dhelix1.GetMomentum(phase[1][0],fMCPdr);
496     mhelix.GetMomentum(phase[1][1],fMCPm);
497     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
498     fMCRr = TMath::Sqrt(radius[1]);
499   }     
500   //            
501   //
502   fMCDist1 = TMath::Sqrt(distance1);
503   fMCDist2 = TMath::Sqrt(distance2);      
504     
505 }
506
507
508 Float_t AliGenKinkInfo::GetQt(){
509   //
510   //
511   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
512   return TMath::Sin(fMCAngle[2])*momentumd;
513 }
514
515
516
517   
518 ////////////////////////////////////////////////////////////////////////
519
520 ////////////////////////////////////////////////////////////////////////
521 //
522 // End of implementation of the class AliMCInfo
523 //
524 ////////////////////////////////////////////////////////////////////////
525
526
527
528 ////////////////////////////////////////////////////////////////////////
529 AliTPCdigitRow::AliTPCdigitRow()
530 {
531   Reset();
532 }
533 ////////////////////////////////////////////////////////////////////////
534 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
535 {
536   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
537   return (*this);
538 }
539 ////////////////////////////////////////////////////////////////////////
540 void AliTPCdigitRow::SetRow(Int_t row) 
541 {
542   if (row >= 8*kgRowBytes) {
543     cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
544     return;
545   }
546   Int_t iC = row/8;
547   Int_t iB = row%8;
548   SETBIT(fDig[iC],iB);
549 }
550
551 ////////////////////////////////////////////////////////////////////////
552 Bool_t AliTPCdigitRow::TestRow(Int_t row)
553 {
554 //
555 // return kTRUE if row is on
556 //
557   Int_t iC = row/8;
558   Int_t iB = row%8;
559   return TESTBIT(fDig[iC],iB);
560 }
561 ////////////////////////////////////////////////////////////////////////
562 Int_t AliTPCdigitRow::RowsOn(Int_t upto)
563 {
564 //
565 // returns number of rows with a digit  
566 // count only rows less equal row number upto
567 //
568   Int_t total = 0;
569   for (Int_t i = 0; i<kgRowBytes; i++) {
570     for (Int_t j = 0; j < 8; j++) {
571       if (i*8+j > upto) return total;
572       if (TESTBIT(fDig[i],j))  total++;
573     }
574   }
575   return total;
576 }
577 ////////////////////////////////////////////////////////////////////////
578 void AliTPCdigitRow::Reset()
579 {
580 //
581 // resets all rows to zero
582 //
583   for (Int_t i = 0; i<kgRowBytes; i++) {
584     fDig[i] <<= 8;
585   }
586 }
587 ////////////////////////////////////////////////////////////////////////
588 Int_t AliTPCdigitRow::Last()
589 {
590 //
591 // returns the last row number with a digit
592 // returns -1 if now digits 
593 //
594   for (Int_t i = kgRowBytes-1; i>=0; i--) {
595     for (Int_t j = 7; j >= 0; j--) {
596       if TESTBIT(fDig[i],j) return i*8+j;
597     }
598   }
599   return -1;
600 }
601 ////////////////////////////////////////////////////////////////////////
602 Int_t AliTPCdigitRow::First()
603 {
604 //
605 // returns the first row number with a digit
606 // returns -1 if now digits 
607 //
608   for (Int_t i = 0; i<kgRowBytes; i++) {
609     for (Int_t j = 0; j < 8; j++) {
610       if (TESTBIT(fDig[i],j)) return i*8+j;
611     }
612   }
613   return -1;
614 }
615
616 ////////////////////////////////////////////////////////////////////////
617 //
618 // end of implementation of a class AliTPCdigitRow
619 //
620 ////////////////////////////////////////////////////////////////////////
621   
622 ////////////////////////////////////////////////////////////////////////
623 AliGenInfoMaker::AliGenInfoMaker()
624 {
625   Reset();
626 }
627
628 ////////////////////////////////////////////////////////////////////////
629 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
630                                    Int_t nEvents, Int_t firstEvent)
631 {
632   Reset();
633   fFirstEventNr = firstEvent;
634   fEventNr = firstEvent;
635   fNEvents = nEvents;
636   //   fFnRes = fnRes;
637   sprintf(fFnRes,"%s",fnRes);
638   //
639   fLoader = AliRunLoader::Open(fnGalice);
640   if (gAlice){
641     delete gAlice->GetRunLoader();
642     delete gAlice;
643     gAlice = 0x0;
644   }
645   if (fLoader->LoadgAlice()){
646     cerr<<"Error occured while l"<<endl;
647   }
648   Int_t nall = fLoader->GetNumberOfEvents();
649   if (nEvents==0) {
650     nEvents =nall;
651     fNEvents=nall;
652     fFirstEventNr=0;
653   }    
654
655   if (nall<=0){
656     cerr<<"no events available"<<endl;
657     fEventNr = 0;
658     return;
659   }
660   if (firstEvent+nEvents>nall) {
661     fEventNr = nall-firstEvent;
662     cerr<<"restricted number of events availaible"<<endl;
663   }
664   AliMagF * magf = gAlice->Field();
665   AliTracker::SetFieldMap(magf,0);
666 }
667
668
669 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
670 {
671   // 
672   if (i<fNParticles) {
673     if (fGenInfo[i]) return  fGenInfo[i];
674     fGenInfo[i] = new AliMCInfo;  
675     fNInfos++;
676     return fGenInfo[i];
677   }
678   else 
679     return 0;  
680 }
681
682 ////////////////////////////////////////////////////////////////////////
683 void AliGenInfoMaker::Reset()
684 {
685   fEventNr = 0;
686   fNEvents = 0;
687   fTreeGenTracks = 0;
688   fFileGenTracks = 0;
689   fGenInfo = 0;
690   fNInfos  = 0;
691   //
692   //
693   fDebug = 0;
694   fVPrim[0] = -1000.;
695   fVPrim[1] = -1000.;
696   fVPrim[2] = -1000.;
697   fParamTPC = 0;
698 }
699 ////////////////////////////////////////////////////////////////////////
700 AliGenInfoMaker::~AliGenInfoMaker()
701 {
702   
703   if (fLoader){
704     fLoader->UnloadgAlice();
705     gAlice = 0;
706     delete fLoader;
707   }
708 }
709
710 Int_t  AliGenInfoMaker::SetIO()
711 {
712   //
713   // 
714   CreateTreeGenTracks();
715   if (!fTreeGenTracks) return 1;
716   //  AliTracker::SetFieldFactor(); 
717  
718   fParamTPC = GetTPCParam();
719   //
720   return 0;
721 }
722
723 ////////////////////////////////////////////////////////////////////////
724 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
725 {
726   //
727   // 
728   // SET INPUT
729   fLoader->SetEventNumber(eventNr);
730   //
731   fLoader->LoadHeader();
732   fLoader->LoadKinematics();  
733   fStack = fLoader->Stack();
734   //
735   fLoader->LoadTrackRefs();
736   fLoader->LoadHits();
737   fTreeTR = fLoader->TreeTR();
738   //
739   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
740   tpcl->LoadDigits();
741   fTreeD = tpcl->TreeD();  
742   return 0;
743 }
744
745 Int_t AliGenInfoMaker::CloseIOEvent()
746 {
747   fLoader->UnloadHeader();
748   fLoader->UnloadKinematics();
749   fLoader->UnloadTrackRefs();
750   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
751   tpcl->UnloadDigits();
752   return 0;
753 }
754
755 Int_t AliGenInfoMaker::CloseIO()
756 {
757   fLoader->UnloadgAlice();
758   return 0;
759 }
760
761
762
763 ////////////////////////////////////////////////////////////////////////
764 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
765 {
766   fNEvents = nEvents;
767   fFirstEventNr = firstEventNr;
768   return Exec();
769 }
770
771 ////////////////////////////////////////////////////////////////////////
772 Int_t AliGenInfoMaker::Exec()  
773 {
774   TStopwatch timer;
775   timer.Start();
776   Int_t status =SetIO();
777   if (status>0) return status;
778   //
779
780   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
781        fEventNr++) {
782     SetIO(fEventNr);
783     fNParticles = fStack->GetNtrack();
784     //
785     fGenInfo = new AliMCInfo*[fNParticles];
786     for (UInt_t i = 0; i<fNParticles; i++) {
787       fGenInfo[i]=0; 
788     }
789     //
790     cout<<"Start to process event "<<fEventNr<<endl;
791     cout<<"\tfNParticles = "<<fNParticles<<endl;
792     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
793     if (TreeTRLoop()>0) return 1;
794     //
795     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
796     if (TreeDLoop()>0) return 1;
797     //
798     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
799     if (TreeKLoop()>0) return 1;
800     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
801     //
802     if (BuildKinkInfo()>0) return 1;
803     if (BuildV0Info()>0) return 1;
804     //if (BuildHitLines()>0) return 1;
805     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
806     //
807     for (UInt_t i = 0; i<fNParticles; i++) {
808       if (fGenInfo[i]) delete fGenInfo[i]; 
809     }
810     delete []fGenInfo;
811     CloseIOEvent();
812   }
813   //
814   CloseIO();
815   CloseOutputFile();
816
817   cerr<<"Exec finished"<<endl;
818
819   timer.Stop();
820   timer.Print();
821   return 0;
822 }
823
824 ////////////////////////////////////////////////////////////////////////
825 void AliGenInfoMaker::CreateTreeGenTracks() 
826 {
827   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
828   if (!fFileGenTracks) {
829     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
830     return;
831   }
832   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
833   AliMCInfo * info = new AliMCInfo;
834   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
835   delete info; 
836   //
837   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
838   fTreeKinks = new TTree("genKinksTree","genKinksTree");  
839   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
840   delete kinkinfo;
841   //
842   AliGenV0Info *v0info = new AliGenV0Info;
843   fTreeV0 = new TTree("genV0Tree","genV0Tree");  
844   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
845   delete v0info;
846   //
847   //
848   AliTrackPointArray * points0 = new AliTrackPointArray;  
849   AliTrackPointArray * points1 = new AliTrackPointArray;  
850   AliTrackPointArray * points2 = new AliTrackPointArray;  
851   fTreeHitLines = new TTree("HitLines","HitLines");
852   fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
853   fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
854   fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
855   Int_t label=0;
856   fTreeHitLines->Branch("Label",&label,"label/I");
857   //
858   fTreeGenTracks->AutoSave();
859   fTreeKinks->AutoSave();
860   fTreeV0->AutoSave();
861   fTreeHitLines->AutoSave();
862 }
863 ////////////////////////////////////////////////////////////////////////
864 void AliGenInfoMaker::CloseOutputFile() 
865 {
866   if (!fFileGenTracks) {
867     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
868     return;
869   }
870   fFileGenTracks->cd();
871   fTreeGenTracks->Write();  
872   delete fTreeGenTracks;
873   fTreeKinks->Write();
874   delete fTreeKinks;
875   fTreeV0->Write();
876   delete fTreeV0;
877
878   fFileGenTracks->Close();
879   delete fFileGenTracks;
880   return;
881 }
882
883 ////////////////////////////////////////////////////////////////////////
884 Int_t AliGenInfoMaker::TreeKLoop()
885 {
886 //
887 // open the file with treeK
888 // loop over all entries there and save information about some tracks
889 //
890
891   AliStack * stack = fStack;
892   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
893   
894   if (fDebug > 0) {
895     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
896         <<fEventNr<<endl;
897   }  
898   Int_t  ipdg = 0;
899   TParticlePDG *ppdg = 0;
900   // not all generators give primary vertex position. Take the vertex
901   // of the particle 0 as primary vertex.
902   TDatabasePDG  pdg; //get pdg table  
903   //thank you very much root for this
904   TBranch * br = fTreeGenTracks->GetBranch("MC");
905   TParticle *particle = stack->ParticleFromTreeK(0);
906   fVPrim[0] = particle->Vx();
907   fVPrim[1] = particle->Vy();
908   fVPrim[2] = particle->Vz();
909   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
910     // load only particles with TR
911     AliMCInfo * info = GetInfo(iParticle);
912     if (!info) continue;
913     //////////////////////////////////////////////////////////////////////
914     info->fLabel = iParticle;
915     //
916     info->fParticle = *(stack->Particle(iParticle));
917     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
918     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
919     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; 
920     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
921                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
922     //
923     //
924     ipdg = info->fParticle.GetPdgCode();
925     info->fPdg = ipdg;
926     ppdg = pdg.GetParticle(ipdg);          
927     info->fEventNr = fEventNr;
928     info->Update();
929     //////////////////////////////////////////////////////////////////////    
930     br->SetAddress(&info);    
931     fTreeGenTracks->Fill();    
932   }
933   fTreeGenTracks->AutoSave();
934   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
935   return 0;
936 }
937
938
939 ////////////////////////////////////////////////////////////////////////////////////
940 Int_t  AliGenInfoMaker::BuildKinkInfo()
941 {
942   //
943   // Build tree with Kink Information
944   //
945   AliStack * stack = fStack;
946   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
947   
948   if (fDebug > 0) {
949     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
950         <<fEventNr<<endl;
951   }  
952   //  Int_t  ipdg = 0;
953   //TParticlePDG *ppdg = 0;
954   // not all generators give primary vertex position. Take the vertex
955   // of the particle 0 as primary vertex.
956   TDatabasePDG  pdg; //get pdg table  
957   //thank you very much root for this
958   TBranch * br = fTreeKinks->GetBranch("MC");
959   //
960   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
961   //
962   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
963     // load only particles with TR
964     AliMCInfo * info = GetInfo(iParticle);
965     if (!info) continue;
966     if (info->fCharge==0) continue;  
967     if (info->fTRdecay.P()<0.13) continue;  //momenta cut 
968     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
969     TParticle & particle = info->fParticle;
970     Int_t first = particle.GetDaughter(0);
971     if (first ==0) continue;
972
973     Int_t last = particle.GetDaughter(1);
974     if (last ==0) last=first;
975     AliMCInfo * info2 =  0;
976     AliMCInfo * dinfo =  0;
977     
978     for (Int_t id2=first;id2<=last;id2++){
979       info2 = GetInfo(id2);
980       if (!info2) continue;
981       TParticle &p2 = info2->fParticle;
982       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
983       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
984       if (!(p2.GetPDG())) continue;
985       dinfo =info2;
986      
987       kinkinfo->fMCm = (*info);
988       kinkinfo->fMCd = (*dinfo);
989       if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
990       kinkinfo->Update();
991       br->SetAddress(&kinkinfo);    
992       fTreeKinks->Fill();
993     }
994     /*
995     if (dinfo){
996       kinkinfo->fMCm = (*info);
997       kinkinfo->fMCd = (*dinfo);
998       kinkinfo->Update();
999       br->SetAddress(&kinkinfo);    
1000       fTreeKinks->Fill();
1001     }
1002     */
1003   }
1004   fTreeGenTracks->AutoSave();
1005   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1006   return 0;
1007 }
1008
1009
1010 ////////////////////////////////////////////////////////////////////////////////////
1011 Int_t  AliGenInfoMaker::BuildV0Info()
1012 {
1013   //
1014   // Build tree with V0 Information
1015   //
1016   AliStack * stack = fStack;
1017   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1018   
1019   if (fDebug > 0) {
1020     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1021         <<fEventNr<<endl;
1022   }  
1023   //  Int_t  ipdg = 0;
1024   //TParticlePDG *ppdg = 0;
1025   // not all generators give primary vertex position. Take the vertex
1026   // of the particle 0 as primary vertex.
1027   TDatabasePDG  pdg; //get pdg table  
1028   //thank you very much root for this
1029   TBranch * br = fTreeV0->GetBranch("MC");
1030   //
1031   AliGenV0Info * v0info = new AliGenV0Info;
1032   //
1033   //
1034   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1035     // load only particles with TR
1036     AliMCInfo * info = GetInfo(iParticle);
1037     if (!info) continue;
1038     if (info->fCharge==0) continue;  
1039     //
1040     //
1041     TParticle & particle = info->fParticle;
1042     Int_t mother = particle.GetMother(0);
1043     if (mother <=0) continue;
1044     //
1045     TParticle * motherparticle = stack->Particle(mother);
1046     if (!motherparticle) continue;
1047     //
1048     Int_t last = motherparticle->GetDaughter(1);
1049     if (last==(int)iParticle) continue;
1050     AliMCInfo * info2 =  info;
1051     AliMCInfo * dinfo =  GetInfo(last);
1052     if (!dinfo) continue;
1053     if (!dinfo->fParticle.GetPDG()) continue;
1054     if (!info2->fParticle.GetPDG()) continue;
1055     if (dinfo){
1056       v0info->fMCm = (*info);
1057       v0info->fMCd = (*dinfo);
1058       v0info->fMotherP = (*motherparticle);
1059       v0info->Update(fVPrim);
1060       br->SetAddress(&v0info);    
1061       fTreeV0->Fill();
1062     }
1063   }
1064   fTreeV0->AutoSave();
1065   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1066   return 0;
1067 }
1068
1069
1070
1071
1072 ////////////////////////////////////////////////////////////////////////
1073 Int_t AliGenInfoMaker::BuildHitLines()
1074 {
1075
1076 //
1077 // open the file with treeK
1078 // loop over all entries there and save information about some tracks
1079 //
1080
1081   AliStack * stack = fStack;
1082   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1083   
1084   if (fDebug > 0) {
1085     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1086         <<fEventNr<<endl;
1087   }  
1088 //   Int_t  ipdg = 0;
1089 //   // TParticlePDG *ppdg = 0;
1090 //   // not all generators give primary vertex position. Take the vertex
1091 //   // of the particle 0 as primary vertex.
1092 //   TDatabasePDG  pdg; //get pdg table  
1093 //   //thank you very much root for this
1094 //   AliTrackPointArray *tpcp = new AliTrackPointArray;
1095 //   AliTrackPointArray *trdp = new AliTrackPointArray;
1096 //   AliTrackPointArray *itsp = new AliTrackPointArray;
1097 //   Int_t label =0;
1098 //   //
1099 //   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1100 //   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");  
1101 //   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1102 //   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1103 //   brlabel->SetAddress(&label);
1104 //   brtpc->SetAddress(&tpcp);
1105 //   brtrd->SetAddress(&trdp);
1106 //   brits->SetAddress(&itsp);
1107 //   //
1108 //   AliDetector *dtpc = gAlice->GetDetector("TPC");
1109 //   AliDetector *dtrd = gAlice->GetDetector("TRD");
1110 //   AliDetector *dits = gAlice->GetDetector("ITS");
1111  
1112 //   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1113 //     // load only particles with TR
1114 //     AliMCInfo * info = GetInfo(iParticle);
1115 //     if (!info) continue;
1116 //     Int_t primpart = info->fPrimPart;
1117 //     ipdg = info->fParticle.GetPdgCode();
1118 //     label = iParticle;
1119 //     //
1120 //     gAlice->ResetHits();
1121 //     tpcp->Reset();
1122 //     itsp->Reset();
1123 //     trdp->Reset();
1124 //     tpcp->fLabel1 = ipdg;
1125 //     trdp->fLabel1 = ipdg;
1126 //     itsp->fLabel1 = ipdg;
1127 //     if (dtpc->TreeH()->GetEvent(primpart)){
1128 //       dtpc->LoadPoints(primpart);
1129 //       tpcp->Reset(dtpc,iParticle);
1130 //     }
1131 //     if (dtrd->TreeH()->GetEvent(primpart)){
1132 //       dtrd->LoadPoints(primpart);
1133 //       trdp->Reset(dtrd,iParticle);
1134 //     }
1135 //     if (dits->TreeH()->GetEvent(primpart)){
1136 //       dits->LoadPoints(primpart);
1137 //       itsp->Reset(dits,iParticle);
1138 //     }    
1139 //     //    
1140 //     fTreeHitLines->Fill();
1141 //     dtpc->ResetPoints();
1142 //     dtrd->ResetPoints();
1143 //     dits->ResetPoints();
1144 //   }
1145 //   delete tpcp;
1146 //   delete trdp;
1147 //   delete itsp;
1148 //   fTreeHitLines->AutoSave();
1149 //   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1150   return 0;
1151 }
1152
1153
1154 ////////////////////////////////////////////////////////////////////////
1155 Int_t AliGenInfoMaker::TreeDLoop()
1156 {
1157   //
1158   // open the file with treeD
1159   // loop over all entries there and save information about some tracks
1160   //
1161   
1162   Int_t nInnerSector = fParamTPC->GetNInnerSector();
1163   Int_t rowShift = 0;
1164   Int_t zero=fParamTPC->GetZeroSup()+6;  
1165   //
1166   //
1167   AliSimDigits digitsAddress, *digits=&digitsAddress;
1168   fTreeD->GetBranch("Segment")->SetAddress(&digits);
1169   
1170   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1171   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1172   for (Int_t i=0; i<sectorsByRows; i++) {
1173     if (!fTreeD->GetEvent(i)) continue;
1174     Int_t sec,row;
1175     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1176     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
1177     // here I expect that upper sectors follow lower sectors
1178     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1179     //
1180     digits->ExpandTrackBuffer();
1181     digits->First();        
1182     do {
1183       Int_t iRow=digits->CurrentRow();
1184       Int_t iColumn=digits->CurrentColumn();
1185       Short_t digitValue = digits->CurrentDigit();
1186       if (digitValue >= zero) {
1187         Int_t label;
1188         for (Int_t j = 0; j<3; j++) {
1189           //      label = digits->GetTrackID(iRow,iColumn,j);
1190           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;       
1191           if (label >= (Int_t)fNParticles) { //don't label from bakground event
1192             continue;
1193           }
1194           if (label >= 0 && label <= (Int_t)fNParticles) {
1195             if (fDebug > 6 ) {
1196               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1197                   <<sec<<" "
1198                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1199                   <<" "<<row<<endl;
1200             }   
1201             AliMCInfo * info = GetInfo(label);
1202             if (info){
1203               info->fTPCRow.SetRow(row+rowShift);
1204             }
1205           }
1206         }
1207       }
1208     } while (digits->Next());
1209   }
1210   
1211   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;  
1212   return 0;
1213 }
1214
1215
1216 ////////////////////////////////////////////////////////////////////////
1217 Int_t AliGenInfoMaker::TreeTRLoop()
1218 {
1219   //
1220   // loop over TrackReferences and store the first one for each track
1221   //  
1222   TTree * treeTR = fTreeTR;
1223   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1224   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1225   //
1226   //
1227   //track references for TPC
1228   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
1229   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
1230   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
1231   TClonesArray* TOFArrayTR = new TClonesArray("AliTrackReference");
1232   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
1233   //
1234   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
1235   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
1236   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
1237   if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
1238   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
1239   //
1240   //
1241   //
1242   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1243     treeTR->GetEntry(iPrimPart);
1244     //
1245     // Loop over TPC references
1246     //
1247     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
1248       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);            
1249       //
1250       if (trackRef->TestBit(BIT(2))){  
1251         //if decay 
1252         if (trackRef->P()<fgTPCPtCut) continue;
1253         Int_t label = trackRef->GetTrack(); 
1254         AliMCInfo * info = GetInfo(label);
1255         if (!info) info = MakeInfo(label);
1256         info->fTRdecay = *trackRef;      
1257       }
1258       //
1259       if (trackRef->P()<fgTPCPtCut) continue;
1260       Int_t label = trackRef->GetTrack();      
1261       AliMCInfo * info = GetInfo(label);
1262       if (!info) info = MakeInfo(label);
1263       if (!info) continue;
1264       info->fPrimPart =  iPrimPart;
1265       TClonesArray & arr = *(info->fTPCReferences);
1266       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1267     }
1268     //
1269     // Loop over ITS references
1270     //
1271     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
1272       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);            
1273       // 
1274       //
1275       if (trackRef->P()<fgTPCPtCut) continue;
1276       Int_t label = trackRef->GetTrack();      
1277       AliMCInfo * info = GetInfo(label);
1278       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
1279       if (!info) info = MakeInfo(label);
1280       if (!info) continue;
1281       info->fPrimPart =  iPrimPart;
1282       TClonesArray & arr = *(info->fITSReferences);
1283       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1284     }
1285     //
1286     // Loop over TRD references
1287     //
1288     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
1289       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);            
1290       //
1291       if (trackRef->P()<fgTPCPtCut) continue;
1292       Int_t label = trackRef->GetTrack();      
1293       AliMCInfo * info = GetInfo(label);
1294       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
1295       if (!info) info = MakeInfo(label);
1296       if (!info) continue;
1297       info->fPrimPart =  iPrimPart;
1298       TClonesArray & arr = *(info->fTRDReferences);
1299       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1300     }
1301     //
1302     // Loop over TOF references
1303     //
1304     for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
1305       AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);            
1306       Int_t label = trackRef->GetTrack();      
1307       AliMCInfo * info = GetInfo(label);
1308       if (!info){
1309         if (trackRef->Pt()<fgTPCPtCut) continue;
1310         if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
1311       }
1312       if (!info) info = MakeInfo(label);
1313       if (!info) continue;
1314       info->fPrimPart =  iPrimPart;
1315       TClonesArray & arr = *(info->fTOFReferences);
1316       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
1317     }
1318     //
1319     // get dacay position
1320     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
1321       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);      
1322       //
1323       Int_t label = trackRef->GetTrack();
1324       AliMCInfo * info = GetInfo(label);
1325       if (!info) continue;
1326       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
1327       //      if (TMath::Abs(trackRef.X());
1328       info->fTRdecay = *trackRef;      
1329     }
1330   }
1331   //
1332   TPCArrayTR->Delete();
1333   delete  TPCArrayTR;
1334   TRDArrayTR->Delete();
1335   delete  TRDArrayTR;
1336   TOFArrayTR->Delete();
1337   delete  TOFArrayTR;
1338
1339   ITSArrayTR->Delete();
1340   delete  ITSArrayTR;
1341   RunArrayTR->Delete();
1342   delete  RunArrayTR;
1343   //
1344   return 0;
1345 }
1346
1347 ////////////////////////////////////////////////////////////////////////
1348 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
1349                                     AliTPCParam *paramTPC) {
1350
1351   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1352   Int_t index[4];
1353   paramTPC->Transform0to1(x,index);
1354   paramTPC->Transform1to2(x,index);
1355   return x[0];
1356 }
1357 ////////////////////////////////////////////////////////////////////////
1358
1359
1360