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 "AliAODEvent.h"
  28 #include "AliAODTrack.h"
  29 #include "AliMCEvent.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliStack.h"
  32 #include "AliVEvent.h"
  33 #include "AliVTrack.h"
  34 #include "AliV0vertexer.h"
  35 #include "AliVCluster.h"
  36 #include "AliOADBContainer.h"
  37
  38 #include <iostream>
  39 using std::cout;
  40 using std::endl;
  41
  42 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  46   AliAnalysisTaskSE(),
  47   fESDClusters(0),
  48   fAODClusters(0),
  49   fSelPrimTracks(0),
  50   fTracks(0),
  51   fESDCells(0),
  52   fAODCells(0),
  53   fPrTrCuts(0),
  54   fGeom(0x0),
  55   fGeoName("EMCAL_COMPLETEV1"),
  56   fOADBContainer(0),
  57   fPeriod("LHC11c"),
  58   fTrigBit("kEMC7"),
  59   fIsTrain(0),
  60   fIsMc(0),
  61   fDebug(0),
  62   fPathStrOpt("/"),
  63   fExoticCut(0.97),
  64   fIsoConeR(0.4),
  65   fNDimensions(7),
  66   fECut(3.),
  67   fTrackMult(0),
  68   fMcIdFamily(""),
  69   fNClusForDirPho(0),
  70   fDirPhoPt(0),
  71   fHigherPtCone(0),
  72   fImportGeometryFromFile(0),
  73   fImportGeometryFilePath(""),
  74   fMaxPtTrack(0),
  75   fMaxEClus(0),
  76   fNCells50(0),
  77   fESD(0),
  78   fAOD(0),
  79   fMCEvent(0),
  80   fStack(0),
  81   fOutputList(0),
  82   fEvtSel(0),
  83   fNClusEt10(0),
  84   fRecoPV(0),
  85   fPVtxZ(0),
  86   fTrMultDist(0),
  87   fMCDirPhotonPtEtaPhi(0),
  88   fMCIsoDirPhotonPtEtaPhi(0),
  89   fDecayPhotonPtMC(0),
  90   fCellAbsIdVsAmpl(0),
  91   fNClusHighClusE(0),
  92   fHigherPtConeM02(0),
  93   fClusEtMcPt(0),
  94   fClusMcDetaDphi(0),
  95   fNClusPerPho(0),
  96   fMcPtInConeBG(0),
  97   fMcPtInConeSBG(0),
  98   fMcPtInConeBGnoUE(0),
  99   fMcPtInConeSBGnoUE(0),
 100   fAllIsoEtMcGamma(0),
 101   fAllIsoNoUeEtMcGamma(0),
 102   fMCDirPhotonPtEtaPhiNoClus(0),
 103   fHnOutput(0),
 104   fQAList(0),
 105   fNTracks(0),
 106   fEmcNCells(0),
 107   fEmcNClus(0),
 108   fEmcNClusCut(0),
 109   fNTracksECut(0),
 110   fEmcNCellsCut(0),
 111   fEmcClusEPhi(0),
 112   fEmcClusEPhiCut(0),
 113   fEmcClusEEta(0),
 114   fEmcClusEEtaCut(0),
 115   fTrackPtPhi(0),
 116   fTrackPtPhiCut(0),
 117   fTrackPtEta(0),
 118   fTrackPtEtaCut(0),
 119   fMaxCellEPhi(0)
 120 {
 121   // Default constructor.
 122   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 123 }
 124
 125 //________________________________________________________________________
 126 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
 127   AliAnalysisTaskSE(name),
 128   fESDClusters(0),
 129   fAODClusters(0),
 130   fSelPrimTracks(0),
 131   fTracks(0),
 132   fESDCells(0),
 133   fAODCells(0),
 134   fPrTrCuts(0),
 135   fGeom(0x0),
 136   fGeoName("EMCAL_COMPLETEV1"),
 137   fOADBContainer(0),
 138   fPeriod("LHC11c"),
 139   fTrigBit("kEMC7"),
 140   fIsTrain(0),
 141   fIsMc(0),
 142   fDebug(0),
 143   fPathStrOpt("/"),
 144   fExoticCut(0.97),
 145   fIsoConeR(0.4),
 146   fNDimensions(7),
 147   fECut(3.),
 148   fTrackMult(0),
 149   fMcIdFamily(""),
 150   fNClusForDirPho(0),
 151   fDirPhoPt(0),
 152   fHigherPtCone(0),
 153   fImportGeometryFromFile(0),
 154   fImportGeometryFilePath(""),
 155   fMaxPtTrack(0),
 156   fMaxEClus(0),
 157   fNCells50(0),
 158   fESD(0),
 159   fAOD(0),
 160   fMCEvent(0),
 161   fStack(0),
 162   fOutputList(0),
 163   fEvtSel(0),
 164   fNClusEt10(0),
 165   fRecoPV(0),
 166   fPVtxZ(0),
 167   fTrMultDist(0),
 168   fMCDirPhotonPtEtaPhi(0),
 169   fMCIsoDirPhotonPtEtaPhi(0),
 170   fDecayPhotonPtMC(0),
 171   fCellAbsIdVsAmpl(0),
 172   fNClusHighClusE(0),
 173   fHigherPtConeM02(0),
 174   fClusEtMcPt(0),
 175   fClusMcDetaDphi(0),
 176   fNClusPerPho(0),
 177   fMcPtInConeBG(0),
 178   fMcPtInConeSBG(0),
 179   fMcPtInConeBGnoUE(0),
 180   fMcPtInConeSBGnoUE(0),
 181   fAllIsoEtMcGamma(0),
 182   fAllIsoNoUeEtMcGamma(0),
 183   fMCDirPhotonPtEtaPhiNoClus(0),
 184   fHnOutput(0),
 185   fQAList(0),
 186   fNTracks(0),
 187   fEmcNCells(0),
 188   fEmcNClus(0),
 189   fEmcNClusCut(0),
 190   fNTracksECut(0),
 191   fEmcNCellsCut(0),
 192   fEmcClusEPhi(0),
 193   fEmcClusEPhiCut(0),
 194   fEmcClusEEta(0),
 195   fEmcClusEEtaCut(0),
 196   fTrackPtPhi(0),
 197   fTrackPtPhiCut(0),
 198   fTrackPtEta(0),
 199   fTrackPtEtaCut(0),
 200   fMaxCellEPhi(0)
 201 {
 202   // Constructor
 203
 204   // Define input and output slots here
 205   // Input slot #0 works with a TChain
 206   DefineInput(0, TChain::Class());
 207   // Output slot #0 id reserved by the base class for AOD
 208   // Output slot #1 writes into a TH1 container
 209   DefineOutput(1, TList::Class());
 210   DefineOutput(2, TList::Class());
 211 }
 212
 213 //________________________________________________________________________
 214 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 215 {
 216   // Create histograms, called once.
 217
 218   fESDClusters = new TObjArray();
 219   fSelPrimTracks = new TObjArray();
 220
 221
 222   fOutputList = new TList();
 223   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 224
 225   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 226   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
 227   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
 228
 229   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 230   fOutputList->Add(fEvtSel);
 231
 232   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 233   fOutputList->Add(fNClusEt10);
 234
 235   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 236   fOutputList->Add(fRecoPV);
 237
 238   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 239   fOutputList->Add(fPVtxZ);
 240
 241   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
 242   fOutputList->Add(fTrMultDist);
 243
 244   fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",101,-1.5,99.5,154,-0.77,0.77,130,1.38,3.20);
 245   fMCDirPhotonPtEtaPhi->Sumw2();
 246   fOutputList->Add(fMCDirPhotonPtEtaPhi);
 247
 248   fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",101,-1.5,99.5,154,-0.77,0.77,130,1.38,3.20);
 249   fMCIsoDirPhotonPtEtaPhi->Sumw2();
 250   fOutputList->Add(fMCIsoDirPhotonPtEtaPhi);
 251
 252   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
 253   fDecayPhotonPtMC->Sumw2();
 254   fOutputList->Add(fDecayPhotonPtMC);
 255
 256   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);
 257   fOutputList->Add(fCellAbsIdVsAmpl);
 258
 259   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);
 260   fOutputList->Add(fNClusHighClusE);
 261
 262   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);
 263   fOutputList->Add(fHigherPtConeM02);
 264
 265   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
 266   fOutputList->Add(fClusEtMcPt);
 267
 268   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
 269   fOutputList->Add(fClusMcDetaDphi);
 270
 271   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
 272   fOutputList->Add(fNClusPerPho);
 273
 274   fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 275   fOutputList->Add(fMcPtInConeBG);
 276
 277   fMcPtInConeSBG  = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 278   fOutputList->Add(fMcPtInConeSBG);
 279
 280   fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 281   fOutputList->Add(fMcPtInConeBGnoUE);
 282
 283   fMcPtInConeSBGnoUE  = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100);
 284   fOutputList->Add(fMcPtInConeSBGnoUE);
 285
 286   fAllIsoEtMcGamma  = new TH2F("hAllIsoEtMcGammaE","ISO^{TRK+EMC} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC} (GeV);",1000,0,100,600,-10,50);
 287   fOutputList->Add(fAllIsoEtMcGamma);
 288
 289   fAllIsoNoUeEtMcGamma  = new TH2F("hAllIsoNoUeEtMcGammaE","ISO^{TRK+EMC}_{noue} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC}_{noue} (GeV);",1000,0,100,600,-10,50);
 290   fOutputList->Add(fAllIsoNoUeEtMcGamma);
 291
 292
 293   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);
 294   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
 295
 296   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=100;
 297   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt};
 298   fNDimensions = sizeof(bins)/sizeof(Int_t);
 299   const Int_t ndims =   fNDimensions;
 300   Double_t xmin[] = { -0.5,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-0.5};
 301   Double_t xmax[] = { 99.5, 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5};
 302   if(fPeriod.Contains("11h")){
 303     xmax[12]=3999.5;
 304   }
 305   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);
 306   fOutputList->Add(fHnOutput);
 307
 308   //QA outputs
 309   fQAList = new TList();
 310   fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 311
 312   fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5);
 313   fNTracks->Sumw2();
 314   fQAList->Add(fNTracks);
 315
 316   fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5);
 317   fEmcNCells->Sumw2();
 318   fQAList->Add(fEmcNCells);
 319   fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5);
 320   fEmcNClus->Sumw2();
 321   fQAList->Add(fEmcNClus);
 322   fEmcNClusCut = new TH1F("fEmcNClusCut",";n/event;count",120,-0.5,119.5);
 323   fEmcNClusCut->Sumw2();
 324   fQAList->Add(fEmcNClusCut);
 325   fNTracksECut = new TH1F("fNTracksECut",";n/event;count",120,-0.5,119.5);
 326   fNTracksECut->Sumw2();
 327   fQAList->Add(fNTracksECut);
 328   fEmcNCellsCut = new TH1F("fEmcNCellsCut",";n/event;count",120,-0.5,119.5);
 329   fEmcNCellsCut->Sumw2();
 330   fQAList->Add(fEmcNCellsCut);
 331   fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",100,-0.25,49.75,63,0,6.3);
 332   fEmcClusEPhi->Sumw2();
 333   fQAList->Add(fEmcClusEPhi);
 334   fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",";GeV;#phi",100,-0.25,49.75,63,0,6.3);
 335   fEmcClusEPhiCut->Sumw2();
 336   fQAList->Add(fEmcClusEPhiCut);
 337   fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",100,-0.25,49.75,19,-0.9,0.9);
 338   fEmcClusEEta->Sumw2();
 339   fQAList->Add(fEmcClusEEta);
 340   fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",";GeV;#eta",100,-0.25,49.75,18,-0.9,0.9);
 341   fEmcClusEEtaCut->Sumw2();
 342   fQAList->Add(fEmcClusEEtaCut);
 343   fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3);
 344   fTrackPtPhi->Sumw2();
 345   fQAList->Add(fTrackPtPhi);
 346   fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3);
 347   fTrackPtPhiCut->Sumw2();
 348   fQAList->Add(fTrackPtPhiCut);
 349   fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9);
 350   fTrackPtEta->Sumw2();
 351   fQAList->Add(fTrackPtEta);
 352   fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9);
 353   fTrackPtEtaCut->Sumw2();
 354   fQAList->Add(fTrackPtEtaCut);
 355
 356
 357   fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",100,-0.25,49.75,63,0,6.3);
 358   fMaxCellEPhi->Sumw2();
 359   fQAList->Add(fMaxCellEPhi);
 360
 361   PostData(1, fOutputList);
 362   PostData(2, fQAList);
 363 }
 364
 365 //________________________________________________________________________
 366 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 367 {
 368   // Main loop, called for each event.
 369   fESDClusters = 0;
 370   fESDCells = 0;
 371   fAODClusters = 0;
 372   fAODCells = 0;
 373   // event trigger selection
 374   Bool_t isSelected = 0;
 375   if(fPeriod.Contains("11a")){
 376     if(fTrigBit.Contains("kEMC"))
 377       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 378     else
 379       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 380   }
 381   else{
 382     if(fTrigBit.Contains("kEMC"))
 383       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 384     else
 385       if(fTrigBit.Contains("kMB"))
 386         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 387       else
 388         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 389   }
 390   if(fPeriod.Contains("11h")){
 391     if(fTrigBit.Contains("kAny"))
 392       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
 393     if(fTrigBit.Contains("kAnyINT"))
 394       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
 395   }
 396   if(fIsMc)
 397     isSelected = kTRUE;
 398   if(fDebug)
 399     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 400   if(!isSelected )
 401         return;
 402   if(fIsMc){
 403     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 404     TFile *file = (TFile*)tree->GetCurrentFile();
 405     TString filename = file->GetName();
 406     if(!filename.Contains(fPathStrOpt.Data()))
 407       return;
 408   }
 409   AliVEvent *event = (AliVEvent*)InputEvent();
 410   if (!event) {
 411     printf("ERROR: event not available\n");
 412     return;
 413   }
 414   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 415   if(fDebug)
 416     printf("run number = %d\n",runnumber);
 417
 418   fESD = dynamic_cast<AliESDEvent*>(event);
 419   if(!fESD){
 420     fAOD = dynamic_cast<AliAODEvent*>(event);
 421     if(!fAOD){
 422       printf("ERROR: Invalid type of event!!!\n");
 423       return;
 424     }
 425     else if(fDebug)
 426       printf("AOD EVENT!\n");
 427   }
 428
 429   fEvtSel->Fill(0);
 430   if(fDebug)
 431     printf("event is ok,\n run number=%d\n",runnumber);
 432
 433
 434   AliVVertex *pv = (AliVVertex*)event->GetPrimaryVertex();
 435   Bool_t pvStatus = kTRUE;
 436   if(fESD){
 437     AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
 438     pvStatus = esdv->GetStatus();
 439   }
 440   if(!pv)
 441     return;
 442   if(!pvStatus)
 443     fRecoPV->Fill(0);
 444   else
 445     fRecoPV->Fill(1);
 446   fPVtxZ->Fill(pv->GetZ());
 447   if(TMath::Abs(pv->GetZ())>15)
 448     return;
 449   if(fDebug)
 450     printf("passed vertex cut\n");
 451
 452   fEvtSel->Fill(1);
 453   if(fESD)
 454     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 455   if(fAOD)
 456     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 457
 458   if(!fTracks){
 459     AliError("Track array in event is NULL!");
 460     if(fDebug)
 461       printf("returning due to not finding tracks!\n");
 462     return;
 463   }
 464   // Track loop to fill a pT spectrum
 465   const Int_t Ntracks = fTracks->GetEntriesFast();
 466   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 467     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 468     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 469     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 470     if (!track)
 471       continue;
 472     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 473     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 474     if (esdTrack && fPrTrCuts && fPrTrCuts->IsSelected(track)){
 475       fSelPrimTracks->Add(track);
 476       /*if(fTrackMaxPt<track->Pt())
 477         fTrackMaxPt = track->Pt();*/
 478       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 479     }
 480     else if(aodTrack){
 481       fSelPrimTracks->Add(track);
 482       /*if(fTrackMaxPt<track->Pt())
 483         fTrackMaxPt = track->Pt();*/
 484     }
 485   }
 486
 487   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 488   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 489     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 490       break;
 491     /*if(fESD)
 492       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 493       else*/
 494     // if(event->GetEMCALMatrix(mod))
 495     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 496     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 497   }
 498
 499   if(fESD){
 500     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 501     fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 502     if(fDebug)
 503       printf("ESD Track mult= %d\n",fTrackMult);
 504   }
 505   else if(fAOD){
 506     fTrackMult = Ntracks;
 507     if(fDebug)
 508       printf("AOD Track mult= %d\n",fTrackMult);
 509   }
 510   fTrMultDist->Fill(fTrackMult);
 511
 512   if(fESD){
 513     TList *l = fESD->GetList();
 514     fESDClusters =  dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 515     fESDCells = fESD->GetEMCALCells();
 516     if(fDebug)
 517       printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
 518   }
 519   else if(fAOD){
 520     fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 521     fAODCells = fAOD->GetEMCALCells();
 522     if(fDebug)
 523       printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
 524   }
 525
 526
 527   fMCEvent = MCEvent();
 528   if(fMCEvent){
 529     if(fDebug)
 530       std::cout<<"MCevent exists\n";
 531     fStack = (AliStack*)fMCEvent->Stack();
 532   }
 533   else{
 534     if(fDebug)
 535       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 536   }
 537   LoopOnCells();
 538   FollowGamma();
 539   if(fDebug)
 540     printf("passed calling of FollowGamma\n");
 541   FillClusHists();
 542   if(fDebug)
 543     printf("passed calling of FillClusHists\n");
 544   FillMcHists();
 545   if(fDebug)
 546     printf("passed calling of FillMcHists\n");
 547   FillQA();
 548   if(fDebug)
 549     printf("passed calling of FillQA\n");
 550   /*if(fESD)
 551     fESDClusters->Clear();*/
 552   fSelPrimTracks->Clear();
 553   fNClusForDirPho = 0;
 554   fNCells50 = 0;
 555
 556   PostData(1, fOutputList);
 557   PostData(2, fQAList);
 558 }
 559
 560 //________________________________________________________________________
 561 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 562 {
 563   if(fDebug)
 564     printf("Inside FillClusHists()....\n");
 565   // Fill cluster histograms.
 566   TObjArray *clusters = fESDClusters;
 567
 568   if (!clusters){
 569     clusters = fAODClusters;
 570     if(fDebug)
 571       printf("ESD clusters empty...");
 572   }
 573   if (!clusters){
 574     if(fDebug)
 575       printf("  and AOD clusters as well!!!\n");
 576     return;
 577   }
 578   if(fDebug)
 579     printf("\n");
 580
 581   const Int_t nclus = clusters->GetEntries();
 582   if(nclus==0)
 583     return;
 584   if(fDebug)
 585     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 586   Double_t maxE = 0;
 587   Int_t nclus10 = 0;
 588   Double_t ptmc=-1;
 589   for(Int_t ic=0;ic<nclus;ic++){
 590     maxE=0;
 591     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 592     if(!c)
 593       continue;
 594     if(!c->IsEMCAL())
 595       continue;
 596     if(c->E()<fECut)
 597       continue;
 598     Short_t id;
 599     Double_t Emax = GetMaxCellEnergy( c, id);
 600     Double_t Ecross = GetCrossEnergy( c, id);
 601     if((1-Ecross/Emax)>fExoticCut)
 602       continue;
 603     Float_t clsPos[3] = {0,0,0};
 604     c->GetPosition(clsPos);
 605     TVector3 clsVec(clsPos);
 606     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 607     if(fDebug)
 608       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 609     if(Et>10)
 610       nclus10++;
 611     Float_t ceiso, cephiband, cecore;
 612     Float_t triso, trphiband, trcore;
 613     Float_t alliso, allphiband, allcore;
 614     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 615     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 616     GetCeIso(clsVec, id, ceiso, cephiband, cecore);
 617     GetTrIso(clsVec, triso, trphiband, trcore);
 618     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 619     if(Et>10 && Et<15 && dr>0.025){
 620       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 621       if(fDebug)
 622         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 623     }
 624     alliso = ceiso + triso;
 625     allphiband = cephiband + trphiband;
 626     allcore = cecore + trcore;
 627     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 628     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 629     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 630     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
 631     if(fDebug && Et>10)
 632       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
 633     if(c->GetM02()>0.5 && c->GetM02()<2.0){
 634       fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum);
 635       fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum);
 636     }
 637     if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
 638       fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum);
 639       fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum);
 640       if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
 641         fAllIsoEtMcGamma->Fill(Et, alliso-/*Et*/cecore-allisoue);
 642         fAllIsoNoUeEtMcGamma->Fill(Et, alliso-/*Et*/cecore);
 643       }
 644     }
 645     const Int_t ndims =   fNDimensions;
 646     Double_t outputValues[ndims];
 647     if(mcptsum<2)
 648       ptmc = GetClusSource(c);
 649     else
 650       ptmc = -0.1;
 651     outputValues[0] = Et;
 652     outputValues[1] = c->GetM02();
 653     outputValues[2] = ceiso/*cecore*/-ceisoue;
 654     outputValues[3] = triso-trisoue;
 655     outputValues[4] = alliso/*cecore*/-allisoue - trcore;
 656     outputValues[5] = ceiso;
 657     outputValues[6] = alliso - trcore;
 658     outputValues[7] = c->GetTrackDx();
 659     outputValues[8] = c->GetTrackDz();
 660     outputValues[9] = clsVec.Eta();
 661     outputValues[10] = clsVec.Phi();
 662     if(fESDCells)
 663       outputValues[11] = fESDCells->GetCellTime(id);
 664     else if(fAODCells)
 665       outputValues[11] = fAODCells->GetCellTime(id);
 666     outputValues[12] = fTrackMult;
 667     outputValues[13] = ptmc;
 668     fHnOutput->Fill(outputValues);
 669     if(c->E()>maxE)
 670       maxE = c->E();
 671   }
 672   fNClusHighClusE->Fill(maxE,nclus);
 673   fMaxEClus = maxE;
 674   fNClusEt10->Fill(nclus10);
 675   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 676 }
 677
 678 //________________________________________________________________________
 679 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core)
 680 {
 681   // Get cell isolation.
 682   AliVCaloCells *cells = fESDCells;
 683   if (!cells){
 684     cells = fAODCells;
 685     if(fDebug)
 686       printf("ESD cells empty...");
 687   }
 688   if (!cells){
 689      if(fDebug)
 690       printf("  and AOD cells are empty  as well!!!\n");
 691     return;
 692   }
 693
 694   TObjArray *clusters = fESDClusters;
 695   if (!clusters)
 696     clusters = fAODClusters;
 697   if (!clusters)
 698     return;
 699
 700
 701   const Int_t nclus = clusters->GetEntries();
 702   //const Int_t ncells = cells->GetNumberOfCells();
 703   Float_t totiso=0;
 704   Float_t totband=0;
 705   Float_t totcore=0;
 706   Float_t etacl = vec.Eta();
 707   Float_t phicl = vec.Phi();
 708   Float_t maxtcl = cells->GetCellTime(maxid);
 709   if(phicl<0)
 710     phicl+=TMath::TwoPi();
 711   /*Int_t absid = -1;
 712   Float_t eta=-1, phi=-1;
 713   for(int icell=0;icell<ncells;icell++){
 714     absid = TMath::Abs(cells->GetCellNumber(icell));
 715     Float_t celltime = cells->GetCellTime(absid);
 716     //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
 717     if(TMath::Abs(celltime-maxtcl)>2e-8 )
 718       continue;
 719     if(!fGeom)
 720       return;
 721     fGeom->EtaPhiFromIndex(absid,eta,phi);
 722     Float_t dphi = TMath::Abs(phi-phicl);
 723     Float_t deta = TMath::Abs(eta-etacl);
 724     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/
 725   for(int ic=0;ic<nclus;ic++){
 726     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 727     if(!c)
 728       continue;
 729     if(!c->IsEMCAL())
 730       continue;
 731     Short_t id;
 732     GetMaxCellEnergy( c, id);
 733     Double_t maxct = cells->GetCellTime(id);
 734     if(TMath::Abs(maxtcl-maxct)>2.5e-9)
 735       continue;
 736     Float_t clsPos[3] = {0,0,0};
 737     c->GetPosition(clsPos);
 738     TVector3 cv(clsPos);
 739     Double_t Et = c->E()*TMath::Sin(cv.Theta());
 740     Float_t dphi = TMath::Abs(cv.Phi()-phicl);
 741     Float_t deta = TMath::Abs(cv.Eta()-etacl);
 742     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 743     if(R<0.007)
 744       continue;
 745     Double_t matchedpt =  GetTrackMatchedPt(c->GetTrackMatchedIndex());
 746     Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
 747     if(nEt<0)
 748       printf("nEt=%1.1f\n",nEt);
 749     if(R<fIsoConeR){
 750       totiso += nEt;
 751       if(R<0.04)
 752         totcore += nEt;
 753     }
 754     else{
 755       if(dphi>fIsoConeR)
 756         continue;
 757       if(deta<fIsoConeR)
 758         continue;
 759       totband += nEt;
 760     }
 761   }
 762   iso = totiso;
 763   phiband = totband;
 764   core = totcore;
 765 }
 766 //________________________________________________________________________
 767 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 768 {
 769   // Get track isolation.
 770
 771   if(!fSelPrimTracks)
 772     return;
 773   fHigherPtCone = 0;
 774   const Int_t ntracks = fSelPrimTracks->GetEntries();
 775   Double_t totiso=0;
 776   Double_t totband=0;
 777   Double_t totcore=0;
 778   Double_t etacl = vec.Eta();
 779   Double_t phicl = vec.Phi();
 780   if(phicl<0)
 781     phicl+=TMath::TwoPi();
 782   for(int itrack=0;itrack<ntracks;itrack++){
 783     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 784     if(!track)
 785       continue;
 786     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 787     Double_t deta = TMath::Abs(track->Eta()-etacl);
 788     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 789     Double_t pt = track->Pt();
 790     if(pt>fHigherPtCone)
 791       fHigherPtCone = pt;
 792     if(R<fIsoConeR){
 793       totiso += track->Pt();
 794       if(R<0.04)
 795         totcore += pt;
 796     }
 797     else{
 798       if(dphi>fIsoConeR)
 799         continue;
 800       if(deta<fIsoConeR)
 801         continue;
 802       totband += track->Pt();
 803     }
 804   }
 805   iso = totiso;
 806   phiband = totband;
 807   core = totcore;
 808 }
 809
 810 //________________________________________________________________________
 811 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 812 {
 813   // Calculate the energy of cross cells around the leading cell.
 814
 815   AliVCaloCells *cells = 0;
 816   cells = fESDCells;
 817   if (!cells)
 818     cells = fAODCells;
 819   if (!cells)
 820     return 0;
 821
 822   if (!fGeom)
 823     return 0;
 824
 825   Int_t iSupMod = -1;
 826   Int_t iTower  = -1;
 827   Int_t iIphi   = -1;
 828   Int_t iIeta   = -1;
 829   Int_t iphi    = -1;
 830   Int_t ieta    = -1;
 831   Int_t iphis   = -1;
 832   Int_t ietas   = -1;
 833
 834   Double_t crossEnergy = 0;
 835
 836   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 837   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 838
 839   Int_t ncells = cluster->GetNCells();
 840   for (Int_t i=0; i<ncells; i++) {
 841     Int_t cellAbsId = cluster->GetCellAbsId(i);
 842     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 843     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 844     Int_t aphidiff = TMath::Abs(iphi-iphis);
 845     if (aphidiff>1)
 846       continue;
 847     Int_t aetadiff = TMath::Abs(ieta-ietas);
 848     if (aetadiff>1)
 849       continue;
 850     if ( (aphidiff==1 && aetadiff==0) ||
 851         (aphidiff==0 && aetadiff==1) ) {
 852       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 853     }
 854   }
 855
 856   return crossEnergy;
 857 }
 858
 859 //________________________________________________________________________
 860 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 861 {
 862   // Get maximum energy of attached cell.
 863
 864   id = -1;
 865
 866   AliVCaloCells *cells = 0;
 867   cells = fESDCells;
 868   if (!cells)
 869     cells = fAODCells;
 870   if(!cells)
 871     return 0;
 872
 873   Double_t maxe = 0;
 874   Int_t ncells = cluster->GetNCells();
 875   for (Int_t i=0; i<ncells; i++) {
 876     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 877     if (e>maxe) {
 878       maxe = e;
 879       id   = cluster->GetCellAbsId(i);
 880     }
 881   }
 882   return maxe;
 883 }
 884
 885 //________________________________________________________________________
 886 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
 887 {
 888   if(!fStack)
 889     return;
 890   Int_t nTracks = fStack->GetNtrack();
 891   if(fDebug)
 892     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
 893   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 894     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 895     if(!mcp)
 896       continue;
 897     Int_t pdg = mcp->GetPdgCode();
 898     if(pdg!=22)
 899       continue;
 900     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
 901       continue;
 902     Int_t imom = mcp->GetMother(0);
 903     if(imom<0 || imom>nTracks)
 904       continue;
 905     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
 906     if(!mcmom)
 907       continue;
 908     Int_t pdgMom = mcmom->GetPdgCode();
 909     if((imom==6 || imom==7) && pdgMom==22) {
 910       fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 911       Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), 0.4);
 912       if(ptsum<2)
 913         fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 914       if(fNClusForDirPho==0)
 915         fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 916       if(fDebug){
 917         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());
 918         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());
 919       }
 920     }
 921     else{
 922       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
 923         fDecayPhotonPtMC->Fill(mcp->Pt());
 924     }
 925   }
 926 }
 927 //________________________________________________________________________
 928 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
 929 {
 930   if(!c)
 931     return -0.1;
 932   if(!fStack)
 933     return -0.1;
 934   Int_t nmcp = fStack->GetNtrack();
 935   Int_t clabel = c->GetLabel();
 936   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
 937     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
 938   if(!fMcIdFamily.Contains(Form("%d",clabel)))
 939     return -0.1;
 940   fNClusForDirPho++;
 941   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
 942   Int_t partonpos = partonposstr.Atoi();
 943   if(fDebug)
 944     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
 945   if(clabel<0 || clabel>nmcp)
 946     return 0;
 947   Float_t clsPos[3] = {0,0,0};
 948   c->GetPosition(clsPos);
 949   TVector3 clsVec(clsPos);
 950   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 951   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
 952   if(!mcp)
 953     return -0.1;
 954   if(fDebug){
 955     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);
 956   }
 957   Int_t lab1 =  mcp->GetFirstDaughter();
 958   if(lab1<0 || lab1>nmcp)
 959     return -0.1;
 960   TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
 961   if(!mcd)
 962     return -0.1;
 963   if(fDebug)
 964     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);
 965   if(mcd->GetPdgCode()==22){
 966     fClusEtMcPt->Fill(mcd->Pt(), Et);
 967     fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
 968   }
 969   else{
 970     printf("Warning: daughter of photon parton is not a photon\n");
 971     return -0.1;
 972   }
 973   return fDirPhoPt;
 974 }
 975 //________________________________________________________________________
 976 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
 977 {
 978   if(!fStack)
 979     return;
 980   Int_t selfid = 6;
 981   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
 982   if(!mcp)
 983     return;
 984   if(mcp->GetPdgCode()!=22){
 985     mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
 986     if(!mcp)
 987       return;
 988   }
 989   Int_t daug0f =  mcp->GetFirstDaughter();
 990   Int_t daug0l =  mcp->GetLastDaughter();
 991   Int_t nd0 = daug0l - daug0f;
 992   if(fDebug)
 993     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);
 994   fMcIdFamily = Form("%d,",selfid);
 995   GetDaughtersInfo(daug0f,daug0l, selfid,"");
 996   if(fDebug){
 997     printf("\t---- end of the prompt  gamma's family tree ----\n\n");
 998     printf("Family id string = %s,\n\n",fMcIdFamily.Data());
 999   }
1000   TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1001   if(!md)
1002     return;
1003   fDirPhoPt = md->Pt();
1004 }
1005 //________________________________________________________________________
1006 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
1007 {
1008   int nmcp = fStack->GetNtrack();
1009   if(firstd<0 || firstd>nmcp)
1010     return;
1011   if(lastd<0 || lastd>nmcp)
1012     return;
1013   if(firstd>lastd){
1014     printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1015     return;
1016   }
1017   TString indenter = Form("\t%s",inputind);
1018   TParticle *dp = 0x0;
1019   if(fDebug)
1020     printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1021   for(int id=firstd; id<lastd+1; id++){
1022     dp =  static_cast<TParticle*>(fStack->Particle(id));
1023     if(!dp)
1024       continue;
1025     Int_t dfd = dp->GetFirstDaughter();
1026     Int_t dld = dp->GetLastDaughter();
1027     Int_t dnd =  dld - dfd + 1;
1028     Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1029     if(fDebug)
1030       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);
1031     fMcIdFamily += Form("%d,",id);
1032     GetDaughtersInfo(dfd,dld,id,indenter.Data());
1033   }
1034 }
1035
1036 //________________________________________________________________________
1037 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1038 {
1039   if(!fStack)
1040     return 0;
1041   if(fDebug)
1042     printf("Inside GetMcPtSumInCone!!\n");
1043   Int_t nTracks = fStack->GetNtrack();
1044   Float_t ptsum = 0;
1045   TString addpartlabels = "";
1046   for(Int_t iTrack=8;iTrack<nTracks;iTrack++){
1047     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1048     if(!mcp)
1049       continue;
1050     Int_t status = mcp->GetStatusCode();
1051     if(status!=1)
1052       continue;
1053     Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1054    if(mcvr>1)
1055       continue;
1056     else {
1057       if(fDebug)
1058         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1059     }
1060     Float_t dphi = mcp->Phi() - phiclus;
1061     Float_t deta = mcp->Eta() - etaclus;
1062     if(fDebug)
1063       printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1064     if(deta>R || dphi>R)
1065       continue;
1066     Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1067     if(dR>R)
1068       continue;
1069     addpartlabels += Form("%d,",iTrack);
1070     if(mcp->Pt()<0.2)
1071       continue;
1072     ptsum += mcp->Pt();
1073   }
1074   return ptsum;
1075 }
1076 //________________________________________________________________________
1077 void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1078 {
1079   if(!fSelPrimTracks)
1080     return;
1081   const int ntracks = fSelPrimTracks->GetEntriesFast();
1082   const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
1083   const Int_t nclus = fESDClusters->GetEntries();
1084
1085   fNTracks->Fill(ntracks);
1086   fEmcNCells->Fill(ncells);
1087   fEmcNClus->Fill(nclus);
1088   if(fMaxEClus>fECut){
1089     fNTracksECut->Fill(ntracks);
1090     fEmcNCellsCut->Fill(ncells);
1091     fEmcNClusCut->Fill(nclus);
1092   }
1093   for(int it=0;it<ntracks;it++){
1094     AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1095     if(!t)
1096       continue;
1097     fTrackPtPhi->Fill(t->Pt(),t->Phi());
1098     fTrackPtEta->Fill(t->Pt(),t->Eta());
1099     if(fMaxEClus>fECut){
1100       fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1101       fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1102     }
1103   }
1104   for(int ic=0;ic<nclus;ic++){
1105     AliVCluster *c = dynamic_cast<AliVCluster*>(fESDClusters->At(ic));
1106     //AliESDCaloCluster *c = (AliESDCaloCluster*)fESDClusters->At(ic);
1107     if(!c)
1108       continue;
1109     if(!c->IsEMCAL())
1110       continue;
1111     Float_t clsPos[3] = {0,0,0};
1112     c->GetPosition(clsPos);
1113     TVector3 clsVec(clsPos);
1114     Double_t cphi = clsVec.Phi();
1115     Double_t ceta = clsVec.Eta();
1116     fEmcClusEPhi->Fill(c->E(), cphi);
1117     fEmcClusEEta->Fill(c->E(), ceta);
1118     if(fMaxEClus>fECut){
1119       fEmcClusEPhiCut->Fill(c->E(), cphi);
1120       fEmcClusEEtaCut->Fill(c->E(), ceta);
1121     }
1122   }
1123 }
1124 //________________________________________________________________________
1125 Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1126 {
1127   Double_t pt = 0;
1128   if(!fTracks)
1129     return pt;
1130   if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1131     if(fDebug)
1132       printf("track-matched index out of track array range!!!\n");
1133     return pt;
1134   }
1135   AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1136   if(!track){
1137     if(fDebug)
1138       printf("track-matched pointer does not exist!!!\n");
1139     return pt;
1140   }
1141   if(fESD){
1142     if(!fPrTrCuts)
1143       return pt;
1144     if(!fPrTrCuts->IsSelected(track))
1145       return pt;
1146     pt = track->Pt();
1147   }
1148   return pt;
1149 }
1150 //________________________________________________________________________
1151 void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1152 {
1153   AliVCaloCells *cells = 0;
1154   cells = fESDCells;
1155   if (!cells)
1156     cells = fAODCells;
1157   if(!cells)
1158     return;
1159   Double_t maxe = 0;
1160   Double_t maxphi = -10;
1161   Int_t ncells = cells->GetNumberOfCells();
1162   Double_t eta,phi;
1163   for (Int_t i=0; i<ncells; i++) {
1164     Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1165     Double_t e = cells->GetCellAmplitude(absid);
1166     if(e>0.05)
1167       fNCells50++;
1168     else
1169       continue;
1170     fGeom->EtaPhiFromIndex(absid,eta,phi);
1171     if(maxe<e){
1172       maxe = e;
1173       maxphi = phi;
1174     }
1175   }
1176   fMaxCellEPhi->Fill(maxe,maxphi);
1177
1178 }
1179 //________________________________________________________________________
1180 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1181 {
1182   // Called once at the end of the query.
1183 }