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