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