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