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