PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

   1 // $Id$
   2
   3 #include "AliAnalysisTaskEMCALIsoPhoton.h"
   4
   5 #include <TCanvas.h>
   6 #include <TChain.h>
   7 #include <TFile.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TLorentzVector.h>
  13 #include <TList.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTask.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliEMCALRecoUtils.h"
  19 #include "AliESDCaloCells.h"
  20 #include "AliESDCaloCluster.h"
  21 #include "AliESDEvent.h"
  22 #include "AliESDHeader.h"
  23 #include "AliESDInputHandler.h"
  24 #include "AliESDUtils.h"
  25 #include "AliESDtrack.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliMCEvent.h"
  28 #include "AliMCEventHandler.h"
  29 #include "AliStack.h"
  30 #include "AliV0vertexer.h"
  31 #include "AliVCluster.h"
  32
  33 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  37   AliAnalysisTaskSE(),
  38   fCaloClusters(0),
  39   fSelPrimTracks(0),
  40   fTracks(0),
  41   fEMCalCells(0),
  42   fPrTrCuts(0),
  43   fGeom(0x0),
  44   fGeoName("EMCAL_COMPLETEV1"),
  45   fPeriod("LHC11c"),
  46   fTrigBit("kEMC7"),
  47   fIsTrain(0),
  48   fIsMc(0),
  49   fDebug(0),
  50   fPathStrOpt("/"),
  51   fExoticCut(0.97),
  52   fIsoConeR(0.4),
  53   fNDimensions(7),
  54   fECut(3.),
  55   fTrackMult(0),
  56   fMcIdFamily(""),
  57   fNClusForDirPho(0),
  58   fDirPhoPt(0),
  59   fHigherPtCone(0),
  60   fESD(0),
  61   fMCEvent(0),
  62   fStack(0),
  63   fOutputList(0),
  64   fEvtSel(0),
  65   fNClusEt10(0),
  66   fRecoPV(0),
  67   fPVtxZ(0),
  68   fTrMultDist(0),
  69   fMCDirPhotonPtEtaPhi(0),
  70   fDecayPhotonPtMC(0),
  71   fCellAbsIdVsAmpl(0),
  72   fNClusHighClusE(0),
  73   fHigherPtConeM02(0),
  74   fClusEtMcPt(0),
  75   fClusMcDetaDphi(0),
  76   fNClusPerPho(0),
  77   fMCDirPhotonPtEtaPhiNoClus(0),
  78   fHnOutput(0)
  79 {
  80   // Default constructor.
  81 }
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
  85   AliAnalysisTaskSE(name),
  86   fCaloClusters(0),
  87   fSelPrimTracks(0),
  88   fTracks(0),
  89   fEMCalCells(0),
  90   fPrTrCuts(0),
  91   fGeom(0x0),
  92   fGeoName("EMCAL_COMPLETEV1"),
  93   fPeriod("LHC11c"),
  94   fTrigBit("kEMC7"),
  95   fIsTrain(0),
  96   fIsMc(0),
  97   fDebug(0),
  98   fPathStrOpt("/"),
  99   fExoticCut(0.97),
 100   fIsoConeR(0.4),
 101   fNDimensions(7),
 102   fECut(3.),
 103   fTrackMult(0),
 104   fMcIdFamily(""),
 105   fNClusForDirPho(0),
 106   fDirPhoPt(0),
 107   fHigherPtCone(0),
 108   fESD(0),
 109   fMCEvent(0),
 110   fStack(0),
 111   fOutputList(0),
 112   fEvtSel(0),
 113   fNClusEt10(0),
 114   fRecoPV(0),
 115   fPVtxZ(0),
 116   fTrMultDist(0),
 117   fMCDirPhotonPtEtaPhi(0),
 118   fDecayPhotonPtMC(0),
 119   fCellAbsIdVsAmpl(0),
 120   fNClusHighClusE(0),
 121   fHigherPtConeM02(0),
 122   fClusEtMcPt(0),
 123   fClusMcDetaDphi(0),
 124   fNClusPerPho(0),
 125   fMCDirPhotonPtEtaPhiNoClus(0),
 126   fHnOutput(0)
 127 {
 128   // Constructor
 129
 130   // Define input and output slots here
 131   // Input slot #0 works with a TChain
 132   DefineInput(0, TChain::Class());
 133   // Output slot #0 id reserved by the base class for AOD
 134   // Output slot #1 writes into a TH1 container
 135   DefineOutput(1, TList::Class());
 136 }
 137
 138 //________________________________________________________________________
 139 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 140 {
 141   // Create histograms, called once.
 142
 143   fCaloClusters = new TRefArray();
 144   fSelPrimTracks = new TObjArray();
 145
 146
 147   fOutputList = new TList();
 148   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 149
 150   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 151
 152   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 153   fOutputList->Add(fEvtSel);
 154
 155   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 156   fOutputList->Add(fNClusEt10);
 157
 158   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 159   fOutputList->Add(fRecoPV);
 160
 161   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 162   fOutputList->Add(fPVtxZ);
 163
 164   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
 165   fOutputList->Add(fTrMultDist);
 166
 167   fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
 168   fMCDirPhotonPtEtaPhi->Sumw2();
 169   fOutputList->Add(fMCDirPhotonPtEtaPhi);
 170
 171   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
 172   fDecayPhotonPtMC->Sumw2();
 173   fOutputList->Add(fDecayPhotonPtMC);
 174
 175   fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
 176   fOutputList->Add(fCellAbsIdVsAmpl);
 177
 178   fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
 179   fOutputList->Add(fNClusHighClusE);
 180
 181   fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",1000,0,100,400,0,4);
 182   fOutputList->Add(fHigherPtConeM02);
 183
 184   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
 185   fOutputList->Add(fClusEtMcPt);
 186
 187   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
 188   fOutputList->Add(fClusMcDetaDphi);
 189
 190   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
 191   fOutputList->Add(fNClusPerPho);
 192
 193   fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and  #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
 194   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
 195
 196   Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000,  nCeIsoNoUE=1000,  nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=101;
 197   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt};
 198   fNDimensions = sizeof(bins)/sizeof(Int_t);
 199   const Int_t ndims =   fNDimensions;
 200   Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-1.5};
 201   Double_t xmax[] = { 100., 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5};
 202   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax);
 203   fOutputList->Add(fHnOutput);
 204
 205
 206
 207   PostData(1, fOutputList);
 208 }
 209
 210 //________________________________________________________________________
 211 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 212 {
 213   // Main loop, called for each event.
 214
 215   // event trigger selection
 216   Bool_t isSelected = 0;
 217   if(fPeriod.Contains("11a")){
 218     if(fTrigBit.Contains("kEMC"))
 219       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 220     else
 221       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 222   }
 223   else{
 224     if(fTrigBit.Contains("kEMC"))
 225       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 226     else
 227       if(fTrigBit.Contains("kMB"))
 228         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 229       else
 230         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 231   }
 232   if(fIsMc)
 233     isSelected = kTRUE;
 234   if(fDebug)
 235     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 236   if(!isSelected )
 237         return;
 238   if(fIsMc){
 239     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 240     TFile *file = (TFile*)tree->GetCurrentFile();
 241     TString filename = file->GetName();
 242     if(!filename.Contains(fPathStrOpt.Data()))
 243       return;
 244   }
 245   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 246   if (!fESD) {
 247     printf("ERROR: fESD not available\n");
 248     return;
 249   }
 250
 251   fEvtSel->Fill(0);
 252   if(fDebug)
 253     printf("fESD is ok\n");
 254
 255   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 256   if(!pv) {
 257     fRecoPV->Fill(0);
 258     return;
 259   }
 260   else
 261     fRecoPV->Fill(1);
 262   fPVtxZ->Fill(pv->GetZ());
 263   if(TMath::Abs(pv->GetZ())>15)
 264     return;
 265   if(fDebug)
 266     printf("passed vertex cut\n");
 267
 268   fEvtSel->Fill(1);
 269
 270   fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 271   if(!fTracks){
 272     AliError("Track array in event is NULL!");
 273     return;
 274   }
 275   // Track loop to fill a pT spectrum
 276   const Int_t Ntracks = fTracks->GetEntriesFast();
 277   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 278     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 279     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 280     AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks));
 281     if (!track)
 282       continue;
 283     if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 284       fSelPrimTracks->Add(track);
 285       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 286     }
 287   }
 288
 289   if(!fIsTrain){
 290     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 291       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 292         break;
 293       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 294     }
 295   }
 296   AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 297   fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 298   fTrMultDist->Fill(fTrackMult);
 299
 300   fESD->GetEMCALClusters(fCaloClusters);
 301   fEMCalCells = fESD->GetEMCALCells();
 302
 303
 304
 305   fMCEvent = MCEvent();
 306   if(fMCEvent){
 307     if(fDebug)
 308       std::cout<<"MCevent exists\n";
 309     fStack = (AliStack*)fMCEvent->Stack();
 310   }
 311   else{
 312     if(fDebug)
 313       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 314   }
 315   FollowGamma();
 316   if(fDebug)
 317     printf("passed calling of FollowGamma\n");
 318   FillClusHists();
 319   if(fDebug)
 320     printf("passed calling of FillClusHists\n");
 321   FillMcHists();
 322   if(fDebug)
 323     printf("passed calling of FillMcHists\n");
 324
 325   fCaloClusters->Clear();
 326   fSelPrimTracks->Clear();
 327   fNClusForDirPho = 0;
 328
 329   PostData(1, fOutputList);
 330 }
 331
 332 //________________________________________________________________________
 333 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 334 {
 335   // Fill cluster histograms.
 336
 337   if(!fCaloClusters)
 338     return;
 339   const Int_t nclus = fCaloClusters->GetEntries();
 340   if(nclus==0)
 341     return;
 342   if(fDebug)
 343     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 344   Double_t maxE = 0;
 345   Int_t nclus10 = 0;
 346   Double_t ptmc=-1;
 347   for(Int_t ic=0;ic<nclus;ic++){
 348     maxE=0;
 349     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 350     if(!c)
 351       continue;
 352     if(!c->IsEMCAL())
 353       continue;
 354     if(c->E()<fECut)
 355       continue;
 356     Short_t id;
 357     Double_t Emax = GetMaxCellEnergy( c, id);
 358     Double_t Ecross = GetCrossEnergy( c, id);
 359     if((1-Ecross/Emax)>fExoticCut)
 360       continue;
 361     Float_t clsPos[3] = {0,0,0};
 362     c->GetPosition(clsPos);
 363     TVector3 clsVec(clsPos);
 364     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 365     if(fDebug)
 366       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 367     if(Et>10)
 368       nclus10++;
 369     Float_t ceiso, cephiband, cecore;
 370     Float_t triso, trphiband, trcore;
 371     Float_t alliso, allphiband, allcore;
 372     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 373     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 374     GetCeIso(clsVec, ceiso, cephiband, cecore);
 375     GetTrIso(clsVec, triso, trphiband, trcore);
 376     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 377     if(Et>10 && Et<15 && dr>0.025){
 378       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 379       if(fDebug)
 380         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 381     }
 382     alliso = ceiso + triso;
 383     allphiband = cephiband + trphiband;
 384     allcore = cecore + trcore;
 385     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 386     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 387     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 388     const Int_t ndims =   fNDimensions;
 389     Double_t outputValues[ndims];
 390     ptmc = GetClusSource(c);
 391     outputValues[0] = Et;
 392     outputValues[1] = c->GetM02();
 393     outputValues[2] = ceiso-cecore-ceisoue;
 394     outputValues[3] = triso-trisoue;
 395     outputValues[4] = alliso-cecore-allisoue;
 396     outputValues[5] = ceiso-Et;
 397     outputValues[6] = alliso-Et;
 398     outputValues[7] = c->GetTrackDx();
 399     outputValues[8] = c->GetTrackDz();
 400     outputValues[9] = clsVec.Eta();
 401     outputValues[10] = clsVec.Phi();
 402     outputValues[11] = fEMCalCells->GetCellTime(id);
 403     outputValues[12] = fTrackMult;
 404     outputValues[13] = ptmc;
 405     fHnOutput->Fill(outputValues);
 406     if(c->E()>maxE)
 407       maxE = c->E();
 408   }
 409   fNClusHighClusE->Fill(maxE,nclus);
 410   fNClusEt10->Fill(nclus10);
 411   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 412 }
 413
 414 //________________________________________________________________________
 415 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 416 {
 417   // Get cell isolation.
 418
 419   if(!fEMCalCells)
 420     return;
 421   const Int_t ncells = fEMCalCells->GetNumberOfCells();
 422   Float_t totiso=0;
 423   Float_t totband=0;
 424   Float_t totcore=0;
 425   Float_t etacl = vec.Eta();
 426   Float_t phicl = vec.Phi();
 427   Float_t thetacl = vec.Theta();
 428   if(phicl<0)
 429     phicl+=TMath::TwoPi();
 430   Int_t absid = -1;
 431   Float_t eta=-1, phi=-1;
 432   for(int icell=0;icell<ncells;icell++){
 433     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 434     Float_t celltime = fEMCalCells->GetCellTime(absid);
 435     if(celltime>2e-8 || celltime<-2e-8)
 436       continue;
 437     if(!fGeom)
 438       return;
 439     fGeom->EtaPhiFromIndex(absid,eta,phi);
 440     Float_t dphi = TMath::Abs(phi-phicl);
 441     Float_t deta = TMath::Abs(eta-etacl);
 442     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 443     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 444     if(R<fIsoConeR){
 445       totiso += etcell;
 446       if(R<0.04)
 447         totcore += etcell;
 448     }
 449     else{
 450       if(dphi>fIsoConeR)
 451         continue;
 452       if(deta<fIsoConeR)
 453         continue;
 454       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 455     }
 456   }
 457   iso = totiso;
 458   phiband = totband;
 459   core = totcore;
 460 }
 461
 462 //________________________________________________________________________
 463 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 464 {
 465   // Get track isolation.
 466
 467   if(!fSelPrimTracks)
 468     return;
 469   fHigherPtCone = 0;
 470   const Int_t ntracks = fSelPrimTracks->GetEntries();
 471   Double_t totiso=0;
 472   Double_t totband=0;
 473   Double_t totcore=0;
 474   Double_t etacl = vec.Eta();
 475   Double_t phicl = vec.Phi();
 476   if(phicl<0)
 477     phicl+=TMath::TwoPi();
 478   for(int itrack=0;itrack<ntracks;itrack++){
 479     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 480     if(!track)
 481       continue;
 482     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 483     Double_t deta = TMath::Abs(track->Eta()-etacl);
 484     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 485     Double_t pt = track->Pt();
 486     if(pt>fHigherPtCone)
 487       fHigherPtCone = pt;
 488     if(R<fIsoConeR){
 489       totiso += track->Pt();
 490       if(R<0.04)
 491         totcore += pt;
 492     }
 493     else{
 494       if(dphi>fIsoConeR)
 495         continue;
 496       if(deta<fIsoConeR)
 497         continue;
 498       totband += track->Pt();
 499     }
 500   }
 501   iso = totiso;
 502   phiband = totband;
 503   core = totcore;
 504 }
 505
 506 //________________________________________________________________________
 507 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 508 {
 509   // Calculate the energy of cross cells around the leading cell.
 510
 511   AliVCaloCells *cells = 0;
 512   cells = fESD->GetEMCALCells();
 513   if (!cells)
 514     return 0;
 515
 516   if (!fGeom)
 517     return 0;
 518
 519   Int_t iSupMod = -1;
 520   Int_t iTower  = -1;
 521   Int_t iIphi   = -1;
 522   Int_t iIeta   = -1;
 523   Int_t iphi    = -1;
 524   Int_t ieta    = -1;
 525   Int_t iphis   = -1;
 526   Int_t ietas   = -1;
 527
 528   Double_t crossEnergy = 0;
 529
 530   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 531   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 532
 533   Int_t ncells = cluster->GetNCells();
 534   for (Int_t i=0; i<ncells; i++) {
 535     Int_t cellAbsId = cluster->GetCellAbsId(i);
 536     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 537     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 538     Int_t aphidiff = TMath::Abs(iphi-iphis);
 539     if (aphidiff>1)
 540       continue;
 541     Int_t aetadiff = TMath::Abs(ieta-ietas);
 542     if (aetadiff>1)
 543       continue;
 544     if ( (aphidiff==1 && aetadiff==0) ||
 545         (aphidiff==0 && aetadiff==1) ) {
 546       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 547     }
 548   }
 549
 550   return crossEnergy;
 551 }
 552
 553 //________________________________________________________________________
 554 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 555 {
 556   // Get maximum energy of attached cell.
 557
 558   id = -1;
 559
 560   AliVCaloCells *cells = 0;
 561   cells = fESD->GetEMCALCells();
 562   if (!cells)
 563     return 0;
 564
 565   Double_t maxe = 0;
 566   Int_t ncells = cluster->GetNCells();
 567   for (Int_t i=0; i<ncells; i++) {
 568     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 569     if (e>maxe) {
 570       maxe = e;
 571       id   = cluster->GetCellAbsId(i);
 572     }
 573   }
 574   return maxe;
 575 }
 576
 577 //________________________________________________________________________
 578 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
 579 {
 580   if(!fStack)
 581     return;
 582   Int_t nTracks = fStack->GetNtrack();
 583   if(fDebug)
 584     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
 585   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 586     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 587     if(!mcp)
 588       continue;
 589     Int_t pdg = mcp->GetPdgCode();
 590     if(pdg!=22)
 591       continue;
 592     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
 593       continue;
 594     Int_t imom = mcp->GetMother(0);
 595     if(imom<0 || imom>nTracks)
 596       continue;
 597     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
 598     if(!mcmom)
 599       continue;
 600     Int_t pdgMom = mcmom->GetPdgCode();
 601     if((imom==6 || imom==7) && pdgMom==22) {
 602       fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 603       if(fNClusForDirPho==0)
 604         fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 605       if(fDebug){
 606         printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
 607         printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
 608       }
 609     }
 610     else{
 611       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
 612         fDecayPhotonPtMC->Fill(mcp->Pt());
 613     }
 614   }
 615 }
 616 //________________________________________________________________________
 617 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
 618 {
 619   if(!c)
 620     return -1;
 621   if(!fStack)
 622     return -1;
 623   Int_t nmcp = fStack->GetNtrack();
 624   Int_t clabel = c->GetLabel();
 625   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
 626     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
 627   if(!fMcIdFamily.Contains(Form("%d",clabel)))
 628     return -1;
 629   fNClusForDirPho++;
 630   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
 631   Int_t partonpos = partonposstr.Atoi();
 632   if(fDebug)
 633     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
 634   if(clabel<0 || clabel>nmcp)
 635     return 0;
 636   Float_t clsPos[3] = {0,0,0};
 637   c->GetPosition(clsPos);
 638   TVector3 clsVec(clsPos);
 639   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 640   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
 641   if(!mcp)
 642     return -1;
 643   if(fDebug){
 644     printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),clabel);
 645   }
 646   Int_t lab1 =  mcp->GetFirstDaughter();
 647   if(lab1<0 || lab1>nmcp)
 648     return -1;
 649   TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
 650   if(!mcd)
 651     return -1;
 652   if(fDebug)
 653     printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->Energy(),mcd->GetPdgCode(),lab1);
 654   if(mcd->GetPdgCode()==22){
 655     fClusEtMcPt->Fill(mcd->Pt(), Et);
 656     fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
 657   }
 658   else{
 659     printf("Warning: daughter of photon parton is not a photon\n");
 660     return -1;
 661   }
 662   return fDirPhoPt;
 663 }
 664 //________________________________________________________________________
 665 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
 666 {
 667   if(!fStack)
 668     return;
 669   Int_t selfid = 6;
 670   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
 671   if(!mcp)
 672     return;
 673   if(mcp->GetPdgCode()!=22){
 674     mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
 675     if(!mcp)
 676       return;
 677   }
 678   Int_t daug0f =  mcp->GetFirstDaughter();
 679   Int_t daug0l =  mcp->GetLastDaughter();
 680   Int_t nd0 = daug0l - daug0f;
 681   if(fDebug)
 682     printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),nd0+1);
 683   fMcIdFamily = Form("%d,",selfid);
 684   GetDaughtersInfo(daug0f,daug0l, selfid,"");
 685   if(fDebug){
 686     printf("\t---- end of the prompt  gamma's family tree ----\n\n");
 687     printf("Family id string = %s,\n\n",fMcIdFamily.Data());
 688   }
 689   TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
 690   if(!md)
 691     return;
 692   fDirPhoPt = md->Pt();
 693 }
 694 //________________________________________________________________________
 695 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(const int firstd, const int lastd, const int selfid, const char *inputind)
 696 {
 697   int nmcp = fStack->GetNtrack();
 698   if(firstd<0 || firstd>nmcp)
 699     return;
 700   if(lastd<0 || lastd>nmcp)
 701     return;
 702   if(firstd>lastd){
 703     printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
 704     return;
 705   }
 706   TString indenter = Form("\t%s",inputind);
 707   TParticle *dp = 0x0;
 708   if(fDebug)
 709     printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
 710   for(int id=firstd; id<lastd+1; id++){
 711     dp =  static_cast<TParticle*>(fStack->Particle(id));
 712     if(!dp)
 713       continue;
 714     Int_t dfd = dp->GetFirstDaughter();
 715     Int_t dld = dp->GetLastDaughter();
 716     Int_t dnd =  dld - dfd + 1;
 717     Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
 718     if(fDebug)
 719       printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->Energy(), vr, dp->GetPdgCode(), dnd, dfd, dld);
 720     fMcIdFamily += Form("%d,",id);
 721     GetDaughtersInfo(dfd,dld,id,indenter.Data());
 722   }
 723 }
 724 //________________________________________________________________________
 725 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 726 {
 727   // Called once at the end of the query.
 728 }