PWGGA/GammaConv/ConvCorrelations/AliAnaConvIsolation.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * ALICE Experiment at CERN, All rights reserved.                         *
   4  *                                                                        *
   5  * Primary Author: Svein Lindal <slindal@fys.uio.no>                      *
   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 /// @file   AliAnalysisTaskGammaJet.cxx
  17 /// @author Svein Lindal
  18 /// @brief  Class used to run isolation studies of conversion gamma/pions
  19
  20
  21 #include "AliAnaConvIsolation.h"
  22 #include "AliAODTrack.h"
  23
  24 #include "TObject.h"
  25 #include "TClonesArray.h"
  26 #include "TH2F.h"
  27 #include "TList.h"
  28 #include "AliAODConversionPhoton.h"
  29
  30
  31 #include <iostream>
  32
  33 // Gamma - jet correlation analysis task
  34 // Author: Svein Lindal
  35
  36
  37 using namespace std;
  38 ClassImp(AliAnaConvIsolation)
  39
  40 //________________________________________________________________________
  41 AliAnaConvIsolation::AliAnaConvIsolation () : TObject(),
  42   fIsoCurve(NULL),
  43   fCurveFunction("0.1*x"),
  44   fConeSize(0),
  45   fMinPt(0.1),
  46   fMaxPtThreshold(0),
  47   fSumPtThreshold(0),
  48   fMaxPtFraction(0),
  49   fSumPtFraction(0)
  50   // fHistograms(NULL),
  51   // fHistogramMaxPt(50)
  52 {
  53   //Constructor
  54   // for(Int_t i = 0; i < 2; i++){
  55   //   fhMaxPtInCone[i] = NULL;
  56   //   fhSumPtInCone[i] = NULL;
  57   //   fhSumPtVsMaxPt[i] = NULL;
  58   //   fhPtCandidates[i] = NULL;
  59   //   fhTrackMult[i] = NULL;
  60   // }
  61
  62
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnaConvIsolation::AliAnaConvIsolation(Float_t coneSize, Float_t maxPtThreshold, Float_t sumPtThreshold, Float_t maxPtFraction, Float_t sumPtFraction) :
  67   TObject(),
  68   fIsoCurve(NULL),
  69   fCurveFunction("0.1*x"),
  70   fConeSize(coneSize),
  71   fMinPt(0.1),
  72   fMaxPtThreshold(maxPtThreshold),
  73   fSumPtThreshold(sumPtThreshold),
  74   fMaxPtFraction(maxPtFraction),
  75   fSumPtFraction(sumPtFraction)
  76   // fHistograms(NULL),
  77   // fHistogramMaxPt(50)
  78 {
  79   //Constructor
  80   // for(Int_t i = 0; i < 2; i++){
  81   //   fhMaxPtInCone[i] = NULL;
  82   //   fhSumPtInCone[i] = NULL;
  83   //   fhSumPtVsMaxPt[i] = NULL;
  84   //   fhPtCandidates[i] = NULL;
  85   //   fhTrackMult[i] = NULL;
  86   // }
  87
  88   fIsoCurve = new TF1("Isolation_curve", fCurveFunction.Data(), 0, 100);
  89
  90 }
  91
  92
  93 //________________________________________________________________________________
  94 AliAnaConvIsolation::~AliAnaConvIsolation() {
  95   //Destructor
  96
  97 }
  98
  99
 100 //________________________________________________________________________
 101 // void AliAnaConvIsolation::CreateHistograms()
 102 // {
 103         //Create histograms
 104   // if(!fHistograms) fHistograms = new TList();
 105   // fHistograms->SetName(Form("Isolation_histo_cone_%f_maxPt_%f_sumPt_%f", fConeSize, fSumPtThreshold, fMaxPtThreshold));
 106
 107   //fHistograms->Add(fIsoCurve);
 108
 109   //cout << "Creatin isolation histograms conesize :" << fConeSize << endl;
 110
 111   //Create histograms add, to outputlis
 112   // for(Int_t i = 0; i < 2; i++){
 113   //   fhMaxPtInCone[i] = new TH2F(Form("fhMaxPtInCone_%s_%f", (i==0)? "nonIso" : "isolated", fConeSize),
 114   //                            Form("Max pt nonIso particle in cone %f vs candidate pt", fConeSize),
 115   //                            200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 116   //   fHistograms->Add(fhMaxPtInCone[i]);
 117
 118   //   fhSumPtInCone[i] = new TH2F(Form("fhSumPtInCone_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 119   //                         Form("Sum pt in cone %f vs candidate pt %s", fConeSize,  (i==0)? "nonIsoay" : "isolated"),
 120   //                         200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 121   //   fHistograms->Add(fhSumPtInCone[i]);
 122
 123   //   fhSumPtVsMaxPt[i] = new TH2F(Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 124   //                             Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 125   //                             200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 126   //   fHistograms->Add(fhSumPtVsMaxPt[i]);
 127   // }
 128
 129   // for(Int_t iIso = 0; iIso < 2; iIso++){
 130   //   fhPtCandidates[iIso] = new TH1F(Form("fhPtCandidates_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
 131   //                                Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
 132   //                             200, 0 , fHistogramMaxPt);
 133   //   fhPtCandidates[iIso]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 134   //   fhPtCandidates[iIso]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 135   //   fhPtCandidates[iIso]->SetMarkerStyle(kFullCircle);
 136   //   fHistograms->Add(fhPtCandidates[iIso]);
 137   // }
 138
 139   // for(int iDec = 0; iDec < 2; iDec++) {
 140   //   for(int iIso = 0; iIso < 2; iIso++) {
 141
 142   //     fHistSumPt[iIso][iDec] = new TH1F(Form("fHistSumPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 143   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 144   //                                 150, 0.1, 50);
 145   //     fHistSumPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 146   //     fHistSumPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 147   //     fHistSumPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
 148   //     fHistograms->Add(fHistSumPt[iIso][iDec]);
 149
 150
 151   //     fHistMaxPt[iIso][iDec] = new TH1F(Form("fHistMaxPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 152   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 153   //                                 150, 0.1, 50);
 154   //     fHistMaxPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 155   //     fHistMaxPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 156   //     fHistMaxPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
 157   //     fHistograms->Add(fHistMaxPt[iIso][iDec]);
 158
 159
 160   //   }
 161   // }
 162
 163
 164   // for(Int_t iIso = 0; iIso < 2; iIso++){
 165   //   fhTrackMult[iIso] = new TH1F(Form("fhTrackMult_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
 166   //                                Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
 167   //                             150, 0 , 150);
 168   //   fhTrackMult[iIso]->GetXaxis()->SetTitle("n tracks in event");
 169   //   fhTrackMult[iIso]->GetYaxis()->SetTitle("dN/dNTracks");
 170   //   fhTrackMult[iIso]->SetMarkerStyle(kFullCircle);
 171   //   fHistograms->Add(fhTrackMult[iIso]);
 172   // }
 173
 174 //}
 175
 176 //________________________________________________________________________
 177 Bool_t AliAnaConvIsolation::IsLeading(AliAODConversionParticle * particle, const TObjArray * tracks, const TObjArray * aodParticles) {
 178   //Check if particle is highest pt particle in cone
 179   for(Int_t ip = 0; ip < aodParticles->GetEntriesFast(); ip++) {
 180         AliAODConversionParticle * oparticle = static_cast<AliAODConversionParticle*>(aodParticles->UncheckedAt(ip));
 181         if(oparticle == particle) continue;
 182         if(oparticle->Pt() > particle->Pt() &&
 183            IsInCone(particle->Eta() - oparticle->Eta(),
 184                                 particle->Phi() - oparticle->Phi(),
 185                                 fConeSize) ) {
 186           return kFALSE;
 187         }
 188   }
 189
 190  for(Int_t it = 0; it < tracks->GetEntriesFast(); it++) {
 191         AliAODTrack * track = static_cast<AliAODTrack*>(tracks->UncheckedAt(it));
 192         if(track->Pt() > particle->Pt()  &&
 193            IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) return kFALSE;
 194   }
 195
 196   return kTRUE;
 197 }
 198
 199 ///________________________________________________________________________
 200 Bool_t AliAnaConvIsolation::IsLeading(const AliAODConversionParticle * particle, const TObjArray * tracks, const Int_t tIDs[4]) const {
 201   //Is there a higher pt particle within cone ?
 202   Float_t tPt = particle->Pt();
 203   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 204     AliVTrack * track = static_cast<AliVTrack*>(tracks->UncheckedAt(it));
 205         Int_t tid = track->GetID();
 206         if(tid == tIDs[0] || tid == tIDs[1] || tid == tIDs[2] || tid == tIDs[3]) continue;
 207         if(track->Pt() > tPt && IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), 0.8)) return kFALSE;
 208   }
 209   return kTRUE;
 210 }
 211
 212 //_________________________________________________________________________
 213 Bool_t AliAnaConvIsolation::IsIsolated(AliAODConversionPhoton * particle, const TClonesArray * const tracks, const Int_t nSpawn, const Int_t * const spawn, Bool_t &leading) {
 214   //See header file for documentation
 215
 216   leading = kTRUE;
 217
 218   Float_t ptSum = 0.0;
 219   Float_t ptMax = 0.0;
 220
 221   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 222     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
 223
 224     if (track) {
 225
 226       if(track->Pt() < fMinPt) continue;
 227
 228                         if(nSpawn && spawn) { ;}
 229
 230       ///Ignore tracks that are grandchildren of pion
 231       // if ( particle->IsMySpawn(track->GetID(), nSpawn, spawn))
 232       //        continue;
 233
 234
 235       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
 236         ptSum += track->Pt();
 237         if(track->Pt() > ptMax) ptMax = track->Pt();
 238         if(track->Pt() > particle->Pt()) leading = kFALSE;
 239       }
 240     } else {
 241       cout << "Bad track"<<endl;
 242     }
 243   }
 244
 245
 246   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
 247
 248   //FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
 249
 250   return isolated;
 251
 252 }
 253
 254
 255 //_________________________________________________________________________
 256 Bool_t AliAnaConvIsolation::IsIsolated(const AliAODConversionPhoton * const particle, const TClonesArray * const tracks, Bool_t &leading ) {
 257   //See header file for documentation
 258
 259   leading = kTRUE;
 260
 261   Float_t ptSum = 0.0;
 262   Float_t ptMax = 0.0;
 263
 264   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 265
 266     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
 267
 268     if (track) {
 269
 270       if(track->Pt() < fMinPt) continue;
 271
 272       if ( (track->GetID() == particle->GetTrackLabel(0)) || track->GetID() == particle->GetTrackLabel(1) )
 273         continue;
 274
 275       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
 276         ptSum += track->Pt();
 277         if(track->Pt() > ptMax) ptMax = track->Pt();
 278         if(track->Pt() > particle->Pt()) leading = kFALSE;
 279       }
 280     } else {
 281       cout << "Bad track"<<endl;
 282     }
 283   }
 284
 285
 286   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
 287
 288   //FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
 289
 290
 291   return isolated;
 292
 293 }
 294
 295 // ///___________________________________________________________________________
 296 // void AliAnaConvIsolation::FillHistograms(Float_t pt, Float_t ptMax, Float_t ptSum, Bool_t isolated, Int_t nTracks) {
 297 //   //Fill histograms
 298 //   fhMaxPtInCone[isolated]->Fill(pt, ptMax);
 299 //   fhSumPtInCone[isolated]->Fill(pt, ptSum);
 300 //   fhSumPtVsMaxPt[isolated]->Fill(ptMax, ptSum);
 301 //   fhPtCandidates[isolated]->Fill(pt);
 302 //   fhTrackMult[isolated]->Fill(nTracks);
 303 // }
 304
 305
 306 ///_____________________________________________________________________________
 307 Bool_t AliAnaConvIsolation::EvaluateIsolationCriteria(Float_t ptSum, Float_t pt) const {
 308   //Evaluate isolation criteria
 309   return (ptSum < fIsoCurve->Eval(pt));
 310
 311 }