2 // This class is the task for connecting together
3 // MC information and the RC information
5 // The task is a wrapper over two components
7 // AliESDRecInfoMaker.h
14 #include <TTreeStream.h>
15 #include <AliAnalysisManager.h>
16 #include <AliESDInputHandler.h>
18 #include "AliMCEvent.h"
19 #include "AliMCEventHandler.h"
22 #include "AliGenInfoTask.h"
23 #include "AliGenInfoMaker.h"
27 #include "AliMCInfo.h"
28 #include "AliComparisonObject.h"
29 #include "AliESDRecInfo.h"
30 #include "AliTPCParamSR.h"
32 #include "TTimeStamp.h"
34 #include "AliTPCseed.h"
41 ClassImp(AliGenInfoTask)
43 //________________________________________________________________________
44 AliGenInfoTask::AliGenInfoTask() :
46 fMCinfo(0), //! MC event handler
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
60 // Default constructor (should not be used)
64 AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
66 fMCinfo(0), //! MC event handler
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
81 // Default constructor
87 //________________________________________________________________________
88 AliGenInfoTask::AliGenInfoTask(const char *name) :
89 AliAnalysisTask(name, "AliGenInfoTask"),
90 fMCinfo(0), //! MC event handler
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
104 // Normal constructor
106 // Input slot #0 works with a TChain
107 DefineInput(0, TChain::Class());
108 // Output slot #0 writes into a TList
109 DefineOutput(0, TObjArray::Class());
112 fCompList = new TObjArray;
115 AliGenInfoTask::~AliGenInfoTask(){
119 if (fDebugLevel>0) printf("AliGenInfoTask::~AliGenInfoTask\n");
120 if (fDebugStreamer) delete fDebugStreamer;
122 if(fCompList) delete fCompList;
127 //________________________________________________________________________
128 void AliGenInfoTask::ConnectInputData(Option_t *)
131 // Connect the input data
134 cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
136 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
138 //Printf("ERROR: Could not read chain from input slot 0");
141 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
143 //Printf("ERROR: Could not get ESDInputHandler");
146 fESD = esdH->GetEvent();
147 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
148 fESD->SetESDfriend(fESDfriend);
149 //Printf("*** CONNECTED NEW EVENT ****");
152 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
153 mcinfo->SetReadTR(kTRUE);
155 fMCinfo = mcinfo->MCEvent();
161 //_____________________________________________________________________________
162 Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj)
164 // add comparison object to the list
166 Printf("ERROR: Could not add comparison object");
169 // add object to the list
170 fCompList->AddLast(pObj);
177 //________________________________________________________________________
178 void AliGenInfoTask::CreateOutputObjects()
181 // Connect the output objects
184 cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
189 //________________________________________________________________________
190 void AliGenInfoTask::Exec(Option_t *) {
192 // Execute analysis for current event
196 cout << "AliGenInfoTask::Exec()" << endl;
199 // If MC has been connected
200 if (fGenTracksArray) fGenTracksArray->Delete();
201 if (fRecTracksArray) fRecTracksArray->Delete();
204 cout << "Not MC info\n" << endl;
218 //________________________________________________________________________
219 void AliGenInfoTask::Terminate(Option_t *) {
224 printf("AliGenInfoTask: Terminate() \n");
226 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
227 if (fDebugStreamer) delete fDebugStreamer;
236 TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){
238 // Get Debug streamer
239 // In case debug streamer not yet initialized and StreamLevel>0 create new one
241 if (fStreamLevel==0) return 0;
242 if (fDebugStreamer) return fDebugStreamer;
245 dsName+="Debug.root";
246 dsName.ReplaceAll(" ","");
247 fDebugStreamer = new TTreeSRedirector(dsName.Data());
248 return fDebugStreamer;
253 AliMCInfo* AliGenInfoTask::GetTrack(Int_t index, Bool_t force){
255 // Get the MC info for given track
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;
262 info = new ((*fGenTracksArray)[index]) AliMCInfo;
267 AliESDRecInfo* AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){
269 // Get the MC info for given track
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;
276 info = new ((*fRecTracksArray)[index]) AliESDRecInfo;
284 void AliGenInfoTask::ProcessMCInfo(){
286 // Dump information from MC to the array
289 TParticle * particle= new TParticle;
290 TClonesArray * trefs = new TClonesArray("AliTrackReference");
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);
300 vertex[0]=particle->Vx();
301 vertex[1]=particle->Vy();
302 vertex[2]=particle->Vz();
303 vertex[3]=particle->R();
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);
316 TTreeSRedirector *pcstream = GetDebugStreamer();
326 void AliGenInfoTask::ProcessESDInfo(){
330 static AliTPCParamSR param;
331 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
333 //Printf("ERROR: fESDfriend not available");
336 fESD->SetESDfriend(fESDfriend);
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);
356 void AliGenInfoTask::ProcessComparison(){
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;
368 for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){
369 AliComparisonObject *pObj= (AliComparisonObject *)fCompList->At(icomp);
371 pObj->Exec(mcinfo,recInfo);
375 PostData(0, fCompList);
378 void AliGenInfoTask::DumpInfo(){
382 TParticle * particle= new TParticle;
383 TClonesArray * trefs = new TClonesArray("AliTrackReference");
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();
395 fMCinfo->GetParticleAndTR(ipart, particle, trefs);
398 if (trefs!=0 && particle!=0){
399 length = GetTPCTrackLength(*trefs, particle , fESD->GetMagneticField(), counter, 3.0);
407 "counter="<<counter<<
408 // "tr.="<<tpctrack<<
416 Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
420 TCut cutPt("p.Pt()>0.1");
421 TCut cutZ("abs(p.Vz())<250");
422 TCut cutR("abs(p.R())<250");
427 if (part->Pt()<0.1) return kFALSE;
428 if (TMath::Abs(part->Vz())>250) return kFALSE;
429 if (part->R()>360) return kFALSE;
431 if (part->GetPDG()->Charge()==0) return kFALSE;
438 void AliGenInfoTask::FinishTaskOutput()
441 // According description in AliAnalisysTask this method is call
442 // on the slaves before sending data
445 RegisterDebugOutput(fDebugOutputPath.Data());
451 void AliGenInfoTask::RegisterDebugOutput(const char */*path*/){
453 // store - copy debug output to the destination position
454 // currently ONLY for local copy
457 dsName+="Debug.root";
458 dsName.ReplaceAll(" ","");
459 TString dsName2=fDebugOutputPath.Data();
460 gSystem->MakeDirectory(dsName2.Data());
462 dsName2+=gSystem->HostName();
463 gSystem->MakeDirectory(dsName2.Data());
465 dsName2+=gSystem->BaseName(gSystem->pwd());
467 gSystem->MakeDirectory(dsName2.Data());
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());
475 Float_t AliGenInfoTask::GetTPCTrackLength(const TClonesArray& trackRefs, TParticle*part, Float_t bz, Int_t &counter, Float_t deadWidth){
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
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;
491 Int_t nrefs = trackRefs.GetEntriesFast();
492 AliExternalTrackParam param;
494 for (Int_t i=0; i<21;i++) cv[i]=0;
498 AliTrackReference *ref0 = (AliTrackReference*)trackRefs.At(0);
499 Float_t direction = 0;
501 for (Int_t iref=1; iref<nrefs;iref++){
502 AliTrackReference *ref = (AliTrackReference*)trackRefs.At(iref);
504 if (!ref0 || ref0->DetectorId()!= AliTrackReference::kTPC){
506 direction = ((ref0->X()*ref0->Px()+ref0->Y()*ref0->Py())>0)? 1.:-1.;
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;
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()};
522 param.Set(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
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();
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++;
537 if (ref->DetectorId()!= AliTrackReference::kTPC) break;