PWGLF/totEt/AliAnalysisEtSelector.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Selector Base class
   4 //  -
   5 // implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen)
   8 //_________________________________________________________________________
   9 #include "AliAnalysisEtSelector.h"
  10 #include "AliAnalysisEtCuts.h"
  11 #include "AliESDCaloCluster.h"
  12 #include "AliStack.h"
  13 #include "TParticle.h"
  14 #include "TParticlePDG.h"
  15 #include "AliAnalysisEtCommon.h"
  16 #include <iostream>
  17
  18 ClassImp(AliAnalysisEtSelector);
  19
  20 AliAnalysisEtSelector::AliAnalysisEtSelector(AliAnalysisEtCuts *cuts) : AliAnalysisEtCommon()
  21 ,fEvent(0)
  22 ,fClusterArray(0)
  23 ,fCuts(cuts)
  24 ,fRunNumber(0)
  25 ,fInitialized(kFALSE)
  26 {
  27 }
  28 AliAnalysisEtSelector::AliAnalysisEtSelector() : AliAnalysisEtCommon()
  29 ,fEvent(0)
  30 ,fClusterArray(0)
  31 ,fCuts(0)
  32 ,fRunNumber(0)
  33 ,fInitialized(kFALSE)
  34 {
  35   fCuts = new AliAnalysisEtCuts();
  36 }
  37
  38 AliAnalysisEtSelector::~AliAnalysisEtSelector()
  39 { // dtor
  40   if(fClusterArray)
  41   {
  42     delete fClusterArray;
  43   }
  44 }
  45 void AliAnalysisEtSelector::SetEvent(const AliESDEvent* event)
  46 {
  47     fEvent = event;
  48     if(!fInitialized) Init(event);
  49 }
  50
  51 Bool_t AliAnalysisEtSelector::IsNeutralMcParticle(Int_t pIdx, AliStack& s, const TParticlePDG& pdg) const
  52 {
  53   return s.IsPhysicalPrimary(pIdx) &&(TMath::Abs(TMath::Abs(pdg.Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3);
  54 }
  55
  56 Bool_t AliAnalysisEtSelector::IsEmEtParticle(const Int_t pdgCode) const
  57 {
  58   return pdgCode == fgGammaCode || pdgCode == fgPi0Code || pdgCode == fgEtaCode || pdgCode == fgEPlusCode || pdgCode == fgEMinusCode;
  59 }
  60
  61
  62 Bool_t AliAnalysisEtSelector::PrimaryIsEmEtParticle(const Int_t pIdx, AliStack& stack) const
  63 {
  64   return IsEmEtParticle(stack.Particle(GetPrimary(pIdx, stack))->GetPdgCode());
  65 }
  66 Int_t AliAnalysisEtSelector::GetPrimary(const Int_t partIdx, AliStack& stack) const
  67 { // get primary
  68   if(partIdx >= 0)
  69   {
  70     Int_t mothIdx = stack.Particle(partIdx)->GetMother(0);
  71     if(mothIdx < 0) return -1;
  72     TParticle *mother = stack.Particle(mothIdx);
  73     if(mother)
  74     {
  75       if(stack.IsPhysicalPrimary(mothIdx)) return mothIdx;
  76       else return GetPrimary(mothIdx, stack);
  77     }
  78     else
  79     {
  80       std::cout << "WAT!" << std::endl;
  81       return -1;
  82     }
  83   }
  84   return -1;
  85 }
  86 //Bool_t AliAnalysisEtSelector::FromSecondaryInteraction(const TParticle& part, AliStack &stack) const
  87 Bool_t AliAnalysisEtSelector::FromSecondaryInteraction(Int_t partID, AliStack &stack) const
  88 {
  89 //   Bool_t partVtxSecondary = (
  90 //                            TMath::Sqrt(part.Vx()*part.Vx() + part.Vy()*part.Vy()) > fCuts->GetPrimaryVertexCutXY()
  91 //                            || TMath::Abs(part.Vz()) > fCuts->GetPrimaryVertexCutZ()
  92 //                          )
  93 //                          && TMath::Sqrt(part.Vx()*part.Vx()+part.Vy()*part.Vy() + part.Vz()*part.Vz())<(fCuts->GetGeometryPhosDetectorRadius()-10);
  94
  95
  96   //Bool_t partVtxSecondary = (TMath::Sqrt(part.Vx()*part.Vx() + part.Vy()*part.Vy()) <420);
  97   //if(partVtxSecondary) return kFalse;
  98 //   //Let's find suspect decay (typical for secondary interaction)...
  99 //   if(partVtxSecondary){
 100   // return SuspiciousDecayInChain(211, 111, part, stack);
 101   //return stack.IsSecondaryFromMaterial(part.GetUniqueID());
 102   return stack.IsSecondaryFromMaterial(partID);
 103 //   }
 104 //   else{
 105 //     return kFALSE;
 106 //   }
 107
 108
 109
 110
 111 }
 112
 113 Int_t AliAnalysisEtSelector::GetMother(Int_t partID, AliStack& stack) const {
 114   if(partID>0){
 115     TParticle *particle = stack.Particle(partID);
 116       if(particle){
 117
 118         return particle->GetMother(0);
 119       }
 120   }
 121   return -1;
 122 }
 123 Bool_t AliAnalysisEtSelector::IsFromDetectorCover(Int_t partID, AliStack& stack) const{
 124   if(partID>0){
 125     TParticle *particle = stack.Particle(partID);
 126     if(particle){//particle exists
 127       if(stack.IsSecondaryFromMaterial(partID)){//particle is from an interaction with the material
 128         //say that it's from the detector cover if its vertex is larger than 400 cm because this is where the cover is
 129         return (TMath::Sqrt(particle->Vx()*particle->Vx() + particle->Vy()*particle->Vy()) >400);
 130
 131       }
 132     }
 133   }
 134   return kFALSE;
 135 }
 136
 137 Int_t AliAnalysisEtSelector::GetFirstMotherNotFromDetectorCover(Int_t partID, AliStack& stack) const{
 138   Int_t targetParticle = partID;
 139   Int_t testParticle = partID;
 140   Int_t iteration  = 0;
 141   while(IsFromDetectorCover(targetParticle,stack) && testParticle>0){//if the particle is from the detector cover
 142     //cout<<"Particle "<<targetParticle<<" is from detector cover.  Getting mother ";
 143     testParticle = GetMother(targetParticle,stack);
 144     if(testParticle>0) targetParticle = testParticle;
 145     //cout<<targetParticle<<" iteration "<<iteration<<endl;
 146     iteration++;
 147   }
 148   //if(iteration>0) cout<<"iterations "<<iteration<<endl;
 149   return targetParticle;
 150 }
 151 Bool_t AliAnalysisEtSelector::SuspiciousDecayInChain(const UInt_t suspectMotherPdg, const UInt_t suspectDaughterPdg, const TParticle &part, AliStack& stack) const
 152 {
 153   UInt_t partPdg = TMath::Abs(part.GetPdgCode());
 154   if(part.GetFirstMother() == -1)
 155   {
 156     return kFALSE;
 157   }
 158   TParticle *mother = stack.Particle(part.GetFirstMother());
 159   UInt_t motherPdg = TMath::Abs(mother->GetPdgCode());
 160   if((suspectDaughterPdg==partPdg || 2112 == partPdg )&& suspectMotherPdg == motherPdg)
 161   {
 162     return kTRUE;
 163   }
 164   return SuspiciousDecayInChain(suspectMotherPdg, suspectDaughterPdg, *mother, stack);
 165 }
 166
 167 Float_t AliAnalysisEtSelector::ShiftAngle(Float_t phi){//Always returns an angle in radians between -pi<phi<pi
 168   float myphi = phi;
 169   while(myphi>TMath::Pi()){//angle is too high, decrease the angle
 170     myphi = myphi - 2*TMath::Pi();
 171   }
 172   while(myphi<-TMath::Pi()){//angle is too low, increase the angle
 173     myphi = myphi + 2*TMath::Pi();
 174   }
 175   return myphi;
 176 }
 177 Bool_t AliAnalysisEtSelector::PassMinEnergyCut(Double_t e) const
 178 {
 179   return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
 180 }
 181