PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx

   1 #include <iostream>
   2 #include "TChain.h"
   3 #include "TTree.h"
   4 #include "TH1F.h"
   5 #include "TCanvas.h"
   6 #include "TString.h"
   7
   8 #include "AliAnalysisTask.h"
   9 #include "AliAnalysisManager.h"
  10 #include "AliAnalysisTaskGammaJet.h"
  11
  12 #include "AliESDEvent.h"
  13 #include "AliESDCaloCluster.h"
  14 #include "AliESDInputHandler.h"
  15
  16 #include "AliAODPWG4ParticleCorrelation.h"
  17 #include "AliAODEvent.h"
  18 #include "AliAODHandler.h"
  19 #include "AliAODCaloCluster.h"
  20 #include "AliGammaConversionAODObject.h"
  21 #include "AliAODJet.h"
  22
  23 // Gamma - jet correlation analysis task
  24 // Authors: Svein Lindal
  25
  26
  27 using namespace std;
  28
  29 ClassImp(AliAnalysisTaskGammaJet)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet() : AliAnalysisTaskSE(),
  33   fOutputList(NULL),
  34   fHistPt(NULL),
  35   fHistPtPhos(NULL),
  36   fHistPtEmcal(NULL),
  37   fHistPtJets(NULL),
  38   fHistGammaJets(NULL),
  39   fHistGammaJetsIso(NULL),
  40   fMinPt(5.0),
  41   fConeSize(0.3),
  42   fPtThreshold(2.0),
  43   fDeltaAODFileName(""),
  44   fPhotons(NULL)
  45 {
  46   // Dummy Constructor
  47 }
  48
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet(const char *name) :
  52   AliAnalysisTaskSE(name),
  53   fOutputList(0),
  54   fHistPt(0),
  55   fHistPtPhos(0),
  56   fHistPtEmcal(0),
  57   fHistPtJets(0),
  58   fHistGammaJets(NULL),
  59   fHistGammaJetsIso(NULL),
  60   fMinPt(0.0),
  61   fConeSize(0.0),
  62   fPtThreshold(0.0),
  63   fDeltaAODFileName(""),
  64   fPhotons(NULL)
  65 {
  66   // Constructor
  67   // Define input and output slots here
  68   DefineInput(0, TChain::Class());
  69   // Output slot #0 id reserved by the base class for AOD
  70
  71   // Output slot #1 writes into a TH1 container
  72   DefineOutput(1, TList::Class());
  73 }
  74
  75 //________________________________________________________________________
  76 void AliAnalysisTaskGammaJet::UserCreateOutputObjects()
  77 {
  78   //Create histograms add, to outputlist
  79   fOutputList = new TList();
  80
  81   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 150, 0.1, 50);
  82   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  83   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  84   fHistPt->SetMarkerStyle(kFullCircle);
  85   fOutputList->Add(fHistPt);
  86
  87   fHistPtPhos = new TH1F("fHistPtPhos", "P_{T} distribution", 150, 0.1, 50);
  88   fHistPtPhos->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  89   fHistPtPhos->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  90   fHistPtPhos->SetMarkerStyle(kFullCircle);
  91   fOutputList->Add(fHistPtPhos);
  92
  93   fHistPtEmcal = new TH1F("fHistPtEmcal", "P_{T} distribution", 150, 0.1, 50);
  94   fHistPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  95   fHistPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  96   fHistPtEmcal->SetMarkerStyle(kFullCircle);
  97   fOutputList->Add(fHistPtEmcal);
  98
  99
 100   fHistPtJets = new TH1F("fHistPtJets", "P_{T} distribution", 150, 0.1, 50);
 101   fHistPtJets->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 102   fHistPtJets->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 103   fHistPtJets->SetMarkerStyle(kFullCircle);
 104   fOutputList->Add(fHistPtJets);
 105
 106   fHistGammaJets = new TH1F("fHistGammaJets", "fHistGammaJets", 200, -TMath::Pi(), 2*TMath::Pi());
 107   fOutputList->Add(fHistGammaJets);
 108
 109   fHistGammaJetsIso = new TH1F("fHistGammaJetsIso", "fHistGammaJetsIso", 200, -TMath::Pi(), 2*TMath::Pi());
 110   fOutputList->Add(fHistGammaJetsIso);
 111
 112   //TNtuple * tuple = new TNtuple("fNtuple", "fNtuple", dPhi,
 113
 114
 115   ///Create AOD branch
 116   fPhotons = new TClonesArray("AliAODPWG4ParticleCorrelation", 0);
 117   fPhotons->SetName("fPhotons");
 118   AddAODBranch("TClonesArray", &fPhotons);
 119
 120   ///Isolation class
 121   // fIsolation = new AliAnaParticleIsolation();
 122   // fIsolation->SetInputAODName("fPhotons");
 123
 124
 125 }
 126
 127 //________________________________________________________________________
 128 void AliAnalysisTaskGammaJet::UserExec(Option_t *)
 129 {
 130   // Main loop
 131   // Called for each event
 132
 133
 134   //Clear stuff for new event
 135   CleanUp();
 136
 137
 138   ///Get AOD event
 139   AliAODEvent * aodEvent = GetAODEvent();
 140   if(!aodEvent) {
 141     AliError("No AOD event!!");
 142     return;
 143   }
 144
 145   //FillPWG4PartCorrBranch(convGamma, fPhotons, "ConvGamma");
 146   //fIsolation->MakeAnalysisFillAOD();
 147
 148   ProcessConvGamma(aodEvent);
 149   ProcessCalorimeters(aodEvent);
 150
 151
 152   PostData(1, fOutputList);
 153
 154 }
 155 //_____________________________________________________________________
 156 void AliAnalysisTaskGammaJet::Terminate(Option_t *) {
 157   // Draw result to the screen
 158   // Called once at the end of the query
 159 }
 160
 161 //_____________________________________________________________________
 162 AliAODEvent * AliAnalysisTaskGammaJet::GetAODEvent() {
 163   //Get the AOD event from whereever it might be
 164   AliAODEvent * aodEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 165   if(!aodEvent) {
 166     aodEvent = AODEvent();
 167   }
 168
 169   return aodEvent;
 170
 171 }
 172
 173
 174 //_____________________________________________________________________
 175 TClonesArray * AliAnalysisTaskGammaJet::GetConversionGammas(const AliAODEvent * aodEvent) {
 176
 177   //Get Conversion gamma branch of AOD. First try standard AOD
 178   TClonesArray * convGamma = dynamic_cast<TClonesArray*>(aodEvent->FindListObject("GammaConv_gamma"));
 179
 180   //If it's there, send it back
 181   if(convGamma)  return convGamma;
 182
 183   //If AOD not in standard file have to locate it in delta AOD
 184   if( !(fDeltaAODFileName.Length() > 0)  ) return NULL;
 185
 186   AliAODHandler * aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 187   if(aodHandler) {
 188     AliAODExtension * gExt = dynamic_cast<AliAODExtension*>(aodHandler->GetExtensions()->FindObject(fDeltaAODFileName));
 189     if(gExt) {
 190       AliAODEvent * gcEvent = gExt->GetAOD();
 191       return dynamic_cast<TClonesArray*>(gcEvent->FindListObject("GammaConv_gamma"));
 192     }
 193   }
 194   return NULL;
 195 }
 196
 197
 198 //_____________________________________________________________________
 199 void AliAnalysisTaskGammaJet::FillPWG4PartCorrBranch( TClonesArray * gcBranch, TClonesArray * partCorrBranch , TString detector ) {
 200
 201   for(int i = 0; i < gcBranch->GetEntriesFast(); i++) {
 202     AliGammaConversionAODObject * gcObject = dynamic_cast<AliGammaConversionAODObject*>(gcBranch->At(i));
 203     if ( gcObject ) {
 204       AliAODPWG4ParticleCorrelation pc(gcObject->Px(), gcObject->Py(), gcObject->Pz(), gcObject->E());
 205       pc.SetTagged(gcObject->IsTagged());
 206       pc.SetTrackLabel(gcObject->GetLabel1(), gcObject->GetLabel2());
 207       pc.SetDetector(detector);
 208       new((*partCorrBranch)[i]) AliAODPWG4ParticleCorrelation(pc);
 209
 210     } else {
 211       AliError(Form("Couldn't get gamma conversion aod object"));
 212     }
 213
 214   }
 215 }
 216
 217
 218 //_____________________________________________________________________
 219 AliAODPWG4ParticleCorrelation * AliAnalysisTaskGammaJet::PWG4PartFromGammaConvAODObject(AliGammaConversionAODObject * gcObject, TString detector ) {
 220
 221   AliAODPWG4ParticleCorrelation * pc = new AliAODPWG4ParticleCorrelation(gcObject->Px(), gcObject->Py(), gcObject->Pz(), gcObject->E());
 222   pc->SetTagged(gcObject->IsTagged());
 223   pc->SetTrackLabel(gcObject->GetLabel1(), gcObject->GetLabel2());
 224   pc->SetDetector(detector);
 225   return pc;
 226 }
 227
 228
 229 //_________________________________________________________________________
 230 void AliAnalysisTaskGammaJet::CleanUp() {
 231   fPhotons->Delete();
 232 }
 233
 234 //_________________________________________________________________________
 235 Bool_t AliAnalysisTaskGammaJet::IsIsolated( AliAODPWG4ParticleCorrelation * particle, TClonesArray * tracks, Float_t coneSize, Float_t ptThreshold ) {
 236   //See header file for documentation
 237   for(int it = 0; it < tracks->GetEntriesFast(); it++) {
 238     if ( (it == particle->GetTrackLabel(0)) || it == particle->GetTrackLabel(1) )
 239       continue;
 240
 241     //BALLE Svein:How are you checking the calorimeters for whether they are decay particles ?
 242
 243     AliAODTrack * track = dynamic_cast<AliAODTrack*>(tracks->At(it));
 244     if (track) {
 245       if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), coneSize) ) {
 246         if (track->Pt() > ptThreshold) {
 247           return kFALSE;
 248         }
 249       }
 250     } else {
 251       AliError(Form("Bad track!!!! "));
 252     }
 253   }
 254
 255   //No particle above threshold, it's isolated
 256   return kTRUE;
 257 }
 258
 259
 260 //______________________________________________________________________________________________
 261 void AliAnalysisTaskGammaJet::ProcessCalorimeters( const AliAODEvent * const aodEvent ) {
 262
 263   TClonesArray * clusters = aodEvent->GetCaloClusters();
 264
 265
 266   for(int ic = 0; ic < clusters->GetEntriesFast(); ic++) {
 267     AliAODCaloCluster * cluster = dynamic_cast<AliAODCaloCluster*>(clusters->At(ic));
 268     if (!cluster) {
 269       AliError(Form("Error getting cluster"));
 270       continue;
 271     }
 272
 273
 274     if (cluster->GetNCells() < 6) continue;
 275     if (cluster->GetEmcCpvDistance() < 15) continue;
 276
 277     TLorentzVector tlvec;
 278
 279     AliAODVertex * vertex = aodEvent->GetPrimaryVertex();
 280     Double_t vertexPosition[3];
 281     vertex->GetXYZ(vertexPosition);
 282     cluster->GetMomentum(tlvec, vertexPosition);
 283     if (tlvec.Pt() < GetMinPt()) continue;
 284
 285     AliAODPWG4ParticleCorrelation * photon = new AliAODPWG4ParticleCorrelation(tlvec);
 286
 287     photon->SetIsolated( IsIsolated(photon, aodEvent->GetTracks(), GetConeSize(), GetPtThreshold()) );
 288     CorrelateWithJets(photon, aodEvent->GetJets());
 289   }
 290
 291 }
 292 //___________________________________________________________________________________________
 293 void AliAnalysisTaskGammaJet::ProcessConvGamma( const AliAODEvent * const aodEvent ) {
 294
 295   TClonesArray * convGamma = GetConversionGammas(aodEvent);
 296   if(!convGamma) {
 297     AliError(Form("No convgamma"));
 298     return;
 299   }
 300
 301   for (Int_t iPhot = 0; iPhot < convGamma->GetEntriesFast(); iPhot++) {
 302     AliGammaConversionAODObject * aodO = dynamic_cast<AliGammaConversionAODObject*>(convGamma->At(iPhot));
 303
 304     if (!aodO) {
 305       AliError(Form("ERROR: Could not receive ga %d\n", iPhot));
 306       continue;
 307     }
 308
 309     if(aodO->Pt() < GetMinPt()) continue;
 310
 311
 312     //Use the AODPWG4PartCorr shit!
 313     AliAODPWG4ParticleCorrelation * photon = PWG4PartFromGammaConvAODObject(aodO, "ConvGamma");
 314     photon->SetIsolated( IsIsolated(photon, aodEvent->GetTracks(), GetConeSize(), GetPtThreshold()) );
 315
 316
 317     // if ( (aodO->Phi()) < 0 )
 318     //   cout << aodO->Phi() << endl;
 319
 320     CorrelateWithJets(photon, aodEvent->GetJets());
 321
 322     fHistPt->Fill(photon->Pt());
 323     delete photon;
 324
 325   }
 326 }
 327
 328 void AliAnalysisTaskGammaJet::CorrelateWithJets(AliAODPWG4ParticleCorrelation * photon, const TClonesArray * const jets) {
 329   //See header file for documentation
 330   if (jets) {
 331     for(int ij = 0; ij < jets->GetEntriesFast(); ij++) {
 332       AliAODJet * jet = dynamic_cast<AliAODJet*>(jets->At(ij));
 333       if(jet) {
 334         fHistPtJets->Fill(jet->Pt());
 335
 336         Float_t dPhi = TMath::Abs(photon->Phi() - jet->Phi());
 337         if (photon->IsIsolated())
 338           fHistGammaJetsIso->Fill(dPhi, jet->Pt()/photon->Pt());
 339         else
 340           fHistGammaJets->Fill(dPhi);
 341
 342       }
 343     }
 344   }
 345 }