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