]>
Commit | Line | Data |
---|---|---|
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 | ||
20 | Origin: marian.ivanov@cern.ch | |
21 | ||
22 | Generate complex MC information - used for Comparison later on | |
23 | How to use it? | |
24 | ||
d390cc7e | 25 | |
26 | ---Usage outside of the analysis framework | |
27 | ||
28 | gSystem->Load("libANALYSIS.so") | |
d92975ba | 29 | gSystem->Load("libPWG1.so") |
30 | AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",0,0) | |
31 | t->Exec(); | |
32 | ||
33 | */ | |
34 | ||
35 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
36 | #include <stdio.h> | |
37 | #include <string.h> | |
38 | //ROOT includes | |
39 | #include "TROOT.h" | |
40 | #include "Rtypes.h" | |
41 | #include "TFile.h" | |
42 | #include "TTree.h" | |
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 | ||
68 | ClassImp(AliGenInfoMaker) | |
69 | ||
70 | ||
71 | ||
72 | ||
73 | ||
74 | ||
75 | ||
76 | //////////////////////////////////////////////////////////////////////// | |
d390cc7e | 77 | AliGenInfoMaker::AliGenInfoMaker(): |
78 | fGenTracksArray(0), //clones array with filtered particles | |
79 | fGenKinkArray(0), //clones array with filtered Kinks | |
80 | fGenV0Array(0), //clones array with filtered V0s | |
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 | //////////////////////////////////////////////////////////////////////// |
110 | AliGenInfoMaker::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 | //////////////////////////////////////////////////////////////////////// | |
178 | AliGenInfoMaker::~AliGenInfoMaker() | |
179 | { | |
180 | // | |
181 | // Destructor | |
182 | // | |
183 | ||
184 | if (fLoader){ | |
185 | fLoader->UnloadgAlice(); | |
186 | gAlice = 0; | |
187 | delete fLoader; | |
188 | } | |
d92975ba | 189 | } |
190 | ||
191 | ||
192 | AliMCInfo * 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 | |
208 | Int_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 | ||
231 | Int_t AliGenInfoMaker::ProcessEvent(){ | |
232 | // | |
233 | // Process Event | |
234 | // | |
235 | fNParticles = fStack->GetNtrack(); | |
236 | // | |
237 | fGenInfo = new AliMCInfo*[fNParticles]; | |
238 | for (UInt_t i = 0; i<fNParticles; i++) { | |
239 | fGenInfo[i]=0; | |
240 | } | |
241 | // | |
242 | cout<<"Start to process event "<<fEventNr<<endl; | |
243 | cout<<"\tfNParticles = "<<fNParticles<<endl; | |
244 | if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl; | |
245 | if (TreeTRLoop()>0) return 1; | |
246 | // | |
247 | if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl; | |
248 | if (TreeDLoop()>0) return 1; | |
249 | // | |
250 | if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl; | |
251 | if (TreeKLoop()>0) return 1; | |
252 | if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl; | |
253 | // | |
254 | // if (BuildKinkInfo()>0) return 1; | |
255 | if (BuildV0Info()>0) return 1; | |
256 | if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl; | |
257 | // | |
258 | for (UInt_t i = 0; i<fNParticles; i++) { | |
259 | if (fGenInfo[i]) delete fGenInfo[i]; | |
260 | } | |
261 | delete []fGenInfo; | |
262 | return 0; | |
263 | } | |
264 | ||
265 | ||
266 | ||
267 | ||
268 | ||
269 | ||
270 | ||
d92975ba | 271 | Int_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 | //////////////////////////////////////////////////////////////////////// | |
285 | Int_t AliGenInfoMaker::SetIO(Int_t eventNr) | |
286 | { | |
287 | // | |
288 | // | |
289 | // SET INPUT | |
290 | fLoader->SetEventNumber(eventNr); | |
291 | // | |
292 | fLoader->LoadHeader(); | |
293 | fLoader->LoadKinematics(); | |
294 | fStack = fLoader->Stack(); | |
295 | // | |
296 | fLoader->LoadTrackRefs(); | |
297 | fLoader->LoadHits(); | |
298 | fTreeTR = fLoader->TreeTR(); | |
299 | // | |
300 | AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader"); | |
301 | tpcl->LoadDigits(); | |
302 | fTreeD = tpcl->TreeD(); | |
303 | return 0; | |
304 | } | |
305 | ||
306 | Int_t AliGenInfoMaker::CloseIOEvent() | |
307 | { | |
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 | ||
319 | Int_t AliGenInfoMaker::CloseIO() | |
320 | { | |
321 | fLoader->UnloadgAlice(); | |
322 | return 0; | |
323 | } | |
324 | ||
325 | ||
326 | ||
d92975ba | 327 | //////////////////////////////////////////////////////////////////////// |
328 | Int_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 | //////////////////////////////////////////////////////////////////////// |
360 | void 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 | //////////////////////////////////////////////////////////////////////// |
395 | void 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 | //////////////////////////////////////////////////////////////////////// | |
418 | Int_t AliGenInfoMaker::TreeKLoop() | |
419 | { | |
420 | // | |
421 | // open the file with treeK | |
422 | // loop over all entries there and save information about some tracks | |
423 | // | |
424 | ||
425 | AliStack * stack = fStack; | |
426 | if (!stack) {cerr<<"Stack was not found!\n"; return 1;} | |
427 | ||
428 | if (fDebug > 0) { | |
429 | cout<<"There are "<<fNParticles<<" primary and secondary particles in event " | |
430 | <<fEventNr<<endl; | |
431 | } | |
432 | Int_t ipdg = 0; | |
433 | TParticlePDG *ppdg = 0; | |
434 | // not all generators give primary vertex position. Take the vertex | |
435 | // of the particle 0 as primary vertex. | |
436 | TDatabasePDG pdg; //get pdg table | |
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 | //////////////////////////////////////////////////////////////////////////////////// | |
489 | Int_t AliGenInfoMaker::BuildKinkInfo() | |
490 | { | |
491 | // | |
492 | // Build tree with Kink Information | |
493 | // | |
494 | AliStack * stack = fStack; | |
495 | if (!stack) {cerr<<"Stack was not found!\n"; return 1;} | |
496 | ||
497 | if (fDebug > 0) { | |
498 | cout<<"There are "<<fNParticles<<" primary and secondary particles in event " | |
499 | <<fEventNr<<endl; | |
500 | } | |
501 | // Int_t ipdg = 0; | |
502 | //TParticlePDG *ppdg = 0; | |
503 | // not all generators give primary vertex position. Take the vertex | |
504 | // of the particle 0 as primary vertex. | |
505 | TDatabasePDG pdg; //get pdg table | |
506 | //thank you very much root for this | |
507 | TBranch * br = fTreeKinks->GetBranch("MC"); | |
508 | // | |
509 | AliGenKinkInfo * kinkinfo = new AliGenKinkInfo; | |
510 | // | |
511 | for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) { | |
512 | // load only particles with TR | |
513 | AliMCInfo * info = GetInfo(iParticle); | |
514 | if (!info) continue; | |
515 | if (info->fCharge==0) continue; | |
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 | //////////////////////////////////////////////////////////////////////////////////// | |
550 | Int_t AliGenInfoMaker::BuildV0Info() | |
551 | { | |
552 | // | |
553 | // Build tree with V0 Information | |
554 | // | |
555 | AliStack * stack = fStack; | |
556 | if (!stack) {cerr<<"Stack was not found!\n"; return 1;} | |
557 | ||
558 | if (fDebug > 0) { | |
559 | cout<<"There are "<<fNParticles<<" primary and secondary particles in event " | |
560 | <<fEventNr<<endl; | |
561 | } | |
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 | //////////////////////////////////////////////////////////////////////// | |
604 | Int_t AliGenInfoMaker::TreeDLoop() | |
605 | { | |
606 | // | |
607 | // open the file with treeD | |
608 | // loop over all entries there and save information about some tracks | |
609 | // | |
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 | 668 | Int_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; |
cd875161 | 695 | //AliMCInfo * info = ; |
696 | MakeInfo(ipart); | |
999d8278 | 697 | } |
698 | // | |
699 | // | |
700 | ||
de6c9df4 | 701 | for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) { |
702 | treeTR->GetEntry(iPrimPart); | |
703 | // | |
704 | // gettrack references | |
705 | // | |
706 | for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) { | |
707 | AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef); | |
708 | // | |
709 | Int_t label = trackRef->GetTrack(); | |
d390cc7e | 710 | AliMCInfo * cinfo = GetInfo(label); |
711 | if (cinfo) cinfo->CalcTPCrows(runArrayTR); | |
de6c9df4 | 712 | // |
713 | // TPC | |
714 | // | |
715 | if (trackRef->DetectorId()== AliTrackReference::kTPC){ | |
716 | // | |
de6c9df4 | 717 | Int_t label = trackRef->GetTrack(); |
718 | AliMCInfo * info = GetInfo(label); | |
999d8278 | 719 | if (!info && trackRef->Pt()<fTPCPtCut) continue; |
d390cc7e | 720 | if (!info) { |
721 | info = MakeInfo(label); | |
722 | } | |
de6c9df4 | 723 | if (!info) continue; |
724 | info->fPrimPart = iPrimPart; | |
725 | TClonesArray & arr = *(info->fTPCReferences); | |
726 | new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); | |
727 | } | |
728 | // | |
729 | // ITS | |
730 | // | |
731 | if (trackRef->DetectorId()== AliTrackReference::kITS){ | |
732 | // | |
de6c9df4 | 733 | Int_t label = trackRef->GetTrack(); |
734 | AliMCInfo * info = GetInfo(label); | |
999d8278 | 735 | if (!info && trackRef->Pt()<fITSPtCut) continue; |
de6c9df4 | 736 | if (!info) info = MakeInfo(label); |
737 | if (!info) continue; | |
738 | info->fPrimPart = iPrimPart; | |
739 | TClonesArray & arr = *(info->fITSReferences); | |
740 | new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); | |
741 | } | |
742 | // | |
743 | // TRD | |
744 | // | |
745 | if (trackRef->DetectorId()== AliTrackReference::kTRD){ | |
746 | // | |
de6c9df4 | 747 | Int_t label = trackRef->GetTrack(); |
748 | AliMCInfo * info = GetInfo(label); | |
999d8278 | 749 | if (!info&&trackRef->P()<fTRDPtCut) continue; |
de6c9df4 | 750 | if (!info) info = MakeInfo(label); |
751 | if (!info) continue; | |
752 | info->fPrimPart = iPrimPart; | |
753 | TClonesArray & arr = *(info->fTRDReferences); | |
754 | new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); | |
755 | } | |
756 | // | |
757 | // TOF | |
758 | // | |
759 | if (trackRef->DetectorId()== AliTrackReference::kTOF){ | |
760 | // | |
de6c9df4 | 761 | Int_t label = trackRef->GetTrack(); |
762 | AliMCInfo * info = GetInfo(label); | |
999d8278 | 763 | if (!info&&trackRef->P()<fTOFPtCut) continue; |
de6c9df4 | 764 | if (!info) info = MakeInfo(label); |
765 | if (!info) continue; | |
766 | info->fPrimPart = iPrimPart; | |
767 | TClonesArray & arr = *(info->fTOFReferences); | |
768 | new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef); | |
769 | } | |
770 | // | |
771 | // decay case | |
772 | // | |
773 | if (trackRef->DetectorId()== AliTrackReference::kDisappeared){ | |
774 | // | |
de6c9df4 | 775 | AliMCInfo * info = GetInfo(label); |
776 | if (!info) continue; | |
777 | info->fTRdecay = *trackRef; | |
778 | } | |
de6c9df4 | 779 | } |
780 | } | |
781 | // | |
782 | // | |
d92975ba | 783 | return 0; |
784 | } | |
785 | ||
786 | //////////////////////////////////////////////////////////////////////// | |
787 | Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef, | |
788 | AliTPCParam *paramTPC) const { | |
789 | ||
9f3282ed | 790 | // |
791 | // Transformation from Gloabal to local tacking system | |
792 | // | |
d92975ba | 793 | Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()}; |
794 | Int_t index[4]; | |
795 | paramTPC->Transform0to1(x,index); | |
796 | paramTPC->Transform1to2(x,index); | |
797 | return x[0]; | |
798 | } | |
799 | //////////////////////////////////////////////////////////////////////// | |
800 | ||
801 | ||
802 | ||
803 | AliTPCParam * AliGenInfoMaker::GetTPCParam(){ | |
9f3282ed | 804 | // |
805 | // create default TPC parameters | |
806 | // | |
d92975ba | 807 | AliTPCParamSR * par = new AliTPCParamSR; |
808 | par->Update(); | |
809 | return par; | |
810 | } | |
811 | ||
812 |