PWG4/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
  89
  90 //________________________________________________________________________________
  91 AliAnaConvIsolation::~AliAnaConvIsolation() {
  92   //Destructor
  93
  94 }
  95
  96
  97 //________________________________________________________________________
  98 void AliAnaConvIsolation::CreateHistograms()
  99 {
 100         //Create histograms
 101   if(!fHistograms) fHistograms = new TList();
 102   fHistograms->SetName(Form("Isolation_histo_cone_%f_maxPt_%f_sumPt_%f", fConeSize, fSumPtThreshold, fMaxPtThreshold));
 103
 104   fIsoCurve = new TF1("Isolation_curve", fCurveFunction.Data(), 0, 100);
 105   fHistograms->Add(fIsoCurve);
 106
 107   cout << "Creatin isolation histograms conesize :" << fConeSize << endl;
 108
 109   //Create histograms add, to outputlis
 110   for(Int_t i = 0; i < 2; i++){
 111     fhMaxPtInCone[i] = new TH2F(Form("fhMaxPtInCone_%s_%f", (i==0)? "nonIso" : "isolated", fConeSize),
 112                                 Form("Max pt nonIso particle in cone %f vs candidate pt", fConeSize),
 113                                 200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 114     fHistograms->Add(fhMaxPtInCone[i]);
 115
 116     fhSumPtInCone[i] = new TH2F(Form("fhSumPtInCone_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 117                              Form("Sum pt in cone %f vs candidate pt %s", fConeSize,  (i==0)? "nonIsoay" : "isolated"),
 118                              200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 119     fHistograms->Add(fhSumPtInCone[i]);
 120
 121     fhSumPtVsMaxPt[i] = new TH2F(Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 122                                  Form("fhSumPtVsMaxPt_%s_%f",  (i==0)? "nonIso" : "isolated", fConeSize),
 123                                  200, 0, fHistogramMaxPt, 200, 0, fHistogramMaxPt);
 124     fHistograms->Add(fhSumPtVsMaxPt[i]);
 125   }
 126
 127   for(Int_t iIso = 0; iIso < 2; iIso++){
 128     fhPtCandidates[iIso] = new TH1F(Form("fhPtCandidates_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
 129                                     Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
 130                                  200, 0 , fHistogramMaxPt);
 131     fhPtCandidates[iIso]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 132     fhPtCandidates[iIso]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 133     fhPtCandidates[iIso]->SetMarkerStyle(kFullCircle);
 134     fHistograms->Add(fhPtCandidates[iIso]);
 135   }
 136
 137   // for(int iDec = 0; iDec < 2; iDec++) {
 138   //   for(int iIso = 0; iIso < 2; iIso++) {
 139
 140   //     fHistSumPt[iIso][iDec] = new TH1F(Form("fHistSumPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 141   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 142   //                                 150, 0.1, 50);
 143   //     fHistSumPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 144   //     fHistSumPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 145   //     fHistSumPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
 146   //     fHistograms->Add(fHistSumPt[iIso][iDec]);
 147
 148
 149   //     fHistMaxPt[iIso][iDec] = new TH1F(Form("fHistMaxPt_%f_%s_%s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 150   //                                    Form("P_{T} distribution cone %f %s %s", fConeSize, (iIso==0)?"nonIso":"isolated" , (iDec==0)?"noDec":"decay" ),
 151   //                                 150, 0.1, 50);
 152   //     fHistMaxPt[iIso][iDec]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 153   //     fHistMaxPt[iIso][iDec]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 154   //     fHistMaxPt[iIso][iDec]->SetMarkerStyle(kFullCircle);
 155   //     fHistograms->Add(fHistMaxPt[iIso][iDec]);
 156
 157
 158   //   }
 159   // }
 160
 161
 162   for(Int_t iIso = 0; iIso < 2; iIso++){
 163     fhTrackMult[iIso] = new TH1F(Form("fhTrackMult_%s_%f", (iIso==0)? "nonIso" : "isolated", fConeSize),
 164                                     Form("Pt of %s candidates, cone size %f", (iIso==0)? "nonIsolated" : "isolated", fConeSize),
 165                                  150, 0 , 150);
 166     fhTrackMult[iIso]->GetXaxis()->SetTitle("n tracks in event");
 167     fhTrackMult[iIso]->GetYaxis()->SetTitle("dN/dNTracks");
 168     fhTrackMult[iIso]->SetMarkerStyle(kFullCircle);
 169     fHistograms->Add(fhTrackMult[iIso]);
 170   }
 171
 172 }
 173
 174 //________________________________________________________________________
 175 Bool_t AliAnaConvIsolation::IsLeading(AliAODConversionParticle * particle, const TObjArray * tracks, const TObjArray * aodParticles) {
 176   //Check if particle is highest pt particle in cone
 177   for(Int_t ip = 0; ip < aodParticles->GetEntriesFast(); ip++) {
 178         AliAODConversionParticle * oparticle = static_cast<AliAODConversionParticle*>(aodParticles->UncheckedAt(ip));
 179         if(oparticle == particle) continue;
 180         if(oparticle->Pt() > particle->Pt() &&
 181            IsInCone(particle->Eta() - oparticle->Eta(),
 182                                 particle->Phi() - oparticle->Phi(),
 183                                 fConeSize) ) {
 184           return kFALSE;
 185         }
 186   }
 187
 188  for(Int_t it = 0; it < tracks->GetEntriesFast(); it++) {
 189         AliAODTrack * track = static_cast<AliAODTrack*>(tracks->UncheckedAt(it));
 190         if(track->Pt() > particle->Pt()  &&
 191            IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) return kFALSE;
 192   }
 193
 194   return kTRUE;
 195 }
 196
 197 //_________________________________________________________________________
 198 Bool_t AliAnaConvIsolation::IsIsolated(AliAODConversionPhoton * particle, const TClonesArray * const tracks, const Int_t nSpawn, const Int_t * const spawn, Bool_t &leading) {
 199   //See header file for documentation
 200
 201   leading = kTRUE;
 202
 203   Float_t ptSum = 0.0;
 204   Float_t ptMax = 0.0;
 205
 206   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 207     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
 208
 209     if (track) {
 210
 211       if(track->Pt() < fMinPt) continue;
 212
 213                         if(nSpawn && spawn) { ;}
 214
 215       ///Ignore tracks that are grandchildren of pion
 216       // if ( particle->IsMySpawn(track->GetID(), nSpawn, spawn))
 217       //        continue;
 218
 219
 220       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
 221         ptSum += track->Pt();
 222         if(track->Pt() > ptMax) ptMax = track->Pt();
 223         if(track->Pt() > particle->Pt()) leading = kFALSE;
 224       }
 225     } else {
 226       cout << "Bad track"<<endl;
 227     }
 228   }
 229
 230
 231   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
 232
 233   FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
 234
 235   return isolated;
 236
 237 }
 238
 239
 240 //_________________________________________________________________________
 241 Bool_t AliAnaConvIsolation::IsIsolated(const AliAODConversionPhoton * const particle, const TClonesArray * const tracks, Bool_t &leading ) {
 242   //See header file for documentation
 243
 244   leading = kTRUE;
 245
 246   Float_t ptSum = 0.0;
 247   Float_t ptMax = 0.0;
 248
 249   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 250
 251     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
 252
 253     if (track) {
 254
 255       if(track->Pt() < fMinPt) continue;
 256
 257       if ( (track->GetID() == particle->GetTrackLabel(0)) || track->GetID() == particle->GetTrackLabel(1) )
 258         continue;
 259
 260       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), fConeSize) ) {
 261         ptSum += track->Pt();
 262         if(track->Pt() > ptMax) ptMax = track->Pt();
 263         if(track->Pt() > particle->Pt()) leading = kFALSE;
 264       }
 265     } else {
 266       cout << "Bad track"<<endl;
 267     }
 268   }
 269
 270
 271   Bool_t isolated = EvaluateIsolationCriteria( ptSum, particle->Pt());
 272
 273   FillHistograms(particle->Pt(), ptMax, ptSum, isolated, tracks->GetEntriesFast());
 274
 275
 276   return isolated;
 277
 278 }
 279
 280 ///___________________________________________________________________________
 281 void AliAnaConvIsolation::FillHistograms(Float_t pt, Float_t ptMax, Float_t ptSum, Bool_t isolated, Int_t nTracks) {
 282   //Fill histograms
 283   fhMaxPtInCone[isolated]->Fill(pt, ptMax);
 284   fhSumPtInCone[isolated]->Fill(pt, ptSum);
 285   fhSumPtVsMaxPt[isolated]->Fill(ptMax, ptSum);
 286   fhPtCandidates[isolated]->Fill(pt);
 287   fhTrackMult[isolated]->Fill(nTracks);
 288 }
 289
 290
 291 ///_____________________________________________________________________________
 292 Bool_t AliAnaConvIsolation::EvaluateIsolationCriteria(Float_t ptSum, Float_t pt) const {
 293   //Evaluate isolation criteria
 294   return (ptSum < fIsoCurve->Eval(pt));
 295
 296 }