PWG1/AliGenInfoTask.cxx

   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
  14 #include <TTreeStream.h>
  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
  26 //
  27 #include "AliMCInfo.h"
  28 #include "AliComparisonObject.h"
  29 #include "AliESDRecInfo.h"
  30 #include "AliTPCParamSR.h"
  31
  32 // STL includes
  33 #include <iostream>
  34
  35 using namespace std;
  36
  37 ClassImp(AliGenInfoTask)
  38
  39 //________________________________________________________________________
  40 AliGenInfoTask::AliGenInfoTask() :
  41   AliAnalysisTask(),
  42   fMCinfo(0),     //! MC event handler
  43   fESD(0),
  44   fCompList(0),         //array of comparison objects
  45   fGenTracksArray(0),  //clones array with filtered particles
  46   fGenKinkArray(0),    //clones array with filtered Kinks
  47   fGenV0Array(0),      //clones array with filtered V0s
  48   fRecTracksArray(0),  //clones array with filtered particles
  49   fDebugStreamer(0),
  50   fStreamLevel(0),
  51   fDebugLevel(0)
  52 {
  53   //
  54   // Default constructor (should not be used)
  55   //
  56 }
  57
  58 AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
  59   AliAnalysisTask(),
  60   fMCinfo(0),     //! MC event handler
  61   fESD(0),
  62   fCompList(0),
  63   fGenTracksArray(0),  //clones array with filtered particles
  64   fGenKinkArray(0),    //clones array with filtered Kinks
  65   fGenV0Array(0),      //clones array with filtered V0s
  66   fRecTracksArray(0),  //clones array with filtered particles
  67   //
  68   fDebugStreamer(0),
  69   fStreamLevel(0),
  70   fDebugLevel()
  71 {
  72   //
  73   // Default constructor
  74   //
  75 }
  76
  77
  78
  79 //________________________________________________________________________
  80 AliGenInfoTask::AliGenInfoTask(const char *name) :
  81   AliAnalysisTask(name, "AliGenInfoTask"),
  82   fMCinfo(0),     //! MC event handler
  83   fESD(0),
  84   fCompList(0),
  85   fGenTracksArray(0),  //clones array with filtered particles
  86   fGenKinkArray(0),    //clones array with filtered Kinks
  87   fGenV0Array(0),      //clones array with filtered V0s
  88   fRecTracksArray(0),  //clones array with filtered particles
  89   fDebugStreamer(0),
  90   fStreamLevel(0),
  91   fDebugLevel(0)
  92 {
  93   //
  94   // Normal constructor
  95   //
  96   // Input slot #0 works with a TChain
  97   DefineInput(0, TChain::Class());
  98   // Output slot #0 writes into a TList
  99   DefineOutput(0, TObjArray::Class());
 100   //
 101   //
 102   fCompList = new TObjArray;
 103 }
 104
 105 AliGenInfoTask::~AliGenInfoTask(){
 106   //
 107   //
 108   //
 109   if (fDebugLevel>0)  printf("AliGenInfoTask::~AliGenInfoTask\n");
 110   if (fDebugStreamer) delete fDebugStreamer;
 111   fDebugStreamer=0;
 112   if(fCompList)   delete fCompList;
 113   fCompList =0;
 114 }
 115
 116
 117 //________________________________________________________________________
 118 void AliGenInfoTask::ConnectInputData(Option_t *)
 119 {
 120   //
 121   // Connect the input data
 122   //
 123   if(fDebugLevel>3)
 124     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 125
 126   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 127   if (!tree) {
 128     //Printf("ERROR: Could not read chain from input slot 0");
 129   }
 130   else {
 131     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 132     if (!esdH) {
 133       //Printf("ERROR: Could not get ESDInputHandler");
 134     }
 135     else {
 136       fESD = esdH->GetEvent();
 137       //Printf("*** CONNECTED NEW EVENT ****");
 138     }
 139   }
 140   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 141   mcinfo->SetReadTR(kTRUE);
 142
 143   fMCinfo = mcinfo->MCEvent();
 144
 145
 146 }
 147
 148
 149 //_____________________________________________________________________________
 150 Bool_t AliGenInfoTask::AddComparisonObject(AliComparisonObject *pObj)
 151 {
 152   // add comparison object to the list
 153   if(pObj == 0) {
 154       Printf("ERROR: Could not add comparison object");
 155           return kFALSE;
 156   }
 157   // add object to the list
 158   fCompList->AddLast(pObj);
 159   return kTRUE;
 160 }
 161
 162
 163
 164
 165 //________________________________________________________________________
 166 void AliGenInfoTask::CreateOutputObjects()
 167 {
 168   //
 169   // Connect the output objects
 170   //
 171   if(fDebugLevel>3)
 172     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 173
 174 }
 175
 176
 177 //________________________________________________________________________
 178 void AliGenInfoTask::Exec(Option_t *) {
 179   //
 180   // Execute analysis for current event
 181   //
 182
 183   if(fDebugLevel>3)
 184     cout << "AliGenInfoTask::Exec()" << endl;
 185
 186
 187   // If MC has been connected
 188   if (fGenTracksArray) fGenTracksArray->Delete();
 189   if (fRecTracksArray) fRecTracksArray->Delete();
 190
 191   if (!fMCinfo){
 192     cout << "Not MC info\n" << endl;
 193   }else{
 194     //mcinfo->Print();
 195     ProcessMCInfo();
 196     ProcessESDInfo();
 197     DumpInfo();
 198     ProcessComparison();
 199   }
 200   //
 201 }
 202
 203
 204
 205
 206 //________________________________________________________________________
 207 void AliGenInfoTask::Terminate(Option_t *) {
 208     //
 209     // Terminate loop
 210     //
 211   if(fDebugLevel>3)
 212     printf("AliGenInfoTask: Terminate() \n");
 213   //
 214   if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
 215   if (fDebugStreamer) delete fDebugStreamer;
 216   fDebugStreamer = 0;
 217   return;
 218 }
 219
 220
 221
 222 TTreeSRedirector *AliGenInfoTask::GetDebugStreamer(){
 223   //
 224   // Get Debug streamer
 225   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 226   //
 227   if (fStreamLevel==0) return 0;
 228   if (fDebugStreamer) return fDebugStreamer;
 229   TString dsName;
 230   dsName=GetName();
 231   dsName+="Debug.root";
 232   dsName.ReplaceAll(" ","");
 233   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 234   return fDebugStreamer;
 235 }
 236
 237
 238
 239 AliMCInfo*  AliGenInfoTask::GetTrack(Int_t index, Bool_t force){
 240   //
 241   // Get the MC info for given track
 242   //
 243   if (!fGenTracksArray) fGenTracksArray = new TClonesArray("AliMCInfo",1000);
 244   if (index>fGenTracksArray->GetEntriesFast()) fGenTracksArray->Expand(index*2+10);
 245   AliMCInfo * info = (AliMCInfo*)fGenTracksArray->At(index);
 246   if (!force) return info;
 247   if (!info){
 248     info = new ((*fGenTracksArray)[index]) AliMCInfo;
 249   }
 250   return info;
 251 }
 252
 253 AliESDRecInfo*  AliGenInfoTask::GetRecTrack(Int_t index, Bool_t force){
 254   //
 255   // Get the MC info for given track
 256   //
 257   if (!fRecTracksArray) fRecTracksArray = new TClonesArray("AliESDRecInfo",1000);
 258   if (index>fRecTracksArray->GetEntriesFast()) fRecTracksArray->Expand(index*2+10);
 259   AliESDRecInfo * info = (AliESDRecInfo*)fRecTracksArray->At(index);
 260   if (!force) return info;
 261   if (!info){
 262     info = new ((*fRecTracksArray)[index]) AliESDRecInfo;
 263   }
 264   return info;
 265 }
 266
 267
 268
 269
 270 void  AliGenInfoTask::ProcessMCInfo(){
 271   //
 272   // Dump information from MC to the array
 273   //
 274   //
 275   TParticle * particle= new TParticle;
 276   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 277   //
 278   //
 279   // Process tracks
 280   //
 281   Int_t npart = fMCinfo->GetNumberOfTracks();
 282   if (npart==0) return;
 283   Double_t vertex[4]={0,0,0,0};
 284   fMCinfo->GetParticleAndTR(0, particle, trefs);
 285   if (particle){
 286     vertex[0]=particle->Vx();
 287     vertex[1]=particle->Vy();
 288     vertex[2]=particle->Vz();
 289     vertex[3]=particle->R();
 290   }
 291
 292   for (Int_t ipart=0;ipart<npart;ipart++){
 293     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 294     if (status<0) continue;
 295     if (!particle) continue;
 296     if (!trefs) continue;
 297     if (!AcceptParticle(particle)) continue;
 298     //if (trefs->GetEntries()<1) continue;
 299     AliMCInfo * mcinfo = GetTrack(ipart,kTRUE);
 300     mcinfo->Update(particle,trefs,vertex,ipart);
 301     //
 302     TTreeSRedirector *pcstream = GetDebugStreamer();
 303     if (pcstream){
 304       (*pcstream)<<"MC"<<
 305         "p.="<<particle<<
 306         "MC.="<<mcinfo<<
 307         "\n";
 308     }
 309   }
 310 }
 311
 312 void AliGenInfoTask::ProcessESDInfo(){
 313   //
 314   //
 315   //
 316   static AliTPCParamSR param;
 317   //
 318   //
 319   if (!fESD) return;
 320   Int_t ntracks = fESD->GetNumberOfTracks();
 321   for (Int_t itrack=0; itrack<ntracks; itrack++){
 322     AliESDtrack *track = fESD->GetTrack(itrack);
 323     Int_t label = TMath::Abs(track->GetLabel());
 324     AliMCInfo * mcinfo = GetTrack(label,kFALSE);
 325     if (!mcinfo) continue;
 326     AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
 327     recInfo->AddESDtrack(track,mcinfo);
 328     recInfo->Update(mcinfo,&param,kTRUE);
 329   }
 330   //
 331   //
 332   //
 333   Int_t ntracksMC = fMCinfo->GetNumberOfTracks();
 334   for (Int_t imc=0; imc<ntracksMC; imc++){
 335     AliMCInfo * mcinfo = GetTrack(imc,kFALSE);
 336     if (!mcinfo) continue;
 337     AliESDRecInfo *recInfo= GetRecTrack(imc,kFALSE);
 338     if (recInfo) continue;
 339     if (mcinfo->GetNTPCRef()<2) continue;
 340     //
 341     //
 342     for (Int_t itrack=0; itrack<ntracks; itrack++){
 343       AliESDtrack *track = fESD->GetTrack(itrack);
 344       Int_t label = TMath::Abs(track->GetLabel());
 345       if (label!=mcinfo->GetLabel()) continue;
 346
 347       AliMCInfo * mcinfo2 = GetTrack(label,kFALSE);
 348       if (!mcinfo2) continue;
 349       AliESDRecInfo *recInfo= GetRecTrack(label,kTRUE);
 350       recInfo->AddESDtrack(track,mcinfo2);
 351       recInfo->Update(mcinfo2,&param,kTRUE);
 352     }
 353   }
 354
 355
 356
 357
 358 }
 359
 360
 361 void AliGenInfoTask::ProcessComparison(){
 362   //
 363   //
 364   //
 365   static AliESDRecInfo dummy;
 366   Int_t npart = fMCinfo->GetNumberOfTracks();
 367   for (Int_t ipart=0;ipart<npart;ipart++){
 368     AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
 369     if (!mcinfo) continue;
 370     AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
 371     if (!recInfo) recInfo=&dummy;
 372     //
 373     for (Int_t icomp = 0; icomp<fCompList->GetEntries(); icomp++){
 374       AliComparisonObject *pObj= (AliComparisonObject *)fCompList->At(icomp);
 375       if (pObj){
 376         pObj->Exec(mcinfo,recInfo);
 377       }
 378     }
 379   }
 380   PostData(0, fCompList);
 381 }
 382
 383 void AliGenInfoTask::DumpInfo(){
 384   //
 385   //
 386   //
 387   static AliESDRecInfo dummy;
 388   Int_t npart = fMCinfo->GetNumberOfTracks();
 389   for (Int_t ipart=0;ipart<npart;ipart++){
 390     AliMCInfo * mcinfo = GetTrack(ipart,kFALSE);
 391     if (!mcinfo) continue;
 392     AliESDRecInfo *recInfo= GetRecTrack(ipart,kFALSE);
 393     if (!recInfo) recInfo=&dummy;
 394     TTreeSRedirector *pcstream = GetDebugStreamer();
 395     if (pcstream){
 396       (*pcstream)<<"RC"<<
 397         "MC.="<<mcinfo<<
 398         "RC.="<<recInfo<<
 399         "\n";
 400     }
 401   }
 402 }
 403
 404
 405
 406 Bool_t AliGenInfoTask::AcceptParticle(TParticle *part){
 407   //
 408   /*
 409     MC cuts
 410     TCut cutPt("p.Pt()>0.1");
 411     TCut cutZ("abs(p.Vz())<250");
 412     TCut cutR("abs(p.R())<250");
 413     //
 414   */
 415   //
 416   //
 417   if (part->Pt()<0.1) return kFALSE;
 418   if (TMath::Abs(part->Vz())>250) return kFALSE;
 419   if (part->R()>360)  return kFALSE;
 420   if (part->GetPDG()){
 421     if (part->GetPDG()->Charge()==0) return kFALSE;
 422   }
 423   return kTRUE;
 424 }
 425