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     if(maxid==id)
 746       continue;
 747     Double_t matchedpt =  GetTrackMatchedPt(c->GetTrackMatchedIndex());
 748     Double_t nEt = TMath::Max(Et-matchedpt, 0.0);
 749     if(nEt<0)
 750       printf("nEt=%1.1f\n",nEt);
 751     if(R<fIsoConeR){
 752       totiso += nEt;
 753       if(R<0.04)
 754         totcore += nEt;
 755     }
 756     else{
 757       if(dphi>fIsoConeR)
 758         continue;
 759       if(deta<fIsoConeR)
 760         continue;
 761       totband += nEt;
 762     }
 763   }
 764   iso = totiso;
 765   phiband = totband;
 766   core = totcore;
 767 }
 768 //________________________________________________________________________
 769 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 770 {
 771   // Get track isolation.
 772
 773   if(!fSelPrimTracks)
 774     return;
 775   fHigherPtCone = 0;
 776   const Int_t ntracks = fSelPrimTracks->GetEntries();
 777   Double_t totiso=0;
 778   Double_t totband=0;
 779   Double_t totcore=0;
 780   Double_t etacl = vec.Eta();
 781   Double_t phicl = vec.Phi();
 782   if(phicl<0)
 783     phicl+=TMath::TwoPi();
 784   for(int itrack=0;itrack<ntracks;itrack++){
 785     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 786     if(!track)
 787       continue;
 788     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 789     Double_t deta = TMath::Abs(track->Eta()-etacl);
 790     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 791     Double_t pt = track->Pt();
 792     if(pt>fHigherPtCone)
 793       fHigherPtCone = pt;
 794     if(R<fIsoConeR){
 795       totiso += track->Pt();
 796       if(R<0.04)
 797         totcore += pt;
 798     }
 799     else{
 800       if(dphi>fIsoConeR)
 801         continue;
 802       if(deta<fIsoConeR)
 803         continue;
 804       totband += track->Pt();
 805     }
 806   }
 807   iso = totiso;
 808   phiband = totband;
 809   core = totcore;
 810 }
 811
 812 //________________________________________________________________________
 813 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 814 {
 815   // Calculate the energy of cross cells around the leading cell.
 816
 817   AliVCaloCells *cells = 0;
 818   cells = fESDCells;
 819   if (!cells)
 820     cells = fAODCells;
 821   if (!cells)
 822     return 0;
 823
 824   if (!fGeom)
 825     return 0;
 826
 827   Int_t iSupMod = -1;
 828   Int_t iTower  = -1;
 829   Int_t iIphi   = -1;
 830   Int_t iIeta   = -1;
 831   Int_t iphi    = -1;
 832   Int_t ieta    = -1;
 833   Int_t iphis   = -1;
 834   Int_t ietas   = -1;
 835
 836   Double_t crossEnergy = 0;
 837
 838   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 839   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 840
 841   Int_t ncells = cluster->GetNCells();
 842   for (Int_t i=0; i<ncells; i++) {
 843     Int_t cellAbsId = cluster->GetCellAbsId(i);
 844     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 845     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 846     Int_t aphidiff = TMath::Abs(iphi-iphis);
 847     if (aphidiff>1)
 848       continue;
 849     Int_t aetadiff = TMath::Abs(ieta-ietas);
 850     if (aetadiff>1)
 851       continue;
 852     if ( (aphidiff==1 && aetadiff==0) ||
 853         (aphidiff==0 && aetadiff==1) ) {
 854       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 855     }
 856   }
 857
 858   return crossEnergy;
 859 }
 860
 861 //________________________________________________________________________
 862 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 863 {
 864   // Get maximum energy of attached cell.
 865
 866   id = -1;
 867
 868   AliVCaloCells *cells = 0;
 869   cells = fESDCells;
 870   if (!cells)
 871     cells = fAODCells;
 872   if(!cells)
 873     return 0;
 874
 875   Double_t maxe = 0;
 876   Int_t ncells = cluster->GetNCells();
 877   for (Int_t i=0; i<ncells; i++) {
 878     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 879     if (e>maxe) {
 880       maxe = e;
 881       id   = cluster->GetCellAbsId(i);
 882     }
 883   }
 884   return maxe;
 885 }
 886
 887 //________________________________________________________________________
 888 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
 889 {
 890   if(!fStack)
 891     return;
 892   Int_t nTracks = fStack->GetNtrack();
 893   if(fDebug)
 894     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
 895   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 896     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 897     if(!mcp)
 898       continue;
 899     Int_t pdg = mcp->GetPdgCode();
 900     if(pdg!=22)
 901       continue;
 902     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
 903       continue;
 904     Int_t imom = mcp->GetMother(0);
 905     if(imom<0 || imom>nTracks)
 906       continue;
 907     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
 908     if(!mcmom)
 909       continue;
 910     Int_t pdgMom = mcmom->GetPdgCode();
 911     if((imom==6 || imom==7) && pdgMom==22) {
 912       fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 913       Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), 0.4);
 914       if(ptsum<2)
 915         fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 916       if(fNClusForDirPho==0)
 917         fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 918       if(fDebug){
 919         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());
 920         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());
 921       }
 922     }
 923     else{
 924       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
 925         fDecayPhotonPtMC->Fill(mcp->Pt());
 926     }
 927   }
 928 }
 929 //________________________________________________________________________
 930 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
 931 {
 932   if(!c)
 933     return -0.1;
 934   if(!fStack)
 935     return -0.1;
 936   Int_t nmcp = fStack->GetNtrack();
 937   Int_t clabel = c->GetLabel();
 938   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
 939     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
 940   if(!fMcIdFamily.Contains(Form("%d",clabel)))
 941     return -0.1;
 942   fNClusForDirPho++;
 943   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
 944   Int_t partonpos = partonposstr.Atoi();
 945   if(fDebug)
 946     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
 947   if(clabel<0 || clabel>nmcp)
 948     return 0;
 949   Float_t clsPos[3] = {0,0,0};
 950   c->GetPosition(clsPos);
 951   TVector3 clsVec(clsPos);
 952   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 953   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
 954   if(!mcp)
 955     return -0.1;
 956   if(fDebug){
 957     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);
 958   }
 959   Int_t lab1 =  mcp->GetFirstDaughter();
 960   if(lab1<0 || lab1>nmcp)
 961     return -0.1;
 962   TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
 963   if(!mcd)
 964     return -0.1;
 965   if(fDebug)
 966     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);
 967   if(mcd->GetPdgCode()==22){
 968     fClusEtMcPt->Fill(mcd->Pt(), Et);
 969     fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
 970   }
 971   else{
 972     printf("Warning: daughter of photon parton is not a photon\n");
 973     return -0.1;
 974   }
 975   return fDirPhoPt;
 976 }
 977 //________________________________________________________________________
 978 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
 979 {
 980   if(!fStack)
 981     return;
 982   Int_t selfid = 6;
 983   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
 984   if(!mcp)
 985     return;
 986   if(mcp->GetPdgCode()!=22){
 987     mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
 988     if(!mcp)
 989       return;
 990   }
 991   Int_t daug0f =  mcp->GetFirstDaughter();
 992   Int_t daug0l =  mcp->GetLastDaughter();
 993   Int_t nd0 = daug0l - daug0f;
 994   if(fDebug)
 995     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);
 996   fMcIdFamily = Form("%d,",selfid);
 997   GetDaughtersInfo(daug0f,daug0l, selfid,"");
 998   if(fDebug){
 999     printf("\t---- end of the prompt  gamma's family tree ----\n\n");
1000     printf("Family id string = %s,\n\n",fMcIdFamily.Data());
1001   }
1002   TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
1003   if(!md)
1004     return;
1005   fDirPhoPt = md->Pt();
1006 }
1007 //________________________________________________________________________
1008 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
1009 {
1010   int nmcp = fStack->GetNtrack();
1011   if(firstd<0 || firstd>nmcp)
1012     return;
1013   if(lastd<0 || lastd>nmcp)
1014     return;
1015   if(firstd>lastd){
1016     printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
1017     return;
1018   }
1019   TString indenter = Form("\t%s",inputind);
1020   TParticle *dp = 0x0;
1021   if(fDebug)
1022     printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
1023   for(int id=firstd; id<lastd+1; id++){
1024     dp =  static_cast<TParticle*>(fStack->Particle(id));
1025     if(!dp)
1026       continue;
1027     Int_t dfd = dp->GetFirstDaughter();
1028     Int_t dld = dp->GetLastDaughter();
1029     Int_t dnd =  dld - dfd + 1;
1030     Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
1031     if(fDebug)
1032       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);
1033     fMcIdFamily += Form("%d,",id);
1034     GetDaughtersInfo(dfd,dld,id,indenter.Data());
1035   }
1036 }
1037
1038 //________________________________________________________________________
1039 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
1040 {
1041   if(!fStack)
1042     return 0;
1043   if(fDebug)
1044     printf("Inside GetMcPtSumInCone!!\n");
1045   Int_t nTracks = fStack->GetNtrack();
1046   Float_t ptsum = 0;
1047   TString addpartlabels = "";
1048   for(Int_t iTrack=8;iTrack<nTracks;iTrack++){
1049     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
1050     if(!mcp)
1051       continue;
1052     Int_t status = mcp->GetStatusCode();
1053     if(status!=1)
1054       continue;
1055     Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz());
1056    if(mcvr>1)
1057       continue;
1058     else {
1059       if(fDebug)
1060         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
1061     }
1062     Float_t dphi = mcp->Phi() - phiclus;
1063     Float_t deta = mcp->Eta() - etaclus;
1064     if(fDebug)
1065       printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
1066     if(deta>R || dphi>R)
1067       continue;
1068     Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
1069     if(dR>R)
1070       continue;
1071     addpartlabels += Form("%d,",iTrack);
1072     if(mcp->Pt()<0.2)
1073       continue;
1074     ptsum += mcp->Pt();
1075   }
1076   return ptsum;
1077 }
1078 //________________________________________________________________________
1079 void AliAnalysisTaskEMCALIsoPhoton::FillQA()
1080 {
1081   if(!fSelPrimTracks)
1082     return;
1083   const int ntracks = fSelPrimTracks->GetEntriesFast();
1084   const int ncells = fNCells50;//fESDCells->GetNumberOfCells();
1085   const Int_t nclus = fESDClusters->GetEntries();
1086
1087   fNTracks->Fill(ntracks);
1088   fEmcNCells->Fill(ncells);
1089   fEmcNClus->Fill(nclus);
1090   if(fMaxEClus>fECut){
1091     fNTracksECut->Fill(ntracks);
1092     fEmcNCellsCut->Fill(ncells);
1093     fEmcNClusCut->Fill(nclus);
1094   }
1095   for(int it=0;it<ntracks;it++){
1096     AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it);
1097     if(!t)
1098       continue;
1099     fTrackPtPhi->Fill(t->Pt(),t->Phi());
1100     fTrackPtEta->Fill(t->Pt(),t->Eta());
1101     if(fMaxEClus>fECut){
1102       fTrackPtPhiCut->Fill(t->Pt(), t->Phi());
1103       fTrackPtEtaCut->Fill(t->Pt(), t->Eta());
1104     }
1105   }
1106   for(int ic=0;ic<nclus;ic++){
1107     AliVCluster *c = dynamic_cast<AliVCluster*>(fESDClusters->At(ic));
1108     //AliESDCaloCluster *c = (AliESDCaloCluster*)fESDClusters->At(ic);
1109     if(!c)
1110       continue;
1111     if(!c->IsEMCAL())
1112       continue;
1113     Float_t clsPos[3] = {0,0,0};
1114     c->GetPosition(clsPos);
1115     TVector3 clsVec(clsPos);
1116     Double_t cphi = clsVec.Phi();
1117     Double_t ceta = clsVec.Eta();
1118     fEmcClusEPhi->Fill(c->E(), cphi);
1119     fEmcClusEEta->Fill(c->E(), ceta);
1120     if(fMaxEClus>fECut){
1121       fEmcClusEPhiCut->Fill(c->E(), cphi);
1122       fEmcClusEEtaCut->Fill(c->E(), ceta);
1123     }
1124   }
1125 }
1126 //________________________________________________________________________
1127 Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex)
1128 {
1129   Double_t pt = 0;
1130   if(!fTracks)
1131     return pt;
1132   if(matchIndex<0 || matchIndex>fTracks->GetEntries()){
1133     if(fDebug)
1134       printf("track-matched index out of track array range!!!\n");
1135     return pt;
1136   }
1137   AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex));
1138   if(!track){
1139     if(fDebug)
1140       printf("track-matched pointer does not exist!!!\n");
1141     return pt;
1142   }
1143   if(fESD){
1144     if(!fPrTrCuts)
1145       return pt;
1146     if(!fPrTrCuts->IsSelected(track))
1147       return pt;
1148     pt = track->Pt();
1149   }
1150   return pt;
1151 }
1152 //________________________________________________________________________
1153 void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells()
1154 {
1155   AliVCaloCells *cells = 0;
1156   cells = fESDCells;
1157   if (!cells)
1158     cells = fAODCells;
1159   if(!cells)
1160     return;
1161   Double_t maxe = 0;
1162   Double_t maxphi = -10;
1163   Int_t ncells = cells->GetNumberOfCells();
1164   Double_t eta,phi;
1165   for (Int_t i=0; i<ncells; i++) {
1166     Short_t absid = TMath::Abs(cells->GetCellNumber(i));
1167     Double_t e = cells->GetCellAmplitude(absid);
1168     if(e>0.05)
1169       fNCells50++;
1170     else
1171       continue;
1172     fGeom->EtaPhiFromIndex(absid,eta,phi);
1173     if(maxe<e){
1174       maxe = e;
1175       maxphi = phi;
1176     }
1177   }
1178   fMaxCellEPhi->Fill(maxe,maxphi);
1179
1180 }
1181 //________________________________________________________________________
1182 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
1183 {
1184   // Called once at the end of the query.
1185 }