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