First attempt to use AliAnalisysTask (Marian)
[u/mrichter/AliRoot.git] / PWG1 / AliGenInfoMaker.cxx
CommitLineData
d92975ba 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
20Origin: marian.ivanov@cern.ch
21
22Generate complex MC information - used for Comparison later on
23How to use it?
24
d390cc7e 25
26---Usage outside of the analysis framework
27
28gSystem->Load("libANALYSIS.so")
d92975ba 29gSystem->Load("libPWG1.so")
30AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0)
31t->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"
d92975ba 43#include "TStopwatch.h"
44#include "TParticle.h"
d92975ba 45
46//ALIROOT includes
d390cc7e 47#include "AliMCEvent.h"
48#include "AliMCEventHandler.h"
49
d92975ba 50#include "AliRun.h"
51#include "AliStack.h"
52#include "AliSimDigits.h"
53#include "AliTPCParam.h"
54#include "AliTPC.h"
55#include "AliTPCLoader.h"
d92975ba 56#include "AliTrackReference.h"
57#include "AliTPCParamSR.h"
d92975ba 58#include "AliTrackPointArray.h"
59
60#endif
b9e8c174 61#include "AliMCInfo.h"
62#include "AliGenV0Info.h"
63#include "AliGenKinkInfo.h"
d92975ba 64#include "AliGenInfoMaker.h"
65//
66//
67
68ClassImp(AliGenInfoMaker)
69
70
71
72
73
74
75
76////////////////////////////////////////////////////////////////////////
d390cc7e 77AliGenInfoMaker::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
d92975ba 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
d92975ba 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
d390cc7e 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");
d92975ba 104}
105
d390cc7e 106
107
108
d92975ba 109////////////////////////////////////////////////////////////////////////
110AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
111 Int_t nEvents, Int_t firstEvent):
d390cc7e 112 fGenTracksArray(0), //clones array with filtered particles
113 fGenKinkArray(0), //clones array with filtered Kinks
114 fGenV0Array(0), //clones array with filtered V0s
d92975ba 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
d92975ba 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
d390cc7e 132 fTPCPtCut(0.1),
d92975ba 133 fITSPtCut(0.1),
134 fTRDPtCut(0.1),
135 fTOFPtCut(0.1)
d92975ba 136{
137 //
138 //
139 //
140 fFirstEventNr = firstEvent;
141 fEventNr = firstEvent;
142 fNEvents = nEvents;
143 sprintf(fFnRes,"%s",fnRes);
144 //
d390cc7e 145 //
146 //
d92975ba 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 }
d390cc7e 162
d92975ba 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 }
d390cc7e 172}
173
174
175
176
177////////////////////////////////////////////////////////////////////////
178AliGenInfoMaker::~AliGenInfoMaker()
179{
180 //
181 // Destructor
182 //
183
184 if (fLoader){
185 fLoader->UnloadgAlice();
186 gAlice = 0;
187 delete fLoader;
188 }
d92975ba 189}
190
191
192AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
193{
194 //
9f3282ed 195 // Make info structure for given particle index
196 //
d92975ba 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
d390cc7e 207
208Int_t AliGenInfoMaker::ProcessEvent(AliMCEventHandler* mcinfo){
9f3282ed 209 //
d390cc7e 210 // Process MC info from the task
9f3282ed 211 //
d390cc7e 212 if (!fParamTPC) {
213 SetIO();
214 fParamTPC = GetTPCParam();
d92975ba 215 }
d390cc7e 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;
d92975ba 227}
228
d390cc7e 229
230
231Int_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
d92975ba 271Int_t AliGenInfoMaker::SetIO()
272{
273 //
9f3282ed 274 // Set IO for given event
275 //
d92975ba 276 CreateTreeGenTracks();
277 if (!fTreeGenTracks) return 1;
d92975ba 278
279 fParamTPC = GetTPCParam();
280 //
281 return 0;
282}
283
284////////////////////////////////////////////////////////////////////////
285Int_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
306Int_t AliGenInfoMaker::CloseIOEvent()
307{
9f3282ed 308 //
309 // Close IO for current event
310 //
d92975ba 311 fLoader->UnloadHeader();
312 fLoader->UnloadKinematics();
313 fLoader->UnloadTrackRefs();
314 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
315 tpcl->UnloadDigits();
316 return 0;
317}
318
319Int_t AliGenInfoMaker::CloseIO()
320{
321 fLoader->UnloadgAlice();
322 return 0;
323}
324
325
326
327////////////////////////////////////////////////////////////////////////
d92975ba 328Int_t AliGenInfoMaker::Exec()
329{
9f3282ed 330 //
331 // Make a comparision MC tree
332 // Connect MC information -TPArticle - AliTrackRefernces ...
333 //
d92975ba 334 TStopwatch timer;
335 timer.Start();
336 Int_t status =SetIO();
337 if (status>0) return status;
338 //
d92975ba 339 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
340 fEventNr++) {
341 SetIO(fEventNr);
d390cc7e 342 ProcessEvent();
d92975ba 343 CloseIOEvent();
344 }
345 //
346 CloseIO();
347 CloseOutputFile();
d92975ba 348 cerr<<"Exec finished"<<endl;
d92975ba 349 timer.Stop();
350 timer.Print();
351 return 0;
352}
353
d390cc7e 354
355
356
357
358
d92975ba 359////////////////////////////////////////////////////////////////////////
360void AliGenInfoMaker::CreateTreeGenTracks()
361{
9f3282ed 362 //
363 //
364 //
d92975ba 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 //
d92975ba 386 //
387 fTreeGenTracks->AutoSave();
388 fTreeKinks->AutoSave();
389 fTreeV0->AutoSave();
d92975ba 390}
d390cc7e 391
392
393
d92975ba 394////////////////////////////////////////////////////////////////////////
395void AliGenInfoMaker::CloseOutputFile()
396{
9f3282ed 397 //
398 // Close Output files
399 //
d92975ba 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////////////////////////////////////////////////////////////////////////
418Int_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
d390cc7e 437
d92975ba 438 TParticle *particle = stack->ParticleFromTreeK(0);
439 fVPrim[0] = particle->Vx();
440 fVPrim[1] = particle->Vy();
441 fVPrim[2] = particle->Vz();
d390cc7e 442 Int_t accepted =0;
d92975ba 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();
d390cc7e 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;
d92975ba 476 //////////////////////////////////////////////////////////////////////
477 br->SetAddress(&info);
478 fTreeGenTracks->Fill();
479 }
480 fTreeGenTracks->AutoSave();
d390cc7e 481 //
482 //
d92975ba 483 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
484 return 0;
485}
486
487
488////////////////////////////////////////////////////////////////////////////////////
489Int_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;
d92975ba 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 }
d92975ba 542 }
d390cc7e 543 fTreeKinks->AutoSave();
d92975ba 544 if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
545 return 0;
546}
547
548
549////////////////////////////////////////////////////////////////////////////////////
550Int_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 }
d390cc7e 562 //
d92975ba 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 //
d390cc7e 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
d92975ba 575 //
d390cc7e 576 Int_t first = mParticle->GetDaughter(0);
577 Int_t last = mParticle->GetDaughter(1);
578 if (first-last==0) continue;
d92975ba 579 //
d390cc7e 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);
d92975ba 588 v0info->Update(fVPrim);
d390cc7e 589 if (fGenV0Array) fGenV0Array->AddLast(v0info);
d92975ba 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
d92975ba 602
603////////////////////////////////////////////////////////////////////////
604Int_t AliGenInfoMaker::TreeDLoop()
605{
606 //
607 // open the file with treeD
608 // loop over all entries there and save information about some tracks
609 //
d390cc7e 610 if (!fTreeD) return 0;
d92975ba 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
999d8278 665
666
de6c9df4 667////////////////////////////////////////////////////////////////////////
d390cc7e 668Int_t AliGenInfoMaker::TreeTRLoop()
de6c9df4 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 //
999d8278 685 for (Int_t ipart = 0; ipart<fStack->GetNtrack(); ipart++) {
686 TParticle * part = fStack->Particle(ipart);
687 if (!part) continue;
d390cc7e 688 if (part->GetPDG()==0) continue;
689 if (part->GetPDG()->Charge()==0) continue;
690
999d8278 691 if (part->Pt()<fITSPtCut) continue;
692 if (part->R()>250.) continue;
693 if (TMath::Abs(part->Vz())>250.) continue;
d390cc7e 694 if (TMath::Abs(part->Pz()/part->Pt())>2.5) continue;
695 AliMCInfo * info = MakeInfo(ipart);
999d8278 696 }
697 //
698 //
699
de6c9df4 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();
d390cc7e 709 AliMCInfo * cinfo = GetInfo(label);
710 if (cinfo) cinfo->CalcTPCrows(runArrayTR);
de6c9df4 711 //
712 // TPC
713 //
714 if (trackRef->DetectorId()== AliTrackReference::kTPC){
715 //
de6c9df4 716 Int_t label = trackRef->GetTrack();
717 AliMCInfo * info = GetInfo(label);
999d8278 718 if (!info && trackRef->Pt()<fTPCPtCut) continue;
d390cc7e 719 if (!info) {
720 info = MakeInfo(label);
721 }
de6c9df4 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 //
de6c9df4 732 Int_t label = trackRef->GetTrack();
733 AliMCInfo * info = GetInfo(label);
999d8278 734 if (!info && trackRef->Pt()<fITSPtCut) continue;
de6c9df4 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 //
de6c9df4 746 Int_t label = trackRef->GetTrack();
747 AliMCInfo * info = GetInfo(label);
999d8278 748 if (!info&&trackRef->P()<fTRDPtCut) continue;
de6c9df4 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 //
de6c9df4 760 Int_t label = trackRef->GetTrack();
761 AliMCInfo * info = GetInfo(label);
999d8278 762 if (!info&&trackRef->P()<fTOFPtCut) continue;
de6c9df4 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 //
de6c9df4 774 AliMCInfo * info = GetInfo(label);
775 if (!info) continue;
776 info->fTRdecay = *trackRef;
777 }
de6c9df4 778 }
779 }
780 //
781 //
d92975ba 782 return 0;
783}
784
785////////////////////////////////////////////////////////////////////////
786Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
787 AliTPCParam *paramTPC) const {
788
9f3282ed 789 //
790 // Transformation from Gloabal to local tacking system
791 //
d92975ba 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
802AliTPCParam * AliGenInfoMaker::GetTPCParam(){
9f3282ed 803 //
804 // create default TPC parameters
805 //
d92975ba 806 AliTPCParamSR * par = new AliTPCParamSR;
807 par->Update();
808 return par;
809}
810
811