First attempt to use AliAnalisysTask (Marian)
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfoMaker.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
26 ---Usage outside of the analysis framework
27
28 gSystem->Load("libANALYSIS.so")
29 gSystem->Load("libPWG1.so")
30 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
31 t->Exec();
32
33 */
34
35 #if !defined(__CINT__) || defined(__MAKECINT__)
36 #include <stdio.h>
37 #include <string.h>
38 //ROOT includes
39 #include "TROOT.h"
40 #include "Rtypes.h"
41 #include "TFile.h"
42 #include "TTree.h"
43 #include "TStopwatch.h"
44 #include "TParticle.h"
45
46 //ALIROOT includes
47 #include "AliMCEvent.h"
48 #include "AliMCEventHandler.h"
49
50 #include "AliRun.h"
51 #include "AliStack.h"
52 #include "AliSimDigits.h"
53 #include "AliTPCParam.h"
54 #include "AliTPC.h"
55 #include "AliTPCLoader.h"
56 #include "AliTrackReference.h"
57 #include "AliTPCParamSR.h"
58 #include "AliTrackPointArray.h"
59
60 #endif
61 #include "AliMCInfo.h" 
62 #include "AliGenV0Info.h" 
63 #include "AliGenKinkInfo.h" 
64 #include "AliGenInfoMaker.h" 
65 //
66 // 
67
68 ClassImp(AliGenInfoMaker)
69
70
71
72
73
74
75   
76 ////////////////////////////////////////////////////////////////////////
77 AliGenInfoMaker::AliGenInfoMaker(): 
78   fGenTracksArray(0),          //clones array with filtered particles
79   fGenKinkArray(0),            //clones array with filtered Kinks
80   fGenV0Array(0),              //clones array with filtered V0s
81   fDebug(0),                   //! debug flag  
82   fEventNr(0),                 //! current event number
83   fLabel(0),                   //! track label
84   fNEvents(0),                 //! number of events to process
85   fFirstEventNr(0),            //! first event to process
86   fNParticles(0),              //! number of particles in TreeK
87   fTreeGenTracks(0),           //! output tree with generated tracks
88   fTreeKinks(0),               //!  output tree with Kinks
89   fTreeV0(0),                  //!  output tree with V0
90   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
91   fLoader(0),                  //! pointer to the run loader
92   fTreeD(0),                   //! current tree with digits
93   fTreeTR(0),                  //! current tree with TR
94   fStack(0),                   //! current stack
95   fGenInfo(0),                 //! array with pointers to gen info
96   fNInfos(0),                  //! number of tracks with infos
97   fParamTPC(0),                //! AliTPCParam
98   fTPCPtCut(0.1),              //  TPC pt cut            
99   fITSPtCut(0.1),              //  ITS pt cut
100   fTRDPtCut(0.1),              //  TRD pt cut
101   fTOFPtCut(0.1)               //  TOF pt cut
102 {    
103   sprintf(fFnRes,"%s","genTracks.root");
104 }
105
106
107
108
109 ////////////////////////////////////////////////////////////////////////
110 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
111                                  Int_t nEvents, Int_t firstEvent):
112   fGenTracksArray(0),          //clones array with filtered particles
113   fGenKinkArray(0),            //clones array with filtered Kinks
114   fGenV0Array(0),              //clones array with filtered V0s
115   fDebug(0),                   //! debug flag  
116   fEventNr(0),                 //! current event number
117   fLabel(0),                   //! track label
118   fNEvents(0),                 //! number of events to process
119   fFirstEventNr(0),            //! first event to process
120   fNParticles(0),              //! number of particles in TreeK
121   fTreeGenTracks(0),           //! output tree with generated tracks
122   fTreeKinks(0),               //!  output tree with Kinks
123   fTreeV0(0),                  //!  output tree with V0
124   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
125   fLoader(0),                  //! pointer to the run loader
126   fTreeD(0),                   //! current tree with digits
127   fTreeTR(0),                  //! current tree with TR
128   fStack(0),                   //! current stack
129   fGenInfo(0),                 //! array with pointers to gen info
130   fNInfos(0),                  //! number of tracks with infos
131   fParamTPC(0),                //! AliTPCParam 
132   fTPCPtCut(0.1),  
133   fITSPtCut(0.1),  
134   fTRDPtCut(0.1), 
135   fTOFPtCut(0.1)
136 {
137   //
138   // 
139   //
140   fFirstEventNr = firstEvent;
141   fEventNr = firstEvent;
142   fNEvents = nEvents;
143   sprintf(fFnRes,"%s",fnRes);
144   //
145   //
146   //
147   fLoader = AliRunLoader::Open(fnGalice);
148   if (gAlice){
149     delete gAlice->GetRunLoader();
150     delete gAlice;
151     gAlice = 0x0;
152   }
153   if (fLoader->LoadgAlice()){
154     cerr<<"Error occured while l"<<endl;
155   }
156   Int_t nall = fLoader->GetNumberOfEvents();
157   if (nEvents==0) {
158     nEvents =nall;
159     fNEvents=nall;
160     fFirstEventNr=0;
161   }    
162   
163   if (nall<=0){
164     cerr<<"no events available"<<endl;
165     fEventNr = 0;
166     return;
167   }
168   if (firstEvent+nEvents>nall) {
169     fEventNr = nall-firstEvent;
170     cerr<<"restricted number of events availaible"<<endl;
171   }
172 }
173
174
175
176
177 ////////////////////////////////////////////////////////////////////////
178 AliGenInfoMaker::~AliGenInfoMaker()
179 {
180   //
181   // Destructor
182   //
183   
184   if (fLoader){
185     fLoader->UnloadgAlice();
186     gAlice = 0;
187     delete fLoader;
188   }
189 }
190
191
192 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
193 {
194   // 
195   // Make info structure for given particle index
196   //
197   if (i<fNParticles) {
198     if (fGenInfo[i]) return  fGenInfo[i];
199     fGenInfo[i] = new AliMCInfo;  
200     fNInfos++;
201     return fGenInfo[i];
202   }
203   else 
204     return 0;  
205 }
206
207
208 Int_t  AliGenInfoMaker::ProcessEvent(AliMCEventHandler* mcinfo){
209   //
210   // Process MC info from the task
211   //
212   if (!fParamTPC) {
213     SetIO();
214     fParamTPC = GetTPCParam();
215   }
216   fStack  = mcinfo->MCEvent()->Stack();
217   fTreeTR = mcinfo->TreeTR();
218   fTreeD  = 0;
219   // array with preprocessedprocessed information
220   fGenTracksArray = new TObjArray;
221   fGenKinkArray   = new TObjArray;
222   fGenV0Array     = new TObjArray;
223   //
224   ProcessEvent();
225   fEventNr++;
226   return 0;
227 }
228
229
230
231 Int_t AliGenInfoMaker::ProcessEvent(){
232   //
233   // Process Event 
234   //
235   fNParticles = fStack->GetNtrack();
236   //
237   fGenInfo = new AliMCInfo*[fNParticles];
238   for (UInt_t i = 0; i<fNParticles; i++) {
239     fGenInfo[i]=0; 
240   }
241   //
242   cout<<"Start to process event "<<fEventNr<<endl;
243   cout<<"\tfNParticles = "<<fNParticles<<endl;
244   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
245   if (TreeTRLoop()>0) return 1;
246   //
247   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
248   if (TreeDLoop()>0) return 1;
249   //
250   if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
251   if (TreeKLoop()>0) return 1;
252   if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
253   //
254   //  if (BuildKinkInfo()>0) return 1;
255   if (BuildV0Info()>0) return 1;
256   if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
257   //
258   for (UInt_t i = 0; i<fNParticles; i++) {
259     if (fGenInfo[i]) delete fGenInfo[i]; 
260   }
261   delete []fGenInfo;
262   return 0;
263 }
264
265
266
267
268
269
270
271 Int_t  AliGenInfoMaker::SetIO()
272 {
273   //
274   // Set IO for given event 
275   //
276   CreateTreeGenTracks();
277   if (!fTreeGenTracks) return 1;
278  
279   fParamTPC = GetTPCParam();
280   //
281   return 0;
282 }
283
284 ////////////////////////////////////////////////////////////////////////
285 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
286 {
287   //
288   // 
289   // SET INPUT
290   fLoader->SetEventNumber(eventNr);
291   //
292   fLoader->LoadHeader();
293   fLoader->LoadKinematics();  
294   fStack = fLoader->Stack();
295   //
296   fLoader->LoadTrackRefs();
297   fLoader->LoadHits();
298   fTreeTR = fLoader->TreeTR();
299   //
300   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
301   tpcl->LoadDigits();
302   fTreeD = tpcl->TreeD();  
303   return 0;
304 }
305
306 Int_t AliGenInfoMaker::CloseIOEvent()
307 {
308   //
309   // Close IO for current event
310   //
311   fLoader->UnloadHeader();
312   fLoader->UnloadKinematics();
313   fLoader->UnloadTrackRefs();
314   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
315   tpcl->UnloadDigits();
316   return 0;
317 }
318
319 Int_t AliGenInfoMaker::CloseIO()
320 {
321   fLoader->UnloadgAlice();
322   return 0;
323 }
324
325
326
327 ////////////////////////////////////////////////////////////////////////
328 Int_t AliGenInfoMaker::Exec()  
329 {
330   //
331   // Make a comparision MC tree
332   // Connect MC information -TPArticle - AliTrackRefernces ...
333   //
334   TStopwatch timer;
335   timer.Start();
336   Int_t status =SetIO();
337   if (status>0) return status;
338   //
339   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
340        fEventNr++) {
341     SetIO(fEventNr);
342     ProcessEvent();
343     CloseIOEvent();
344   }
345   //
346   CloseIO();
347   CloseOutputFile();
348   cerr<<"Exec finished"<<endl;
349   timer.Stop();
350   timer.Print();
351   return 0;
352 }
353
354
355
356
357
358
359 ////////////////////////////////////////////////////////////////////////
360 void AliGenInfoMaker::CreateTreeGenTracks() 
361 {
362   //
363   //
364   //
365   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
366   if (!fFileGenTracks) {
367     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
368     return;
369   }
370   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");  
371   AliMCInfo * info = new AliMCInfo;
372   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
373   delete info; 
374   //
375   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
376   fTreeKinks = new TTree("genKinksTree","genKinksTree");  
377   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
378   delete kinkinfo;
379   //
380   AliGenV0Info *v0info = new AliGenV0Info;
381   fTreeV0 = new TTree("genV0Tree","genV0Tree");  
382   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
383   delete v0info;
384   //
385   //
386   //
387   fTreeGenTracks->AutoSave();
388   fTreeKinks->AutoSave();
389   fTreeV0->AutoSave();
390 }
391
392
393
394 ////////////////////////////////////////////////////////////////////////
395 void AliGenInfoMaker::CloseOutputFile() 
396 {
397   //
398   // Close Output files
399   //
400   if (!fFileGenTracks) {
401     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
402     return;
403   }
404   fFileGenTracks->cd();
405   fTreeGenTracks->Write();  
406   delete fTreeGenTracks;
407   fTreeKinks->Write();
408   delete fTreeKinks;
409   fTreeV0->Write();
410   delete fTreeV0;
411
412   fFileGenTracks->Close();
413   delete fFileGenTracks;
414   return;
415 }
416
417 ////////////////////////////////////////////////////////////////////////
418 Int_t AliGenInfoMaker::TreeKLoop()
419 {
420 //
421 // open the file with treeK
422 // loop over all entries there and save information about some tracks
423 //
424
425   AliStack * stack = fStack;
426   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
427   
428   if (fDebug > 0) {
429     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
430         <<fEventNr<<endl;
431   }  
432   Int_t  ipdg = 0;
433   TParticlePDG *ppdg = 0;
434   // not all generators give primary vertex position. Take the vertex
435   // of the particle 0 as primary vertex.
436   TDatabasePDG  pdg; //get pdg table  
437
438   TParticle *particle = stack->ParticleFromTreeK(0);
439   fVPrim[0] = particle->Vx();
440   fVPrim[1] = particle->Vy();
441   fVPrim[2] = particle->Vz();
442   Int_t accepted =0;
443   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
444     // load only particles with TR
445     AliMCInfo * info = GetInfo(iParticle);
446     if (!info) continue;
447     //////////////////////////////////////////////////////////////////////
448     info->fLabel = iParticle;
449     //
450     info->fParticle = *(stack->Particle(iParticle));
451     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
452     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
453     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2]; 
454     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
455                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
456     //
457     //
458     ipdg = info->fParticle.GetPdgCode();
459     info->fPdg = ipdg;
460     ppdg = pdg.GetParticle(ipdg);          
461     info->fEventNr = fEventNr;
462     info->Update();
463     if (fGenTracksArray){
464       fGenTracksArray->AddLast(info->Clone());
465     }
466     accepted++;
467   }
468   //
469   // write the results to the tree - if specified
470   //
471   TBranch * br = fTreeGenTracks->GetBranch("MC");
472   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
473     // load only particles with TR
474     AliMCInfo * info = GetInfo(iParticle);
475     if (!info) continue;
476     //////////////////////////////////////////////////////////////////////    
477     br->SetAddress(&info);    
478     fTreeGenTracks->Fill();    
479   }
480   fTreeGenTracks->AutoSave();
481   //
482   //
483   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
484   return 0;
485 }
486
487
488 ////////////////////////////////////////////////////////////////////////////////////
489 Int_t  AliGenInfoMaker::BuildKinkInfo()
490 {
491   //
492   // Build tree with Kink Information
493   //
494   AliStack * stack = fStack;
495   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
496   
497   if (fDebug > 0) {
498     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
499         <<fEventNr<<endl;
500   }  
501   //  Int_t  ipdg = 0;
502   //TParticlePDG *ppdg = 0;
503   // not all generators give primary vertex position. Take the vertex
504   // of the particle 0 as primary vertex.
505   TDatabasePDG  pdg; //get pdg table  
506   //thank you very much root for this
507   TBranch * br = fTreeKinks->GetBranch("MC");
508   //
509   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
510   //
511   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
512     // load only particles with TR
513     AliMCInfo * info = GetInfo(iParticle);
514     if (!info) continue;
515     if (info->fCharge==0) continue;  
516     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
517     TParticle & particle = info->fParticle;
518     Int_t first = particle.GetDaughter(0);
519     if (first ==0) continue;
520
521     Int_t last = particle.GetDaughter(1);
522     if (last ==0) last=first;
523     AliMCInfo * info2 =  0;
524     AliMCInfo * dinfo =  0;
525     
526     for (Int_t id2=first;id2<=last;id2++){
527       info2 = GetInfo(id2);
528       if (!info2) continue;
529       TParticle &p2 = info2->fParticle;
530       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
531       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
532       if (!(p2.GetPDG())) continue;
533       dinfo =info2;
534      
535       kinkinfo->SetInfoMother(*info);
536       kinkinfo->SetInfoDaughter(*dinfo);
537       if (kinkinfo->GetMinus().GetParticle().GetPDG()==0 ||  kinkinfo->GetPlus().GetParticle().GetPDG()==0) continue;
538       kinkinfo->Update();
539       br->SetAddress(&kinkinfo);    
540       fTreeKinks->Fill();
541     }
542   }
543   fTreeKinks->AutoSave();
544   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
545   return 0;
546 }
547
548
549 ////////////////////////////////////////////////////////////////////////////////////
550 Int_t  AliGenInfoMaker::BuildV0Info()
551 {
552   //
553   // Build tree with V0 Information
554   //
555   AliStack * stack = fStack;
556   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
557   
558   if (fDebug > 0) {
559     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
560         <<fEventNr<<endl;
561   }  
562   //
563   TDatabasePDG  pdg; //get pdg table  
564   //thank you very much root for this
565   TBranch * br = fTreeV0->GetBranch("MC");
566   //
567   AliGenV0Info * v0info = new AliGenV0Info;
568   //
569   //
570   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {    
571     TParticle * mParticle = fStack->Particle(iParticle);
572     if (!mParticle) continue;
573     if (mParticle->GetPDG()==0) continue;
574     if (mParticle->GetPDG()->Charge()!=0) continue;  //only neutral particle
575     //
576     Int_t first = mParticle->GetDaughter(0);
577     Int_t last  = mParticle->GetDaughter(1);
578     if (first-last==0) continue;
579     //
580     AliMCInfo * info0 = GetInfo(first);
581     AliMCInfo * info1 = GetInfo(first+1);
582     if (info0 && info1){
583       if (info0->GetPdg()==0) continue;
584       if (info1->GetPdg()==0) continue;
585       v0info->SetInfoP(*info0);
586       v0info->SetInfoM(*info1);
587       v0info->SetMother(*mParticle);
588       v0info->Update(fVPrim);
589       if (fGenV0Array) fGenV0Array->AddLast(v0info);
590       br->SetAddress(&v0info);    
591       fTreeV0->Fill();
592     }
593   }
594   fTreeV0->AutoSave();
595   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
596   return 0;
597 }
598
599
600
601
602
603 ////////////////////////////////////////////////////////////////////////
604 Int_t AliGenInfoMaker::TreeDLoop()
605 {
606   //
607   // open the file with treeD
608   // loop over all entries there and save information about some tracks
609   //
610   if (!fTreeD) return 0;
611   Int_t nInnerSector = fParamTPC->GetNInnerSector();
612   Int_t rowShift = 0;
613   Int_t zero=fParamTPC->GetZeroSup()+6;  
614   //
615   //
616   AliSimDigits digitsAddress, *digits=&digitsAddress;
617   fTreeD->GetBranch("Segment")->SetAddress(&digits);
618   
619   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
620   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
621   for (Int_t i=0; i<sectorsByRows; i++) {
622     if (!fTreeD->GetEvent(i)) continue;
623     Int_t sec,row;
624     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
625     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
626     // here I expect that upper sectors follow lower sectors
627     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
628     //
629     digits->ExpandTrackBuffer();
630     digits->First();        
631     do {
632       Int_t iRow=digits->CurrentRow();
633       Int_t iColumn=digits->CurrentColumn();
634       Short_t digitValue = digits->CurrentDigit();
635       if (digitValue >= zero) {
636         Int_t label;
637         for (Int_t j = 0; j<3; j++) {
638           //      label = digits->GetTrackID(iRow,iColumn,j);
639           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;       
640           if (label >= (Int_t)fNParticles) { //don't label from bakground event
641             continue;
642           }
643           if (label >= 0 && label <= (Int_t)fNParticles) {
644             if (fDebug > 6 ) {
645               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
646                   <<sec<<" "
647                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
648                   <<" "<<row<<endl;
649             }   
650             AliMCInfo * info = GetInfo(label);
651             if (info){
652               info->fTPCRow.SetRow(row+rowShift);
653             }
654           }
655         }
656       }
657     } while (digits->Next());
658   }
659   
660   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;  
661   return 0;
662 }
663
664
665
666
667 ////////////////////////////////////////////////////////////////////////
668 Int_t AliGenInfoMaker::TreeTRLoop()
669 {
670   //
671   // loop over TrackReferences and store the first one for each track
672   //  
673   TTree * treeTR = fTreeTR;
674   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
675   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
676   //
677   //
678   //
679   TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
680   //
681   if (treeTR->GetBranch("TrackReferences")) treeTR->GetBranch("TrackReferences")->SetAddress(&runArrayTR);
682   //
683   //
684   //
685   for (Int_t ipart = 0; ipart<fStack->GetNtrack(); ipart++) {
686     TParticle * part = fStack->Particle(ipart);
687     if (!part) continue;
688     if (part->GetPDG()==0) continue;
689     if (part->GetPDG()->Charge()==0) continue;
690     
691     if (part->Pt()<fITSPtCut) continue;
692     if (part->R()>250.) continue;
693     if (TMath::Abs(part->Vz())>250.) continue;
694     if (TMath::Abs(part->Pz()/part->Pt())>2.5) continue;     
695     AliMCInfo * info = MakeInfo(ipart); 
696   }
697   //
698   //
699
700   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
701     treeTR->GetEntry(iPrimPart);
702     //
703     // gettrack references
704     //
705     for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
706       AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);      
707       //
708       Int_t label = trackRef->GetTrack();
709       AliMCInfo * cinfo = GetInfo(label);
710       if (cinfo) cinfo->CalcTPCrows(runArrayTR);
711       //
712       // TPC
713       //
714       if (trackRef->DetectorId()== AliTrackReference::kTPC){
715         //
716         Int_t label = trackRef->GetTrack();      
717         AliMCInfo * info = GetInfo(label);
718         if (!info && trackRef->Pt()<fTPCPtCut) continue;
719         if (!info) {
720           info = MakeInfo(label);
721         }
722         if (!info) continue;
723         info->fPrimPart =  iPrimPart;
724         TClonesArray & arr = *(info->fTPCReferences);
725         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
726       }
727       //
728       // ITS
729       //
730       if (trackRef->DetectorId()== AliTrackReference::kITS){
731         //
732         Int_t label = trackRef->GetTrack();      
733         AliMCInfo * info = GetInfo(label);
734         if (!info && trackRef->Pt()<fITSPtCut) continue;
735         if (!info) info = MakeInfo(label);
736         if (!info) continue;
737         info->fPrimPart =  iPrimPart;
738         TClonesArray & arr = *(info->fITSReferences);
739         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
740       }
741       //
742       // TRD
743       //
744       if (trackRef->DetectorId()== AliTrackReference::kTRD){
745         //
746         Int_t label = trackRef->GetTrack();      
747         AliMCInfo * info = GetInfo(label);
748         if (!info&&trackRef->P()<fTRDPtCut) continue;
749         if (!info) info = MakeInfo(label);
750         if (!info) continue;
751         info->fPrimPart =  iPrimPart;
752         TClonesArray & arr = *(info->fTRDReferences);
753         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
754       }
755       //
756       // TOF
757       //
758       if (trackRef->DetectorId()== AliTrackReference::kTOF){
759         //
760         Int_t label = trackRef->GetTrack();      
761         AliMCInfo * info = GetInfo(label);
762         if (!info&&trackRef->P()<fTOFPtCut) continue;
763         if (!info) info = MakeInfo(label);
764         if (!info) continue;
765         info->fPrimPart =  iPrimPart;
766         TClonesArray & arr = *(info->fTOFReferences);
767         new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);     
768       }
769       //
770       // decay case
771       //
772       if (trackRef->DetectorId()== AliTrackReference::kDisappeared){
773         //
774         AliMCInfo * info = GetInfo(label);
775         if (!info) continue;
776         info->fTRdecay = *trackRef;             
777       }
778     }
779   }
780   //
781   //
782   return 0;
783 }
784
785 ////////////////////////////////////////////////////////////////////////
786 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
787                                     AliTPCParam *paramTPC) const {
788
789   //
790   // Transformation from Gloabal to local tacking system
791   //
792   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
793   Int_t index[4];
794   paramTPC->Transform0to1(x,index);
795   paramTPC->Transform1to2(x,index);
796   return x[0];
797 }
798 ////////////////////////////////////////////////////////////////////////
799
800
801
802 AliTPCParam * AliGenInfoMaker::GetTPCParam(){
803   //
804   // create default TPC parameters
805   //
806   AliTPCParamSR * par = new AliTPCParamSR;
807   par->Update();
808   return par;
809 }
810
811