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