ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx

   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 //  Analysis task for Kinematic filtering
  18 //  Fill AODtrackMC tracks from Kinematic stack
  19 //
  20
  21 #include <TChain.h>
  22 #include <TFile.h>
  23 #include "TParticle.h"
  24 #include "TParameter.h"
  25 #include "TString.h"
  26 #include "TList.h"
  27 #include "TProfile.h"
  28 #include "TH1F.h"
  29
  30
  31 #include "AliAnalysisTaskMCParticleFilter.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliAnalysisFilter.h"
  34 #include "AliHeader.h"
  35 #include "AliStack.h"
  36 #include "AliMCEvent.h"
  37 #include "AliMCEventHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODHeader.h"
  40 #include "AliAODMCHeader.h"
  41 #include "AliAODHandler.h"
  42 #include "AliAODVertex.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliCollisionGeometry.h"
  45 #include "AliGenDPMjetEventHeader.h"
  46 #include "AliGenHijingEventHeader.h"
  47 #include "AliGenPythiaEventHeader.h"
  48 #include "AliGenCocktailEventHeader.h"
  49
  50 #include "AliLog.h"
  51
  52 ClassImp(AliAnalysisTaskMCParticleFilter)
  53
  54 ////////////////////////////////////////////////////////////////////////
  55
  56 //____________________________________________________________________
  57 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
  58   AliAnalysisTaskSE(),
  59   fTrackFilterMother(0x0),
  60   fAODMcHeader(0x0),
  61   fAODMcParticles(0x0),
  62   fHistList(0x0)
  63 {
  64   // Default constructor
  65 }
  66
  67 Bool_t AliAnalysisTaskMCParticleFilter::Notify()
  68 {
  69   //
  70   // Implemented Notify() to read the cross sections
  71   // from pyxsec.root
  72   //
  73   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
  74   Double_t xsection = 0;
  75   UInt_t   ntrials  = 0;
  76   if(tree){
  77     TFile *curfile = tree->GetCurrentFile();
  78     if (!curfile) {
  79       Error("Notify","No current file");
  80       return kFALSE;
  81     }
  82
  83     TString fileName(curfile->GetName());
  84     if(fileName.Contains("AliESDs.root")){
  85         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
  86     }
  87     else if(fileName.Contains("AliAOD.root")){
  88         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
  89     }
  90     else if(fileName.Contains("galice.root")){
  91         // for running with galice and kinematics alone...
  92         fileName.ReplaceAll("galice.root", "pyxsec.root");
  93     }
  94
  95
  96     TFile *fxsec = TFile::Open(fileName.Data());
  97     if(!fxsec){
  98       AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
  99       // not a severe condition we just do not have the information...
 100       return kTRUE;
 101     }
 102     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 103     if(!xtree){
 104       AliWarning(Form("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__));
 105       return kTRUE;
 106     }
 107     xtree->SetBranchAddress("xsection",&xsection);
 108     xtree->SetBranchAddress("ntrials",&ntrials);
 109     xtree->GetEntry(0);
 110     ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
 111   }
 112   return kTRUE;
 113 }
 114
 115
 116
 117 //____________________________________________________________________
 118 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
 119     AliAnalysisTaskSE(name),
 120     fTrackFilterMother(0x0),
 121     fAODMcHeader(0x0),
 122     fAODMcParticles(0x0),
 123     fHistList(0x0)
 124 {
 125   // Default constructor
 126   DefineOutput(1, TList::Class());
 127 }
 128
 129 /*
 130 //____________________________________________________________________
 131 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
 132     AliAnalysisTaskSE(obj),
 133     fTrackFilterMother(obj.fTrackFilterMother)
 134 {
 135   // Copy constructor
 136 }
 137 */
 138 //____________________________________________________________________
 139 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
 140 {
 141
 142   if(fAODMcHeader){
 143     delete fAODMcHeader;
 144   }
 145   if(fAODMcParticles){
 146     fAODMcParticles->Delete();
 147     delete fAODMcParticles;
 148   }
 149 }
 150
 151 /*
 152 //____________________________________________________________________
 153 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
 154 {
 155 // Assignment
 156   if(this!=&other) {
 157     AliAnalysisTaskSE::operator=(other);
 158     fTrackFilterMother = other.fTrackFilterMother;
 159   }
 160   return *this;
 161 }
 162 */
 163 //____________________________________________________________________
 164 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
 165 {
 166   // Create the output container
 167
 168
 169     if (OutputTree()&&fTrackFilterMother)
 170         OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
 171
 172     // this part is mainly needed to set the MCEventHandler
 173     // to the AODHandler, this will not be needed when
 174     // AODHandler goes to ANALYSISalice
 175     // setting in the steering macro will not work on proof :(
 176     // so we have to do it in a task
 177
 178     // the branch booking can also go into the AODHandler then
 179
 180
 181     // mcparticles
 182     fAODMcParticles = new TClonesArray("AliAODMCParticle", 0);
 183     fAODMcParticles->SetName(AliAODMCParticle::StdBranchName());
 184     AddAODBranch("TClonesArray",&fAODMcParticles);
 185
 186     // MC header...
 187     fAODMcHeader = new AliAODMCHeader();
 188     fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
 189     AddAODBranch("AliAODMCHeader",&fAODMcHeader);
 190
 191     AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 192     AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 193     if(!aodH){
 194       AliWarning("Could not get AODHandler");
 195       return;
 196     }
 197     aodH->SetMCEventHandler(mcH);
 198
 199
 200     // these are histograms, for averaging and summing
 201     // do it via histograms to be PROOF-proof
 202     // which merges the results from different workers
 203     // histos are not saved reside only in memory
 204
 205
 206
 207     fHistList = new TList();
 208     TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
 209     h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 210     fHistList->Add(h1Xsec);
 211     TH1F* h1Trials = new TH1F("h1Trials","trials from MC header",1,0,1);
 212     h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 213     fHistList->Add(h1Trials);
 214
 215     PostData(1,fHistList);
 216
 217 }
 218
 219 //____________________________________________________________________
 220 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
 221 {
 222   // Execute analysis for current event
 223   //
 224
 225   // Fill AOD tracks from Kinematic stack
 226   PostData(1,fHistList);
 227
 228   // get AliAOD Event
 229   AliAODEvent* aod = AODEvent();
 230   if (!aod) {
 231       AliWarning("No Output Handler connected, doing nothing !") ;
 232       return;
 233   }
 234
 235   AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 236   if(!mcH){
 237     AliWarning("No MC handler Found");
 238     return;
 239   }
 240
 241   AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 242   if(!aodH){
 243     AliWarning("Could not get AODHandler");
 244     return;
 245   }
 246
 247
 248
 249   // fetch the output
 250   // Fill control histos
 251
 252   // tmp array for holding the mctracks
 253   // need to set mother and duagther __before__
 254   // filling in the tree
 255
 256   AliMCEvent *mcE = MCEvent();
 257
 258   if(!mcE){
 259     AliWarning("No MC event Found");
 260     return;
 261   }
 262
 263   Int_t np    = mcE->GetNumberOfTracks();
 264   Int_t nprim = mcE->GetNumberOfPrimaries();
 265   // TODO ADD MC VERTEX
 266
 267   // Get the proper MC Collision Geometry
 268   AliGenEventHeader* mcEH = mcE->GenEventHeader();
 269
 270   AliGenPythiaEventHeader *pyH  = dynamic_cast<AliGenPythiaEventHeader*>(mcEH);
 271   AliGenHijingEventHeader *hiH  = 0;
 272   AliCollisionGeometry    *colG = 0;
 273   AliGenDPMjetEventHeader *dpmH = 0;
 274   // it can be only one save some casts
 275   // assuming PYTHIA and HIJING are the most likely ones...
 276   if(!pyH){
 277       hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
 278       if(!hiH){
 279           dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
 280       }
 281   }
 282
 283   if (hiH || dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);
 284
 285   // fetch the trials on a event by event basis, not from pyxsec.root otherwise
 286   // we will get a problem when running on proof since Notify may be called
 287   // more than once per file
 288   // consider storing this information in the AOD output via AliAODHandler
 289   Float_t ntrials = 0;
 290   if (!colG) {
 291     AliGenCocktailEventHeader *ccEH = dynamic_cast<AliGenCocktailEventHeader *>(mcEH);
 292     if (ccEH) {
 293       TList *genHeaders = ccEH->GetHeaders();
 294       for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
 295         if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
 296         if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
 297         if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
 298         if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
 299       }
 300     }
 301   }
 302
 303   // take the trials from the p+p event
 304   if(hiH)ntrials = hiH->Trials();
 305   if(dpmH)ntrials = dpmH->Trials();
 306   if(pyH)ntrials = pyH->Trials();
 307   if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
 308
 309
 310
 311
 312   if (colG) {
 313     fAODMcHeader->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
 314     AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
 315   }
 316
 317
 318
 319
 320   Int_t j=0;
 321   for (Int_t ip = 0; ip < np; ip++){
 322     AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
 323     TParticle* part = mcpart->Particle();
 324     Float_t xv = part->Vx();
 325     Float_t yv = part->Vy();
 326     Float_t zv = part->Vz();
 327     Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
 328
 329     Bool_t write = kFALSE;
 330
 331     if (ip < nprim) {
 332       // Select the primary event
 333       write = kTRUE;
 334     } else if (part->GetUniqueID() == kPDecay) {
 335       // Particles from decay
 336       // Check that the decay chain ends at a primary particle
 337       AliMCParticle* mother = mcpart;
 338       Int_t imo = mcpart->GetMother();
 339       while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) {
 340         mother =  (AliMCParticle*) mcE->GetTrack(imo);
 341         imo =  mother->GetMother();
 342       }
 343       // Select according to pseudorapidity and production point of primary ancestor
 344       if (imo < nprim)write = kTRUE;
 345       // if(!Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kFALSE; // selection on eta and phi of the mother
 346     } else if (part->GetUniqueID() == kPPair) {
 347       // Now look for pair production
 348       Int_t imo = mcpart->GetMother();
 349       if (imo < nprim) {
 350         // Select, if the gamma is a primary
 351         write = kTRUE;
 352       } else {
 353         // Check if the gamma comes from the decay chain of a primary particle
 354         AliMCParticle* mother =  (AliMCParticle*) mcE->GetTrack(imo);
 355         imo = mother->GetMother();
 356         while((imo >= nprim) && (mother->GetUniqueID() == kPDecay)) {
 357           mother =   (AliMCParticle*) mcE->GetTrack(imo);
 358           imo =  mother->GetMother();
 359         }
 360         // Select according to pseudorapidity and production point
 361         if (imo < nprim && Select(mother->Particle(), rv, zv))
 362           write = kTRUE;
 363       }
 364     }
 365
 366     if (write) {
 367       if(mcH)mcH->SelectParticle(ip);
 368       j++;
 369     }
 370   }
 371
 372   /*
 373
 374   // check varibales for charm need all daughters
 375   static int  iTaken = 0;
 376   static int  iAll = 0;
 377   static int  iCharm = 0;
 378
 379   for (Int_t ip = 0; ip < np; ip++){
 380     AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
 381     TParticle* part = mcpart->Particle();
 382     // debug info to check fro charm daugthers
 383     if((TMath::Abs(part->GetPdgCode()))==411){
 384       iCharm++;
 385       Int_t d0 =  part->GetFirstDaughter();
 386       Int_t d1 =  part->GetLastDaughter();
 387       if(d0>0&&d1>0){
 388         for(int id = d0;id <= d1;id++){
 389           iAll++;
 390           if(mcH->IsParticleSelected(id)){
 391             iTaken++;
 392             Printf("> %d/%d Taken",id,nprim);
 393             PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
 394           }
 395           else{
 396             Printf("> %d/%d NOT Taken",id,nprim);
 397             PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
 398           }
 399         }
 400       }
 401     }// if charm
 402     AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
 403   }
 404   */
 405
 406
 407   aodH->StoreMCParticles();
 408
 409   return;
 410 }
 411
 412 void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
 413   //
 414   // Terminate the execution save the PYTHIA x_section to the UserInfo()
 415   //
 416
 417
 418   fHistList = (TList*)GetOutputData(1);
 419   if(!fHistList){
 420     Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
 421     return;
 422   }
 423
 424   TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
 425   TH1F* hTrials  = ((TH1F*)(fHistList->FindObject("h1Trials")));
 426   if(!hXsec)return;
 427   if(!hTrials)return;
 428
 429   Float_t xsec = hXsec->GetBinContent(1);
 430   Float_t trials = hTrials->Integral();
 431   AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
 432
 433 }
 434
 435 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
 436 {
 437   // Selection accoring to eta of the mother and production point
 438   // This has to be refined !!!!!!
 439
 440   // Esp if we don't have collisison in the central barrel but e.g. beam gas
 441   //  return kTRUE;
 442
 443   Float_t eta = part->Eta();
 444
 445   // central barrel consider everything in the ITS...
 446   // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
 447   // larger for V0s in the TPC
 448   //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
 449
 450   if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
 451
 452   // Andreas' Cuts
 453   //  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
 454
 455
 456
 457   // Muon arm
 458   if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
 459
 460   // PMD acceptance 2.3 <= eta < = 3.5
 461   //  if(eta>2.0&&eta<3.7)return kTRUE;
 462
 463   return kFALSE;
 464
 465 }
 466
 467 void AliAnalysisTaskMCParticleFilter::PrintMCParticle(const AliMCParticle *mcp,Int_t np){
 468
 469   const TParticle *p = mcp->Particle();
 470   Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughte\
 471 r2 %d ",
 472          np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter \
 473          (0),p->GetDaughter(1));
 474   Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
 475   p->Print();
 476   Printf("---------------------------------------");
 477 }