PWGPP/TPC/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 #include "TSystem.h"
  32 #include "TTimeStamp.h"
  33 #include "TFile.h"
  34 #include "AliTPCseed.h"
  35
  36 // STL includes
  37 #include <iostream>
  38
  39 using namespace std;
  40
  41 ClassImp(AliGenInfoTask)
  42
  43 //________________________________________________________________________
  44 AliGenInfoTask::AliGenInfoTask() :
  45   AliAnalysisTask(),
  46   fMCinfo(0),     //! MC event handler
  47   fESD(0),
  48   fESDfriend(0),
  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),
  56   fDebugLevel(0),
  57   fDebugOutputPath()
  58 {
  59   //
  60   // Default constructor (should not be used)
  61   //
  62 }
  63
  64 AliGenInfoTask::AliGenInfoTask(const AliGenInfoTask& /*info*/) :
  65   AliAnalysisTask(),
  66   fMCinfo(0),     //! MC event handler
  67   fESD(0),
  68   fESDfriend(0),
  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),
  77   fDebugLevel(),
  78   fDebugOutputPath()
  79 {
  80   //
  81   // Default constructor
  82   //
  83 }
  84
  85
  86
  87 //________________________________________________________________________
  88 AliGenInfoTask::AliGenInfoTask(const char *name) :
  89   AliAnalysisTask(name, "AliGenInfoTask"),
  90   fMCinfo(0),     //! MC event handler
  91   fESD(0),
  92   fESDfriend(0),
  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),
 100   fDebugLevel(0),
 101   fDebugOutputPath()
 102 {
 103   //
 104   // Normal constructor
 105   //
 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());
 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;
 124 }
 125
 126
 127 //________________________________________________________________________
 128 void AliGenInfoTask::ConnectInputData(Option_t *)
 129 {
 130   //
 131   // Connect the input data
 132   //
 133   if(fDebugLevel>3)
 134     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 135
 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();
 147       fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 148       fESD->SetESDfriend(fESDfriend);
 149       //Printf("*** CONNECTED NEW EVENT ****");
 150     }
 151   }
 152   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 153   mcinfo->SetReadTR(kTRUE);
 154
 155   fMCinfo = mcinfo->MCEvent();
 156
 157
 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;
 168   }
 169   // add object to the list
 170   fCompList->AddLast(pObj);
 171   return kTRUE;
 172 }
 173
 174
 175
 176
 177 //________________________________________________________________________
 178 void AliGenInfoTask::CreateOutputObjects()
 179 {
 180   //
 181   // Connect the output objects
 182   //
 183   if(fDebugLevel>3)
 184     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 185
 186 }
 187
 188
 189 //________________________________________________________________________
 190 void AliGenInfoTask::Exec(Option_t *) {
 191   //
 192   // Execute analysis for current event
 193   //
 194
 195   if(fDebugLevel>3)
 196     cout << "AliGenInfoTask::Exec()" << endl;
 197
 198
 199   // If MC has been connected
 200   if (fGenTracksArray) fGenTracksArray->Delete();
 201   if (fRecTracksArray) fRecTracksArray->Delete();
 202
 203   if (!fMCinfo){
 204     cout << "Not MC info\n" << endl;
 205   }else{
 206     //mcinfo->Print();
 207     ProcessMCInfo();
 208     ProcessESDInfo();
 209     DumpInfo();
 210     ProcessComparison();
 211   }
 212   //
 213 }
 214
 215
 216
 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
 234
 235
 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;
 263   }
 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;
 277   }
 278   return info;
 279 }
 280
 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");
 291   //
 292   //
 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   }
 305
 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   }
 324 }
 325
 326 void AliGenInfoTask::ProcessESDInfo(){
 327   //
 328   //
 329   //
 330   if (!fESD) return;
 331   static AliTPCParamSR param;
 332   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 333   if (!fESDfriend) {
 334     //Printf("ERROR: fESDfriend not available");
 335     return;
 336   }
 337   fESD->SetESDfriend(fESDfriend);
 338   //
 339   //
 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,&param,kTRUE);
 349   }
 350   //
 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   //
 382   TParticle * particle= new TParticle;
 383   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 384   //
 385
 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();
 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
 402     if (pcstream){
 403       (*pcstream)<<"RC"<<
 404         "MC.="<<mcinfo<<
 405         "RC.="<<recInfo<<
 406         "length="<<length<<
 407         "counter="<<counter<<
 408         //      "tr.="<<tpctrack<<
 409         "\n";
 410     }
 411   }
 412 }
 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
 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
 451 void AliGenInfoTask::RegisterDebugOutput(const char */*path*/){
 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 }
 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 }