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 "TString.h"
  25
  26 #include "AliAnalysisTaskMCParticleFilter.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliAnalysisFilter.h"
  29 #include "AliHeader.h"
  30 #include "AliStack.h"
  31 #include "AliMCEvent.h"
  32 #include "AliMCEventHandler.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODHeader.h"
  35 #include "AliAODHandler.h"
  36 #include "AliAODVertex.h"
  37 #include "AliAODMCParticle.h"
  38
  39
  40 #include "AliLog.h"
  41
  42 ClassImp(AliAnalysisTaskMCParticleFilter)
  43
  44 ////////////////////////////////////////////////////////////////////////
  45
  46 //____________________________________________________________________
  47 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
  48 AliAnalysisTaskSE(),
  49   fTrackFilterMother(0x0)
  50 {
  51   // Default constructor
  52 }
  53
  54 //____________________________________________________________________
  55 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
  56     AliAnalysisTaskSE(name),
  57     fTrackFilterMother(0x0)
  58 {
  59   // Default constructor
  60 }
  61
  62
  63 //____________________________________________________________________
  64 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
  65     AliAnalysisTaskSE(obj),
  66     fTrackFilterMother(obj.fTrackFilterMother)
  67 {
  68   // Copy constructor
  69 }
  70
  71 //____________________________________________________________________
  72 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
  73 {
  74   //  if( fTrackFilterMother ) delete fTrackFilterMother;
  75 }
  76
  77
  78 //____________________________________________________________________
  79 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
  80 {
  81 // Assignment
  82   if(this!=&other) {
  83     AliAnalysisTaskSE::operator=(other);
  84     fTrackFilterMother = other.fTrackFilterMother;
  85   }
  86   return *this;
  87 }
  88
  89 //____________________________________________________________________
  90 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
  91 {
  92   // Create the output container
  93     if (OutputTree()&&fTrackFilterMother)
  94         OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
  95
  96     // this part is mainly needed to set the MCEventHandler
  97     // to the AODHandler, this will not be needed when
  98     // AODHandler goes to ANALYSISalice
  99     // setting in the steering macro will not work on proof :(
 100     // so we have to do it in a task
 101
 102     // the branch booking can also go into the AODHandler then
 103
 104     TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
 105     tca->SetName(AliAODMCParticle::StdBranchName());
 106     AddAODBranch("TClonesArray",&tca);
 107     AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 108
 109     AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 110     if(!aodH){
 111       Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
 112       return;
 113     }
 114     aodH->SetMCEventHandler(mcH);
 115     // TODO ADD MC VERTEX
 116
 117 }
 118
 119 //____________________________________________________________________
 120 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
 121 {
 122 // Execute analysis for current event
 123 //
 124
 125 // Fill AOD tracks from Kinematic stack
 126
 127   // get AliAOD Event
 128   AliAODEvent* aod = AODEvent();
 129   if (!aod) {
 130       AliWarning("No Output Handler connected, doing nothing !") ;
 131       return;
 132   }
 133
 134   AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 135   if(!mcH){
 136     AliWarning("No MC handler Found");
 137     return;
 138   }
 139
 140
 141   // fetch the output
 142   // Fill control histos
 143
 144   // tmp array for holding the mctracks
 145   // need to set mother and duagther __before__
 146   // filling in the tree
 147
 148   AliMCEvent *mcE = MCEvent();
 149
 150   if(!mcE){
 151     AliWarning("No MC event Found");
 152     return;
 153   }
 154
 155   AliStack* stack = mcE->Stack();
 156   Int_t np    = mcE->GetNumberOfTracks();
 157   Int_t nprim = stack->GetNprimary();
 158   // TODO ADD MC VERTEX
 159
 160   // We take all real primaries
 161
 162
 163
 164   Int_t j=0;
 165   for (Int_t ip = 0; ip < np; ip++){
 166     AliMCParticle* mcpart =  mcE->GetTrack(ip);
 167     TParticle* part = mcpart->Particle();
 168
 169     Float_t xv = part->Vx();
 170     Float_t yv = part->Vy();
 171     Float_t zv = part->Vz();
 172     Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
 173
 174     Bool_t write = kFALSE;
 175     Int_t flag = 0;
 176
 177
 178     if (ip < nprim) {
 179       // Select the primary event
 180       write = kTRUE;
 181     } else if (part->GetUniqueID() == 4) {
 182       // Particles from decay
 183       // Check that the decay chain ends at a primary particle
 184       TParticle* mother = part;
 185       Int_t imo = part->GetFirstMother();
 186       while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
 187         mother =  mcE->Stack()->Particle(imo);
 188         imo =  mother->GetFirstMother();
 189       }
 190       // Select according to pseudorapidity and production point
 191         if (imo < nprim && Select(mother, rv, zv))
 192           write = kTRUE;
 193     } else if (part->GetUniqueID() == 5) {
 194       // Now look for pair production
 195       continue;
 196       Int_t imo = part->GetFirstMother();
 197       if (imo < nprim) {
 198         // Select, if the gamma is a primary
 199         write = kTRUE;
 200       } else {
 201         // Check if the gamma comes from the decay chain of a primary particle
 202         TParticle* mother = mcE->Stack()->Particle(imo);
 203         imo = mother->GetFirstMother();
 204         while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
 205           mother =  mcE->Stack()->Particle(imo);
 206           imo =  mother->GetFirstMother();
 207         }
 208         // Select according to pseudorapidity and production point
 209         if (imo < nprim && Select(mother, rv, zv))
 210           write = kTRUE;
 211       }
 212     }
 213     /*
 214     else if (part->GetUniqueID() == 13){
 215       // Evaporation
 216       // Check that we end at a primary particle
 217       TParticle* mother = part;
 218       Int_t imo = part->GetFirstMother();
 219       while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
 220         mother =  mcE->Stack()->Particle(imo);
 221         imo =  mother->GetFirstMother();
 222       }
 223       // Select according to pseudorapidity and production point
 224         if (imo < nprim && Select(mother, rv, zv))
 225           write = kTRUE;
 226     }
 227     */
 228     if (write) {
 229       if(mcH)mcH->SelectParticle(ip);
 230       j++;
 231     }
 232   }
 233
 234   return;
 235 }
 236
 237 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
 238 {
 239   // Selection accoring to eta of the mother and production point
 240   // This has to be refined !!!!!!
 241
 242   // Esp if we don't have collisison in the central barrel but e.g. beam gas
 243   //  return kTRUE;
 244
 245   Float_t eta = part->Eta();
 246
 247   // central barrel consider everything in the ITS...
 248   // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
 249   // larger for V0s in the TPC
 250   //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
 251
 252   // Andreas' Cuts
 253   if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
 254
 255   // Muon arm
 256   if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
 257
 258   // PMD acceptance 2.3 <= eta < = 3.5
 259   //  if(eta>2.0&&eta<3.7)return kTRUE;
 260
 261   return kFALSE;
 262
 263 }
 264