PWGLF/totEt/AliAnalysisEtSelectorEmcal.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Selection class for EMCAL
   4 //
   5 //*-- Authors: Oystein Djuvsland (Bergen)
   6 //_________________________________________________________________________
   7
   8
   9 #include "AliAnalysisEtSelectorEmcal.h"
  10 #include "AliAnalysisEtCuts.h"
  11 #include "AliEMCALTrack.h"
  12 #include "TParticle.h"
  13 #include <iostream>
  14 #include "AliStack.h"
  15
  16 AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal(AliAnalysisEtCuts* cuts):AliAnalysisEtSelector(cuts)
  17 {
  18
  19 }
  20 AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal():AliAnalysisEtSelector()
  21 {
  22 }
  23
  24 AliAnalysisEtSelectorEmcal::~AliAnalysisEtSelectorEmcal()
  25 {
  26
  27 }
  28
  29 void AliAnalysisEtSelectorEmcal::Init()
  30 {
  31     AliAnalysisEtSelector::Init();
  32 }
  33
  34 Int_t AliAnalysisEtSelectorEmcal::Init(const AliESDEvent* event)
  35 {   // Init
  36
  37     AliAnalysisEtSelector::Init(event);
  38     Printf("Initializing selector for run: %d", event->GetRunNumber());
  39     fInitialized = kTRUE;
  40     return 0;
  41 }
  42
  43 TRefArray* AliAnalysisEtSelectorEmcal::GetClusters()
  44 {   // Get clusters
  45
  46     if(!fClusterArray) fClusterArray = new TRefArray;
  47
  48     if(fClusterArray)
  49     {
  50         fEvent->GetEMCALClusters(fClusterArray);
  51     }
  52     else
  53     {
  54         Printf("Could not initialize cluster array");
  55     }
  56     return fClusterArray;
  57 }
  58
  59 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const AliESDCaloCluster& cl) const
  60 {
  61   Float_t pos[3];
  62   cl.GetPosition(pos);
  63   TVector3 cp(pos);
  64   return TMath::Sin(cp.Theta())*cl.E() > fCuts->GetReconstructedEmcalClusterEnergyCut();
  65 }
  66
  67 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const TParticle& p) const
  68 {
  69   return TMath::Sin(p.Theta())*p.Energy() > fCuts->GetReconstructedEmcalClusterEnergyCut();
  70 }
  71 Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(Double_t e) const
  72 {
  73   return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
  74 }
  75
  76 Bool_t AliAnalysisEtSelectorEmcal::PassDistanceToBadChannelCut(const AliESDCaloCluster& ) const
  77 {
  78     return kTRUE;
  79 }
  80
  81 Bool_t AliAnalysisEtSelectorEmcal::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
  82 {
  83
  84   Int_t nTracksMatched = cluster.GetNTracksMatched();
  85   if(nTracksMatched == 0){
  86     return kTRUE;
  87   }
  88
  89   Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
  90   if(trackMatchedIndex < 0){
  91     return kTRUE;
  92   }
  93   AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
  94   if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random.  It takes ~460 MeV to reach the EMCal
  95
  96   //Float_t recoE = cluster.E();
  97   //Float_t pos[3];
  98
  99 //     //cluster.GetPosition(pos);
 100 //     Int_t trackMatchIdx = cluster.GetTrackMatchedIndex();
 101 //     //Double_t distance = 9999.0;
 102 //     if(trackMatchIdx>-1)
 103 //     {
 104 //       return kTRUE;
 105 //       //distance = CalcTrackClusterDistance(pos, &trackMatchIdx);
 106 //     }
 107
 108
 109     return kFALSE;
 110     //return distance > fCuts->GetEmcalTrackDistanceCut();
 111 }
 112
 113 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const TParticle& part)
 114 {
 115   float myphi = part.Phi();
 116   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 117     return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
 118     && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
 119     && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
 120 }
 121
 122 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliVTrack& part)
 123 {
 124   float myphi = part.Phi();
 125   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 126   return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
 127     && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
 128     && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
 129 }
 130
 131 Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
 132 {
 133   Float_t pos[3];
 134   cluster.GetPosition(pos);
 135   TVector3 cp(pos);
 136   float myphi = cp.Phi();
 137   myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
 138   return TMath::Abs(cp.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
 139     && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
 140     && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
 141 }
 142
 143
 144 Double_t
 145 AliAnalysisEtSelectorEmcal::CalcTrackClusterDistance(const Float_t clsPos[3],Int_t *trkMatchId) const
 146 {   // calculate distance between cluster and closest track
 147
 148     Double_t trkPos[3] = {0,0,0};
 149
 150     Int_t bestTrkMatchId = -1;
 151     Double_t distance = 9999; // init to a big number
 152
 153     Double_t dist = 0;
 154     Double_t distX = 0, distY = 0, distZ = 0;
 155
 156     for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
 157         AliESDtrack *track = static_cast<AliESDtrack*>( fEvent->GetTrack(iTrack) );
 158         if (!track) {
 159             AliError(Form("ERROR: Could not get track %d", iTrack));
 160             continue;
 161         }
 162
 163         // check for approx. eta and phi range before we propagate..
 164         // TBD
 165
 166         AliEMCALTrack emctrack(*track);
 167         if (!emctrack.PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
 168             continue;
 169         }
 170         emctrack.GetXYZ(trkPos);
 171
 172         distX = clsPos[0]-trkPos[0];
 173         distY = clsPos[1]-trkPos[1];
 174         distZ = clsPos[2]-trkPos[2];
 175         dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);
 176
 177         if (dist < distance) {
 178             distance = dist;
 179             bestTrkMatchId = iTrack;
 180         }
 181     } // iTrack
 182
 183     // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
 184     *trkMatchId = bestTrkMatchId;
 185     return distance;
 186 }
 187
 188 UInt_t AliAnalysisEtSelectorEmcal::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
 189   UInt_t correctLabel = TMath::Abs(cluster->GetLabel());
 190   correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
 191
 192   //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
 193 //   if(!stack){return 0;}
 194 //   if( correctLabel > 0 && stack->IsSecondaryFromMaterial(correctLabel)){//if this is flagged as a secondary then we look to see where it really came from
 195 //     TParticle *hitParticle = stack->Particle(correctLabel);
 196 //     if(hitParticle){
 197 //       Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >420);
 198 //       if(partVtxSecondary){//at this point we have something which converted near the detector.  Let's find the mother particle
 199 //      UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
 200 //      if(mothIdx>0){
 201 //        TParticle *mother = stack->Particle(mothIdx);
 202 //        if(mother){
 203 //          if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
 204 //            if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
 205 //              cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
 206 //              cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
 207 //              //so now we know that the particle originated near the cover and within the acceptance of the detector
 208 //              return mothIdx;
 209 //            }
 210 //          }
 211 //        }
 212 //      }
 213 //       }
 214
 215 //     }
 216 //   }
 217   return  correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
 218
 219 }