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