]>
Commit | Line | Data |
---|---|---|
13db274d | 1 | // |
2 | // This class is the task for connecting together | |
3 | // MC information and the RC information | |
4 | // | |
5 | // The task is a wrapper over two components | |
6 | // AliGenInfoMaker | |
7 | // AliESDRecInfoMaker.h | |
8 | ||
9 | // ROOT includes | |
10 | #include <TChain.h> | |
11 | #include <TMath.h> | |
12 | ||
13 | // ALIROOT includes | |
9c3fd353 | 14 | #include <TTreeStream.h> |
13db274d | 15 | #include <AliAnalysisManager.h> |
16 | #include <AliESDInputHandler.h> | |
17 | #include "AliStack.h" | |
18 | #include "AliMCEvent.h" | |
19 | #include "AliMCEventHandler.h" | |
20 | ||
21 | #include <AliESD.h> | |
22 | #include "AliGenInfoTask.h" | |
23 | #include "AliGenInfoMaker.h" | |
24 | #include "AliHelix.h" | |
25 | ||
9c3fd353 | 26 | // |
27 | #include "AliMCInfo.h" | |
28 | #include "AliComparisonObject.h" | |
29 | #include "AliESDRecInfo.h" | |
30 | #include "AliTPCParamSR.h" | |
95cba04d | 31 | #include "TSystem.h" |
32 | #include "TTimeStamp.h" | |
33 | #include "TFile.h" | |
34 | #include "AliTPCseed.h" | |
9c3fd353 | 35 | |
13db274d | 36 | // STL includes |
37 | #include <iostream> | |
38 | ||
39 | using namespace std; | |
40 | ||
41 | ClassImp(AliGenInfoTask) | |
42 | ||
43 | //________________________________________________________________________ | |
44 | AliGenInfoTask::AliGenInfoTask() : | |
cd875161 | 45 | AliAnalysisTask(), |
9c3fd353 | 46 | fMCinfo(0), //! MC event handler |
47 | fESD(0), | |
252a6266 | 48 | fESDfriend(0), |
9c3fd353 | 49 | fCompList(0), //array of comparison objects |
50 | fGenTracksArray(0), //clones array with filtered particles | |
51 | fGenKinkArray(0), //clones array with filtered Kinks | |
52 | fGenV0Array(0), //clones array with filtered V0s | |
53 | fRecTracksArray(0), //clones array with filtered particles | |
54 | fDebugStreamer(0), | |
55 | fStreamLevel(0), | |
95cba04d | 56 | fDebugLevel(0), |
57 | fDebugOutputPath() | |
13db274d | 58 | { |
59 | // | |
60 | // Default constructor (should not be used) | |
61 | // | |
13db274d | 62 | } |
63 | ||
cd875161 | 64 | AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) : |
65 | AliAnalysisTask(), | |
9c3fd353 | 66 | fMCinfo(0), //! MC event handler |
67 | fESD(0), | |
252a6266 | 68 | fESDfriend(0), |
9c3fd353 | 69 | fCompList(0), |
70 | fGenTracksArray(0), //clones array with filtered particles | |
71 | fGenKinkArray(0), //clones array with filtered Kinks | |
72 | fGenV0Array(0), //clones array with filtered V0s | |
73 | fRecTracksArray(0), //clones array with filtered particles | |
74 | // | |
75 | fDebugStreamer(0), | |
76 | fStreamLevel(0), | |
95cba04d | 77 | fDebugLevel(), |
78 | fDebugOutputPath() | |
cd875161 | 79 | { |
80 | // | |
81 | // Default constructor | |
82 | // | |
cd875161 | 83 | } |
84 | ||
85 | ||
86 | ||
13db274d | 87 | //________________________________________________________________________ |
88 | AliGenInfoTask::AliGenInfoTask(const char *name) : | |
89 | AliAnalysisTask(name, "AliGenInfoTask"), | |
9c3fd353 | 90 | fMCinfo(0), //! MC event handler |
91 | fESD(0), | |
252a6266 | 92 | fESDfriend(0), |
9c3fd353 | 93 | fCompList(0), |
94 | fGenTracksArray(0), //clones array with filtered particles | |
95 | fGenKinkArray(0), //clones array with filtered Kinks | |
96 | fGenV0Array(0), //clones array with filtered V0s | |
97 | fRecTracksArray(0), //clones array with filtered particles | |
98 | fDebugStreamer(0), | |
99 | fStreamLevel(0), | |
95cba04d | 100 | fDebugLevel(0), |
101 | fDebugOutputPath() | |
13db274d | 102 | { |
103 | // | |
104 | // Normal constructor | |
105 | // | |
13db274d | 106 | // Input slot #0 works with a TChain |
107 | DefineInput(0, TChain::Class()); | |
108 | // Output slot #0 writes into a TList | |
9c3fd353 | 109 | DefineOutput(0, TObjArray::Class()); |
110 | // | |
111 | // | |
112 | fCompList = new TObjArray; | |
113 | } | |
114 | ||
115 | AliGenInfoTask::~AliGenInfoTask(){ | |
116 | // | |
117 | // | |
118 | // | |
119 | if (fDebugLevel>0) printf("AliGenInfoTask::~AliGenInfoTask\n"); | |
120 | if (fDebugStreamer) delete fDebugStreamer; | |
121 | fDebugStreamer=0; | |
122 | if(fCompList) delete fCompList; | |
123 | fCompList =0; | |
13db274d | 124 | } |
125 | ||
9c3fd353 | 126 | |
13db274d | 127 | //________________________________________________________________________ |
128 | void AliGenInfoTask::ConnectInputData(Option_t *) | |
129 | { | |
130 | // | |
131 | // Connect the input data | |
132 | // | |
9c3fd353 | 133 | if(fDebugLevel>3) |
13db274d | 134 | cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl; |
135 | ||
9c3fd353 | 136 | TTree* tree=dynamic_cast<TTree*>(GetInputData(0)); |
137 | if (!tree) { | |
138 | //Printf("ERROR: Could not read chain from input slot 0"); | |
139 | } | |
140 | else { | |
141 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
142 | if (!esdH) { | |
143 | //Printf("ERROR: Could not get ESDInputHandler"); | |
144 | } | |
145 | else { | |
146 | fESD = esdH->GetEvent(); | |
95cba04d | 147 | fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend")); |
148 | fESD->SetESDfriend(fESDfriend); | |
9c3fd353 | 149 | //Printf("*** CONNECTED NEW EVENT ****"); |
150 | } | |
151 | } | |
152 | AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
153 | mcinfo->SetReadTR(kTRUE); | |
154 | ||
155 | fMCinfo = mcinfo->MCEvent(); | |
13db274d | 156 | |
13db274d | 157 | |
9c3fd353 | 158 | } |
159 | ||
160 | ||
161 | //_____________________________________________________________________________ | |
162 | Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj) | |
163 | { | |
164 | // add comparison object to the list | |
165 | if(pObj == 0) { | |
166 | Printf("ERROR: Could not add comparison object"); | |
167 | return kFALSE; | |
13db274d | 168 | } |
9c3fd353 | 169 | // add object to the list |
170 | fCompList->AddLast(pObj); | |
171 | return kTRUE; | |
13db274d | 172 | } |
173 | ||
9c3fd353 | 174 | |
175 | ||
176 | ||
13db274d | 177 | //________________________________________________________________________ |
178 | void AliGenInfoTask::CreateOutputObjects() | |
179 | { | |
180 | // | |
181 | // Connect the output objects | |
182 | // | |
9c3fd353 | 183 | if(fDebugLevel>3) |
13db274d | 184 | cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl; |
185 | ||
13db274d | 186 | } |
187 | ||
188 | ||
189 | //________________________________________________________________________ | |
190 | void AliGenInfoTask::Exec(Option_t *) { | |
191 | // | |
192 | // Execute analysis for current event | |
13db274d | 193 | // |
194 | ||
9c3fd353 | 195 | if(fDebugLevel>3) |
13db274d | 196 | cout << "AliGenInfoTask::Exec()" << endl; |
9c3fd353 | 197 | |
13db274d | 198 | |
13db274d | 199 | // If MC has been connected |
9c3fd353 | 200 | if (fGenTracksArray) fGenTracksArray->Delete(); |
201 | if (fRecTracksArray) fRecTracksArray->Delete(); | |
202 | ||
203 | if (!fMCinfo){ | |
13db274d | 204 | cout << "Not MC info\n" << endl; |
9c3fd353 | 205 | }else{ |
206 | //mcinfo->Print(); | |
207 | ProcessMCInfo(); | |
208 | ProcessESDInfo(); | |
209 | DumpInfo(); | |
210 | ProcessComparison(); | |
13db274d | 211 | } |
9c3fd353 | 212 | // |
213 | } | |
13db274d | 214 | |
215 | ||
216 | ||
9c3fd353 | 217 | |
218 | //________________________________________________________________________ | |
219 | void AliGenInfoTask::Terminate(Option_t *) { | |
220 | // | |
221 | // Terminate loop | |
222 | // | |
223 | if(fDebugLevel>3) | |
224 | printf("AliGenInfoTask: Terminate() \n"); | |
225 | // | |
226 | if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n"); | |
227 | if (fDebugStreamer) delete fDebugStreamer; | |
228 | fDebugStreamer = 0; | |
229 | return; | |
230 | } | |
231 | ||
232 | ||
233 | ||
95cba04d | 234 | |
235 | ||
9c3fd353 | 236 | TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){ |
237 | // | |
238 | // Get Debug streamer | |
239 | // In case debug streamer not yet initialized and StreamLevel>0 create new one | |
240 | // | |
241 | if (fStreamLevel==0) return 0; | |
242 | if (fDebugStreamer) return fDebugStreamer; | |
243 | TString dsName; | |
244 | dsName=GetName(); | |
245 | dsName+="Debug.root"; | |
246 | dsName.ReplaceAll(" ",""); | |
247 | fDebugStreamer = new TTreeSRedirector(dsName.Data()); | |
248 | return fDebugStreamer; | |
249 | } | |
250 | ||
251 | ||
252 | ||
253 | AliMCInfo* AliGenInfoTask::GetTrack(Int_t index, Bool_t force){ | |
254 | // | |
255 | // Get the MC info for given track | |
256 | // | |
257 | if (!fGenTracksArray) fGenTracksArray = new TClonesArray("AliMCInfo",1000); | |
258 | if (index>fGenTracksArray->GetEntriesFast()) fGenTracksArray->Expand(index*2+10); | |
259 | AliMCInfo * info = (AliMCInfo*)fGenTracksArray->At(index); | |
260 | if (!force) return info; | |
261 | if (!info){ | |
262 | info = new ((*fGenTracksArray)[index]) AliMCInfo; | |
13db274d | 263 | } |
9c3fd353 | 264 | return info; |
265 | } | |
266 | ||
267 | AliESDRecInfo* AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){ | |
268 | // | |
269 | // Get the MC info for given track | |
270 | // | |
271 | if (!fRecTracksArray) fRecTracksArray = new TClonesArray("AliESDRecInfo",1000); | |
272 | if (index>fRecTracksArray->GetEntriesFast()) fRecTracksArray->Expand(index*2+10); | |
273 | AliESDRecInfo * info = (AliESDRecInfo*)fRecTracksArray->At(index); | |
274 | if (!force) return info; | |
275 | if (!info){ | |
276 | info = new ((*fRecTracksArray)[index]) AliESDRecInfo; | |
13db274d | 277 | } |
9c3fd353 | 278 | return info; |
279 | } | |
13db274d | 280 | |
9c3fd353 | 281 | |
282 | ||
283 | ||
284 | void AliGenInfoTask::ProcessMCInfo(){ | |
285 | // | |
286 | // Dump information from MC to the array | |
287 | // | |
288 | // | |
289 | TParticle * particle= new TParticle; | |
290 | TClonesArray * trefs = new TClonesArray("AliTrackReference"); | |
13db274d | 291 | // |
13db274d | 292 | // |
9c3fd353 | 293 | // Process tracks |
294 | // | |
295 | Int_t npart = fMCinfo->GetNumberOfTracks(); | |
296 | if (npart==0) return; | |
297 | Double_t vertex[4]={0,0,0,0}; | |
298 | fMCinfo->GetParticleAndTR(0, particle, trefs); | |
299 | if (particle){ | |
300 | vertex[0]=particle->Vx(); | |
301 | vertex[1]=particle->Vy(); | |
302 | vertex[2]=particle->Vz(); | |
303 | vertex[3]=particle->R(); | |
304 | } | |
13db274d | 305 | |
9c3fd353 | 306 | for (Int_t ipart=0;ipart<npart;ipart++){ |
307 | Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs); | |
308 | if (status<0) continue; | |
309 | if (!particle) continue; | |
310 | if (!trefs) continue; | |
311 | if (!AcceptParticle(particle)) continue; | |
312 | //if (trefs->GetEntries()<1) continue; | |
313 | AliMCInfo * mcinfo = GetTrack(ipart,kTRUE); | |
314 | mcinfo->Update(particle,trefs,vertex,ipart); | |
315 | // | |
316 | TTreeSRedirector *pcstream = GetDebugStreamer(); | |
317 | if (pcstream){ | |
318 | (*pcstream)<<"MC"<< | |
319 | "p.="<<particle<< | |
320 | "MC.="<<mcinfo<< | |
321 | "\n"; | |
322 | } | |
323 | } | |
13db274d | 324 | } |
325 | ||
9c3fd353 | 326 | void AliGenInfoTask::ProcessESDInfo(){ |
327 | // | |
328 | // | |
329 | // | |
330 | static AliTPCParamSR param; | |
95cba04d | 331 | fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend")); |
332 | if (!fESDfriend) { | |
333 | //Printf("ERROR: fESDfriend not available"); | |
334 | return; | |
335 | } | |
336 | fESD->SetESDfriend(fESDfriend); | |
9c3fd353 | 337 | // |
338 | // | |
339 | if (!fESD) return; | |
340 | Int_t ntracks = fESD->GetNumberOfTracks(); | |
341 | for (Int_t itrack=0; itrack<ntracks; itrack++){ | |
342 | AliESDtrack *track = fESD->GetTrack(itrack); | |
343 | Int_t label = TMath::Abs(track->GetLabel()); | |
344 | AliMCInfo * mcinfo = GetTrack(label,kFALSE); | |
345 | if (!mcinfo) continue; | |
346 | AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE); | |
347 | recInfo->AddESDtrack(track,mcinfo); | |
348 | recInfo->Update(mcinfo,¶m,kTRUE); | |
349 | } | |
350 | // | |
9c3fd353 | 351 | |
352 | ||
353 | } | |
354 | ||
355 | ||
356 | void AliGenInfoTask::ProcessComparison(){ | |
357 | // | |
358 | // | |
359 | // | |
360 | static AliESDRecInfo dummy; | |
361 | Int_t npart = fMCinfo->GetNumberOfTracks(); | |
362 | for (Int_t ipart=0;ipart<npart;ipart++){ | |
363 | AliMCInfo * mcinfo = GetTrack(ipart,kFALSE); | |
364 | if (!mcinfo) continue; | |
365 | AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE); | |
366 | if (!recInfo) recInfo=&dummy; | |
367 | // | |
368 | for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){ | |
369 | AliComparisonObject *pObj= (AliComparisonObject *)fCompList->At(icomp); | |
370 | if (pObj){ | |
371 | pObj->Exec(mcinfo,recInfo); | |
372 | } | |
373 | } | |
374 | } | |
375 | PostData(0, fCompList); | |
376 | } | |
377 | ||
378 | void AliGenInfoTask::DumpInfo(){ | |
379 | // | |
380 | // | |
381 | // | |
76421d44 | 382 | TParticle * particle= new TParticle; |
383 | TClonesArray * trefs = new TClonesArray("AliTrackReference"); | |
384 | // | |
95cba04d | 385 | |
9c3fd353 | 386 | static AliESDRecInfo dummy; |
387 | Int_t npart = fMCinfo->GetNumberOfTracks(); | |
388 | for (Int_t ipart=0;ipart<npart;ipart++){ | |
389 | AliMCInfo * mcinfo = GetTrack(ipart,kFALSE); | |
390 | if (!mcinfo) continue; | |
391 | AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE); | |
392 | if (!recInfo) recInfo=&dummy; | |
393 | TTreeSRedirector *pcstream = GetDebugStreamer(); | |
76421d44 | 394 | |
395 | fMCinfo->GetParticleAndTR(ipart, particle, trefs); | |
396 | Int_t counter=0; | |
397 | Float_t length=0; | |
398 | if (trefs!=0 && particle!=0){ | |
399 | length = GetTPCTrackLength(*trefs, particle , fESD->GetMagneticField(), counter, 3.0); | |
400 | } | |
401 | ||
9c3fd353 | 402 | if (pcstream){ |
95cba04d | 403 | (*pcstream)<<"RC"<< |
9c3fd353 | 404 | "MC.="<<mcinfo<< |
405 | "RC.="<<recInfo<< | |
76421d44 | 406 | "length="<<length<< |
407 | "counter="<<counter<< | |
95cba04d | 408 | // "tr.="<<tpctrack<< |
9c3fd353 | 409 | "\n"; |
410 | } | |
411 | } | |
13db274d | 412 | } |
9c3fd353 | 413 | |
414 | ||
415 | ||
416 | Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){ | |
417 | // | |
418 | /* | |
419 | MC cuts | |
420 | TCut cutPt("p.Pt()>0.1"); | |
421 | TCut cutZ("abs(p.Vz())<250"); | |
422 | TCut cutR("abs(p.R())<250"); | |
423 | // | |
424 | */ | |
425 | // | |
426 | // | |
427 | if (part->Pt()<0.1) return kFALSE; | |
428 | if (TMath::Abs(part->Vz())>250) return kFALSE; | |
429 | if (part->R()>360) return kFALSE; | |
430 | if (part->GetPDG()){ | |
431 | if (part->GetPDG()->Charge()==0) return kFALSE; | |
432 | } | |
433 | return kTRUE; | |
434 | } | |
435 | ||
95cba04d | 436 | |
437 | ||
438 | void AliGenInfoTask::FinishTaskOutput() | |
439 | { | |
440 | // | |
441 | // According description in AliAnalisysTask this method is call | |
442 | // on the slaves before sending data | |
443 | // | |
444 | Terminate("slave"); | |
445 | RegisterDebugOutput(fDebugOutputPath.Data()); | |
446 | ||
447 | } | |
448 | ||
449 | ||
450 | ||
252a6266 | 451 | void AliGenInfoTask::RegisterDebugOutput(const char */*path*/){ |
95cba04d | 452 | // |
453 | // store - copy debug output to the destination position | |
454 | // currently ONLY for local copy | |
455 | TString dsName; | |
456 | dsName=GetName(); | |
457 | dsName+="Debug.root"; | |
458 | dsName.ReplaceAll(" ",""); | |
459 | TString dsName2=fDebugOutputPath.Data(); | |
460 | gSystem->MakeDirectory(dsName2.Data()); | |
461 | dsName2+="/";; | |
462 | dsName2+=gSystem->HostName(); | |
463 | gSystem->MakeDirectory(dsName2.Data()); | |
464 | dsName2+="/"; | |
465 | dsName2+=gSystem->BaseName(gSystem->pwd()); | |
466 | dsName2+="/"; | |
467 | gSystem->MakeDirectory(dsName2.Data()); | |
468 | dsName2+=dsName; | |
469 | AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data())); | |
470 | printf("copy %s\t%s\n",dsName.Data(),dsName2.Data()); | |
471 | TFile::Cp(dsName.Data(),dsName2.Data()); | |
472 | } | |
76421d44 | 473 | |
474 | ||
475 | Float_t AliGenInfoTask::GetTPCTrackLength(const TClonesArray& trackRefs, TParticle*part, Float_t bz, Int_t &counter, Float_t deadWidth){ | |
476 | // | |
477 | // return track length in geometrically active volume of TPC. | |
478 | // z nad rphi acceptance is included | |
479 | // doesn't take into account dead channel and ExB | |
480 | // Intput: | |
481 | // trackRefs | |
482 | // bz - magnetic field | |
483 | // deadWidth - dead zone in r-phi | |
484 | // Additional output: | |
485 | // counter - number of circles | |
486 | const Float_t kRMin = 90; | |
487 | const Float_t kRMax = 245; | |
488 | const Float_t kZMax = 250; | |
489 | const Float_t kMinPt= 0.01; | |
490 | Float_t length =0; | |
491 | Int_t nrefs = trackRefs.GetEntriesFast(); | |
492 | AliExternalTrackParam param; | |
493 | Double_t cv[21]; | |
494 | for (Int_t i=0; i<21;i++) cv[i]=0; | |
495 | counter=0; | |
496 | // | |
497 | // | |
498 | AliTrackReference *ref0 = (AliTrackReference*)trackRefs.At(0); | |
499 | Float_t direction = 0; | |
500 | // | |
501 | for (Int_t iref=1; iref<nrefs;iref++){ | |
502 | AliTrackReference *ref = (AliTrackReference*)trackRefs.At(iref); | |
503 | if (!ref) continue; | |
504 | if (!ref0 || ref0->DetectorId()!= AliTrackReference::kTPC){ | |
505 | ref0 = ref; | |
506 | direction = ((ref0->X()*ref0->Px()+ref0->Y()*ref0->Py())>0)? 1.:-1.; | |
507 | continue; | |
508 | } | |
509 | Float_t newdirection = ((ref->X()*ref->Px()+ref->Y()*ref->Py())>0)? 1.:-1.; | |
510 | if (newdirection*direction<0) { | |
511 | counter++; //circle counter | |
512 | direction = newdirection; | |
513 | continue; | |
514 | } | |
515 | if (counter>0) continue; | |
516 | if (ref0->Pt()<kMinPt) break; | |
517 | Float_t radius0 = TMath::Max(TMath::Min(ref0->R(),kRMax),kRMin);; | |
518 | Float_t radius1 = TMath::Max(TMath::Min(ref->R(),kRMax),kRMin); | |
519 | Double_t xyz[3]={ref0->X(),ref0->Y(),ref0->Z()}; | |
520 | Double_t pxyz[3]={ref0->Px(),ref0->Py(),ref0->Pz()}; | |
521 | Double_t alpha; | |
522 | param.Set(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.)); | |
523 | ||
524 | for (Float_t radius = radius0; radius<radius1; radius+=1){ | |
525 | param.GetXYZAt(radius, bz, xyz); | |
526 | if (TMath::Abs(xyz[2])>kZMax) continue; | |
527 | Float_t gradius = TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]); | |
528 | if (gradius>kRMax) continue; | |
529 | alpha = TMath::ATan2(xyz[1],xyz[0]); | |
530 | if (alpha<0) alpha+=TMath::TwoPi(); | |
531 | // | |
532 | Int_t sector = Int_t(9*alpha/TMath::Pi()); | |
533 | Float_t lalpha = alpha-((sector+0.5)*TMath::Pi()/9.); | |
534 | Float_t dedge = (TMath::Tan(TMath::Pi()/18.)-TMath::Abs(TMath::Tan(lalpha)))*gradius; | |
535 | if (dedge>deadWidth) length++; | |
536 | } | |
537 | if (ref->DetectorId()!= AliTrackReference::kTPC) break; | |
538 | ref0 = ref; | |
539 | } | |
540 | return length; | |
541 | } |