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