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