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   fESD(0),
  75   fAOD(0),
  76   fMCEvent(0),
  77   fStack(0),
  78   fOutputList(0),
  79   fEvtSel(0),
  80   fNClusEt10(0),
  81   fRecoPV(0),
  82   fPVtxZ(0),
  83   fTrMultDist(0),
  84   fMCDirPhotonPtEtaPhi(0),
  85   fDecayPhotonPtMC(0),
  86   fCellAbsIdVsAmpl(0),
  87   fNClusHighClusE(0),
  88   fHigherPtConeM02(0),
  89   fClusEtMcPt(0),
  90   fClusMcDetaDphi(0),
  91   fNClusPerPho(0),
  92   fMcPtInConeBG(0),
  93   fMcPtInConeSBG(0),
  94   fMcPtInConeBGnoUE(0),
  95   fMcPtInConeSBGnoUE(0),
  96   fAllIsoEtMcGamma(0),
  97   fAllIsoNoUeEtMcGamma(0),
  98   fMCDirPhotonPtEtaPhiNoClus(0),
  99   fHnOutput(0)
 100 {
 101   // Default constructor.
 102   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 103 }
 104
 105 //________________________________________________________________________
 106 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
 107   AliAnalysisTaskSE(name),
 108   fESDClusters(0),
 109   fAODClusters(0),
 110   fSelPrimTracks(0),
 111   fTracks(0),
 112   fESDCells(0),
 113   fAODCells(0),
 114   fPrTrCuts(0),
 115   fGeom(0x0),
 116   fGeoName("EMCAL_COMPLETEV1"),
 117   fOADBContainer(0),
 118   fPeriod("LHC11c"),
 119   fTrigBit("kEMC7"),
 120   fIsTrain(0),
 121   fIsMc(0),
 122   fDebug(0),
 123   fPathStrOpt("/"),
 124   fExoticCut(0.97),
 125   fIsoConeR(0.4),
 126   fNDimensions(7),
 127   fECut(3.),
 128   fTrackMult(0),
 129   fMcIdFamily(""),
 130   fNClusForDirPho(0),
 131   fDirPhoPt(0),
 132   fHigherPtCone(0),
 133   fImportGeometryFromFile(0),
 134   fImportGeometryFilePath(""),
 135   fESD(0),
 136   fAOD(0),
 137   fMCEvent(0),
 138   fStack(0),
 139   fOutputList(0),
 140   fEvtSel(0),
 141   fNClusEt10(0),
 142   fRecoPV(0),
 143   fPVtxZ(0),
 144   fTrMultDist(0),
 145   fMCDirPhotonPtEtaPhi(0),
 146   fDecayPhotonPtMC(0),
 147   fCellAbsIdVsAmpl(0),
 148   fNClusHighClusE(0),
 149   fHigherPtConeM02(0),
 150   fClusEtMcPt(0),
 151   fClusMcDetaDphi(0),
 152   fNClusPerPho(0),
 153   fMcPtInConeBG(0),
 154   fMcPtInConeSBG(0),
 155   fMcPtInConeBGnoUE(0),
 156   fMcPtInConeSBGnoUE(0),
 157   fAllIsoEtMcGamma(0),
 158   fAllIsoNoUeEtMcGamma(0),
 159   fMCDirPhotonPtEtaPhiNoClus(0),
 160   fHnOutput(0)
 161 {
 162   // Constructor
 163
 164   // Define input and output slots here
 165   // Input slot #0 works with a TChain
 166   DefineInput(0, TChain::Class());
 167   // Output slot #0 id reserved by the base class for AOD
 168   // Output slot #1 writes into a TH1 container
 169   DefineOutput(1, TList::Class());
 170 }
 171
 172 //________________________________________________________________________
 173 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 174 {
 175   // Create histograms, called once.
 176
 177   fESDClusters = new TObjArray();
 178   fSelPrimTracks = new TObjArray();
 179
 180
 181   fOutputList = new TList();
 182   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 183
 184   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 185   fOADBContainer = new AliOADBContainer("AliEMCALgeo");
 186   fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo");
 187
 188   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 189   fOutputList->Add(fEvtSel);
 190
 191   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 192   fOutputList->Add(fNClusEt10);
 193
 194   fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 195   fOutputList->Add(fRecoPV);
 196
 197   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 198   fOutputList->Add(fPVtxZ);
 199
 200   fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5);
 201   fOutputList->Add(fTrMultDist);
 202
 203   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);
 204   fMCDirPhotonPtEtaPhi->Sumw2();
 205   fOutputList->Add(fMCDirPhotonPtEtaPhi);
 206
 207   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
 208   fDecayPhotonPtMC->Sumw2();
 209   fOutputList->Add(fDecayPhotonPtMC);
 210
 211   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);
 212   fOutputList->Add(fCellAbsIdVsAmpl);
 213
 214   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);
 215   fOutputList->Add(fNClusHighClusE);
 216
 217   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);
 218   fOutputList->Add(fHigherPtConeM02);
 219
 220   fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100);
 221   fOutputList->Add(fClusEtMcPt);
 222
 223   fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7);
 224   fOutputList->Add(fClusMcDetaDphi);
 225
 226   fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5);
 227   fOutputList->Add(fNClusPerPho);
 228
 229   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);
 230   fOutputList->Add(fMcPtInConeBG);
 231
 232   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);
 233   fOutputList->Add(fMcPtInConeSBG);
 234
 235   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);
 236   fOutputList->Add(fMcPtInConeBGnoUE);
 237
 238   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);
 239   fOutputList->Add(fMcPtInConeSBGnoUE);
 240
 241   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);
 242   fOutputList->Add(fAllIsoEtMcGamma);
 243
 244   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);
 245   fOutputList->Add(fAllIsoNoUeEtMcGamma);
 246
 247
 248   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);
 249   fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus);
 250
 251   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;
 252   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt};
 253   fNDimensions = sizeof(bins)/sizeof(Int_t);
 254   const Int_t ndims =   fNDimensions;
 255   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};
 256   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};
 257   if(fPeriod.Contains("11h")){
 258     xmax[12]=3999.5;
 259   }
 260   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);
 261   fOutputList->Add(fHnOutput);
 262
 263
 264
 265   PostData(1, fOutputList);
 266 }
 267
 268 //________________________________________________________________________
 269 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 270 {
 271   // Main loop, called for each event.
 272   fESDClusters = 0;
 273   fESDCells = 0;
 274   fAODClusters = 0;
 275   fAODCells = 0;
 276   // event trigger selection
 277   Bool_t isSelected = 0;
 278   if(fPeriod.Contains("11a")){
 279     if(fTrigBit.Contains("kEMC"))
 280       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 281     else
 282       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 283   }
 284   else{
 285     if(fTrigBit.Contains("kEMC"))
 286       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 287     else
 288       if(fTrigBit.Contains("kMB"))
 289         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 290       else
 291         isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 292   }
 293   if(fPeriod.Contains("11h")){
 294     if(fTrigBit.Contains("kAny"))
 295       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);
 296     if(fTrigBit.Contains("kAnyINT"))
 297       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT);
 298   }
 299   if(fIsMc)
 300     isSelected = kTRUE;
 301   if(fDebug)
 302     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 303   if(!isSelected )
 304         return;
 305   if(fIsMc){
 306     TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
 307     TFile *file = (TFile*)tree->GetCurrentFile();
 308     TString filename = file->GetName();
 309     if(!filename.Contains(fPathStrOpt.Data()))
 310       return;
 311   }
 312   AliVEvent *event = (AliVEvent*)InputEvent();
 313   if (!event) {
 314     printf("ERROR: event not available\n");
 315     return;
 316   }
 317   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 318   if(fDebug)
 319     printf("run number = %d\n",runnumber);
 320
 321   fESD = dynamic_cast<AliESDEvent*>(event);
 322   if(!fESD){
 323     fAOD = dynamic_cast<AliAODEvent*>(event);
 324     if(!fAOD){
 325       printf("ERROR: Invalid type of event!!!\n");
 326       return;
 327     }
 328     else if(fDebug)
 329       printf("AOD EVENT!\n");
 330   }
 331
 332   fEvtSel->Fill(0);
 333   if(fDebug)
 334     printf("event is ok,\n run number=%d\n",runnumber);
 335
 336
 337   AliVVertex *pv = (AliVVertex*)event->GetPrimaryVertex();
 338   Bool_t pvStatus = kTRUE;
 339   if(fESD){
 340     AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex();
 341     pvStatus = esdv->GetStatus();
 342   }
 343   if(!pv)
 344     return;
 345   if(!pvStatus)
 346     fRecoPV->Fill(0);
 347   else
 348     fRecoPV->Fill(1);
 349   fPVtxZ->Fill(pv->GetZ());
 350   if(TMath::Abs(pv->GetZ())>15)
 351     return;
 352   if(fDebug)
 353     printf("passed vertex cut\n");
 354
 355   fEvtSel->Fill(1);
 356   if(fESD)
 357     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 358   if(fAOD)
 359     fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks());
 360
 361   if(!fTracks){
 362     AliError("Track array in event is NULL!");
 363     if(fDebug)
 364       printf("returning due to not finding tracks!\n");
 365     return;
 366   }
 367   // Track loop to fill a pT spectrum
 368   const Int_t Ntracks = fTracks->GetEntriesFast();
 369   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 370     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 371     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 372     AliVTrack *track = (AliVTrack*)fTracks->At(iTracks);
 373     if (!track)
 374       continue;
 375     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 376     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
 377     if (esdTrack && fPrTrCuts && fPrTrCuts->IsSelected(track)){
 378       fSelPrimTracks->Add(track);
 379       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 380     }
 381     else if(aodTrack)
 382       fSelPrimTracks->Add(track);
 383   }
 384
 385   TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
 386   for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 387     if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 388       break;
 389     /*if(fESD)
 390       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 391       else*/
 392     // if(event->GetEMCALMatrix(mod))
 393     fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
 394     fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
 395   }
 396
 397   if(fESD){
 398     AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 399     fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 400     if(fDebug)
 401       printf("ESD Track mult= %d\n",fTrackMult);
 402   }
 403   else if(fAOD){
 404     fTrackMult = Ntracks;
 405     if(fDebug)
 406       printf("AOD Track mult= %d\n",fTrackMult);
 407   }
 408   fTrMultDist->Fill(fTrackMult);
 409
 410   if(fESD){
 411     TList *l = fESD->GetList();
 412     fESDClusters =  dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
 413     fESDCells = fESD->GetEMCALCells();
 414     if(fDebug)
 415       printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast());
 416   }
 417   else if(fAOD){
 418     fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
 419     fAODCells = fAOD->GetEMCALCells();
 420     if(fDebug)
 421       printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast());
 422   }
 423
 424
 425   fMCEvent = MCEvent();
 426   if(fMCEvent){
 427     if(fDebug)
 428       std::cout<<"MCevent exists\n";
 429     fStack = (AliStack*)fMCEvent->Stack();
 430   }
 431   else{
 432     if(fDebug)
 433       std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
 434   }
 435   FollowGamma();
 436   if(fDebug)
 437     printf("passed calling of FollowGamma\n");
 438   FillClusHists();
 439   if(fDebug)
 440     printf("passed calling of FillClusHists\n");
 441   FillMcHists();
 442   if(fDebug)
 443     printf("passed calling of FillMcHists\n");
 444   /*if(fESD)
 445     fESDClusters->Clear();*/
 446   fSelPrimTracks->Clear();
 447   fNClusForDirPho = 0;
 448
 449   PostData(1, fOutputList);
 450 }
 451
 452 //________________________________________________________________________
 453 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 454 {
 455   if(fDebug)
 456     printf("Inside FillClusHists()....\n");
 457   // Fill cluster histograms.
 458   TObjArray *clusters = fESDClusters;
 459
 460   if (!clusters){
 461     clusters = fAODClusters;
 462     if(fDebug)
 463       printf("ESD clusters empty...");
 464   }
 465   if (!clusters){
 466     if(fDebug)
 467       printf("  and AOD clusters as well!!!\n");
 468     return;
 469   }
 470   if(fDebug)
 471     printf("\n");
 472
 473   const Int_t nclus = clusters->GetEntries();
 474   if(nclus==0)
 475     return;
 476   if(fDebug)
 477     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 478   Double_t maxE = 0;
 479   Int_t nclus10 = 0;
 480   Double_t ptmc=-1;
 481   for(Int_t ic=0;ic<nclus;ic++){
 482     maxE=0;
 483     AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic));
 484     if(!c)
 485       continue;
 486     if(!c->IsEMCAL())
 487       continue;
 488     if(c->E()<fECut)
 489       continue;
 490     Short_t id;
 491     Double_t Emax = GetMaxCellEnergy( c, id);
 492     Double_t Ecross = GetCrossEnergy( c, id);
 493     if((1-Ecross/Emax)>fExoticCut)
 494       continue;
 495     Float_t clsPos[3] = {0,0,0};
 496     c->GetPosition(clsPos);
 497     TVector3 clsVec(clsPos);
 498     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 499     if(fDebug)
 500       printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E());
 501     if(Et>10)
 502       nclus10++;
 503     Float_t ceiso, cephiband, cecore;
 504     Float_t triso, trphiband, trcore;
 505     Float_t alliso, allphiband, allcore;
 506     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 507     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 508     GetCeIso(clsVec, id, ceiso, cephiband, cecore);
 509     GetTrIso(clsVec, triso, trphiband, trcore);
 510     Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz());
 511     if(Et>10 && Et<15 && dr>0.025){
 512       fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02());
 513       if(fDebug)
 514         printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02());
 515     }
 516     alliso = ceiso + triso;
 517     allphiband = cephiband + trphiband;
 518     allcore = cecore + trcore;
 519     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 520     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 521     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 522     Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR);
 523     if(fDebug && Et>10)
 524       printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et);
 525     if(c->GetM02()>0.5 && c->GetM02()<2.0){
 526       fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum);
 527       fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum);
 528     }
 529     if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){
 530       fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum);
 531       fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum);
 532       if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){
 533         fAllIsoEtMcGamma->Fill(Et, alliso-/*Et*/cecore-allisoue);
 534         fAllIsoNoUeEtMcGamma->Fill(Et, alliso-/*Et*/cecore);
 535       }
 536     }
 537     const Int_t ndims =   fNDimensions;
 538     Double_t outputValues[ndims];
 539     ptmc = GetClusSource(c);
 540     outputValues[0] = Et;
 541     outputValues[1] = c->GetM02();
 542     outputValues[2] = ceiso-Et/*cecore*/-ceisoue;
 543     outputValues[3] = triso-trisoue;
 544     outputValues[4] = alliso-Et/*cecore*/-allisoue - trcore;
 545     outputValues[5] = ceiso-Et;
 546     outputValues[6] = alliso-Et - trcore;
 547     outputValues[7] = c->GetTrackDx();
 548     outputValues[8] = c->GetTrackDz();
 549     outputValues[9] = clsVec.Eta();
 550     outputValues[10] = clsVec.Phi();
 551     if(fESDCells)
 552       outputValues[11] = fESDCells->GetCellTime(id);
 553     else if(fAODCells)
 554       outputValues[11] = fAODCells->GetCellTime(id);
 555     outputValues[12] = fTrackMult;
 556     outputValues[13] = ptmc;
 557     fHnOutput->Fill(outputValues);
 558     if(c->E()>maxE)
 559       maxE = c->E();
 560   }
 561   fNClusHighClusE->Fill(maxE,nclus);
 562   fNClusEt10->Fill(nclus10);
 563   fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho);
 564 }
 565
 566 //________________________________________________________________________
 567 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core)
 568 {
 569   // Get cell isolation.
 570   AliVCaloCells *cells = fESDCells;
 571   if (!cells)
 572     cells = fAODCells;
 573   if (!cells)
 574     return;
 575
 576   const Int_t ncells = cells->GetNumberOfCells();
 577   Float_t totiso=0;
 578   Float_t totband=0;
 579   Float_t totcore=0;
 580   Float_t etacl = vec.Eta();
 581   Float_t phicl = vec.Phi();
 582   Float_t thetacl = vec.Theta();
 583   Float_t maxtcl = cells->GetCellTime(maxid);
 584   if(phicl<0)
 585     phicl+=TMath::TwoPi();
 586   Int_t absid = -1;
 587   Float_t eta=-1, phi=-1;
 588   for(int icell=0;icell<ncells;icell++){
 589     absid = TMath::Abs(cells->GetCellNumber(icell));
 590     Float_t celltime = cells->GetCellTime(absid);
 591     //if(TMath::Abs(celltime)>2e-8 && (!fIsMc))
 592     if(TMath::Abs(celltime-maxtcl)>2e-8 )
 593       continue;
 594     if(!fGeom)
 595       return;
 596     fGeom->EtaPhiFromIndex(absid,eta,phi);
 597     Float_t dphi = TMath::Abs(phi-phicl);
 598     Float_t deta = TMath::Abs(eta-etacl);
 599     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 600     Float_t etcell = cells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 601     if(R<fIsoConeR){
 602       totiso += etcell;
 603       if(R<0.04)
 604         totcore += etcell;
 605     }
 606     else{
 607       if(dphi>fIsoConeR)
 608         continue;
 609       if(deta<fIsoConeR)
 610         continue;
 611       totband += cells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 612     }
 613   }
 614   iso = totiso;
 615   phiband = totband;
 616   core = totcore;
 617 }
 618
 619 //________________________________________________________________________
 620 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 621 {
 622   // Get track isolation.
 623
 624   if(!fSelPrimTracks)
 625     return;
 626   fHigherPtCone = 0;
 627   const Int_t ntracks = fSelPrimTracks->GetEntries();
 628   Double_t totiso=0;
 629   Double_t totband=0;
 630   Double_t totcore=0;
 631   Double_t etacl = vec.Eta();
 632   Double_t phicl = vec.Phi();
 633   if(phicl<0)
 634     phicl+=TMath::TwoPi();
 635   for(int itrack=0;itrack<ntracks;itrack++){
 636     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 637     if(!track)
 638       continue;
 639     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 640     Double_t deta = TMath::Abs(track->Eta()-etacl);
 641     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 642     Double_t pt = track->Pt();
 643     if(pt>fHigherPtCone)
 644       fHigherPtCone = pt;
 645     if(R<fIsoConeR){
 646       totiso += track->Pt();
 647       if(R<0.04)
 648         totcore += pt;
 649     }
 650     else{
 651       if(dphi>fIsoConeR)
 652         continue;
 653       if(deta<fIsoConeR)
 654         continue;
 655       totband += track->Pt();
 656     }
 657   }
 658   iso = totiso;
 659   phiband = totband;
 660   core = totcore;
 661 }
 662
 663 //________________________________________________________________________
 664 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 665 {
 666   // Calculate the energy of cross cells around the leading cell.
 667
 668   AliVCaloCells *cells = 0;
 669   cells = fESDCells;
 670   if (!cells)
 671     cells = fAODCells;
 672   if (!cells)
 673     return 0;
 674
 675   if (!fGeom)
 676     return 0;
 677
 678   Int_t iSupMod = -1;
 679   Int_t iTower  = -1;
 680   Int_t iIphi   = -1;
 681   Int_t iIeta   = -1;
 682   Int_t iphi    = -1;
 683   Int_t ieta    = -1;
 684   Int_t iphis   = -1;
 685   Int_t ietas   = -1;
 686
 687   Double_t crossEnergy = 0;
 688
 689   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 690   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 691
 692   Int_t ncells = cluster->GetNCells();
 693   for (Int_t i=0; i<ncells; i++) {
 694     Int_t cellAbsId = cluster->GetCellAbsId(i);
 695     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 696     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 697     Int_t aphidiff = TMath::Abs(iphi-iphis);
 698     if (aphidiff>1)
 699       continue;
 700     Int_t aetadiff = TMath::Abs(ieta-ietas);
 701     if (aetadiff>1)
 702       continue;
 703     if ( (aphidiff==1 && aetadiff==0) ||
 704         (aphidiff==0 && aetadiff==1) ) {
 705       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 706     }
 707   }
 708
 709   return crossEnergy;
 710 }
 711
 712 //________________________________________________________________________
 713 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 714 {
 715   // Get maximum energy of attached cell.
 716
 717   id = -1;
 718
 719   AliVCaloCells *cells = 0;
 720   cells = fESDCells;
 721   if (!cells)
 722     cells = fAODCells;
 723   if(!cells)
 724     return 0;
 725
 726   Double_t maxe = 0;
 727   Int_t ncells = cluster->GetNCells();
 728   for (Int_t i=0; i<ncells; i++) {
 729     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 730     if (e>maxe) {
 731       maxe = e;
 732       id   = cluster->GetCellAbsId(i);
 733     }
 734   }
 735   return maxe;
 736 }
 737
 738 //________________________________________________________________________
 739 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
 740 {
 741   if(!fStack)
 742     return;
 743   Int_t nTracks = fStack->GetNtrack();
 744   if(fDebug)
 745     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
 746   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 747     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 748     if(!mcp)
 749       continue;
 750     Int_t pdg = mcp->GetPdgCode();
 751     if(pdg!=22)
 752       continue;
 753     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
 754       continue;
 755     Int_t imom = mcp->GetMother(0);
 756     if(imom<0 || imom>nTracks)
 757       continue;
 758     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
 759     if(!mcmom)
 760       continue;
 761     Int_t pdgMom = mcmom->GetPdgCode();
 762     if((imom==6 || imom==7) && pdgMom==22) {
 763       fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 764       if(fNClusForDirPho==0)
 765         fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
 766       if(fDebug){
 767         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());
 768         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());
 769       }
 770     }
 771     else{
 772       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
 773         fDecayPhotonPtMC->Fill(mcp->Pt());
 774     }
 775   }
 776 }
 777 //________________________________________________________________________
 778 Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c)
 779 {
 780   if(!c)
 781     return -1;
 782   if(!fStack)
 783     return -1;
 784   Int_t nmcp = fStack->GetNtrack();
 785   Int_t clabel = c->GetLabel();
 786   if(fDebug && fMcIdFamily.Contains(Form("%d",clabel)))
 787     printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel);
 788   if(!fMcIdFamily.Contains(Form("%d",clabel)))
 789     return -1;
 790   fNClusForDirPho++;
 791   TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1);
 792   Int_t partonpos = partonposstr.Atoi();
 793   if(fDebug)
 794     printf("\t^^^^ parton position = %d ^^^^\n",partonpos);
 795   if(clabel<0 || clabel>nmcp)
 796     return 0;
 797   Float_t clsPos[3] = {0,0,0};
 798   c->GetPosition(clsPos);
 799   TVector3 clsVec(clsPos);
 800   Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 801   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos));
 802   if(!mcp)
 803     return -1;
 804   if(fDebug){
 805     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);
 806   }
 807   Int_t lab1 =  mcp->GetFirstDaughter();
 808   if(lab1<0 || lab1>nmcp)
 809     return -1;
 810   TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1));
 811   if(!mcd)
 812     return -1;
 813   if(fDebug)
 814     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);
 815   if(mcd->GetPdgCode()==22){
 816     fClusEtMcPt->Fill(mcd->Pt(), Et);
 817     fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi());
 818   }
 819   else{
 820     printf("Warning: daughter of photon parton is not a photon\n");
 821     return -1;
 822   }
 823   return fDirPhoPt;
 824 }
 825 //________________________________________________________________________
 826 void AliAnalysisTaskEMCALIsoPhoton::FollowGamma()
 827 {
 828   if(!fStack)
 829     return;
 830   Int_t selfid = 6;
 831   TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid));
 832   if(!mcp)
 833     return;
 834   if(mcp->GetPdgCode()!=22){
 835     mcp = static_cast<TParticle*>(fStack->Particle(++selfid));
 836     if(!mcp)
 837       return;
 838   }
 839   Int_t daug0f =  mcp->GetFirstDaughter();
 840   Int_t daug0l =  mcp->GetLastDaughter();
 841   Int_t nd0 = daug0l - daug0f;
 842   if(fDebug)
 843     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);
 844   fMcIdFamily = Form("%d,",selfid);
 845   GetDaughtersInfo(daug0f,daug0l, selfid,"");
 846   if(fDebug){
 847     printf("\t---- end of the prompt  gamma's family tree ----\n\n");
 848     printf("Family id string = %s,\n\n",fMcIdFamily.Data());
 849   }
 850   TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f));
 851   if(!md)
 852     return;
 853   fDirPhoPt = md->Pt();
 854 }
 855 //________________________________________________________________________
 856 void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind)
 857 {
 858   int nmcp = fStack->GetNtrack();
 859   if(firstd<0 || firstd>nmcp)
 860     return;
 861   if(lastd<0 || lastd>nmcp)
 862     return;
 863   if(firstd>lastd){
 864     printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd);
 865     return;
 866   }
 867   TString indenter = Form("\t%s",inputind);
 868   TParticle *dp = 0x0;
 869   if(fDebug)
 870     printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid);
 871   for(int id=firstd; id<lastd+1; id++){
 872     dp =  static_cast<TParticle*>(fStack->Particle(id));
 873     if(!dp)
 874       continue;
 875     Int_t dfd = dp->GetFirstDaughter();
 876     Int_t dld = dp->GetLastDaughter();
 877     Int_t dnd =  dld - dfd + 1;
 878     Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy());
 879     if(fDebug)
 880       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);
 881     fMcIdFamily += Form("%d,",id);
 882     GetDaughtersInfo(dfd,dld,id,indenter.Data());
 883   }
 884 }
 885
 886 //________________________________________________________________________
 887 Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R)
 888 {
 889   if(!fStack)
 890     return 0;
 891   if(fDebug)
 892     printf("Inside GetMcPtSumInCone!!\n");
 893   Int_t nTracks = fStack->GetNtrack();
 894   Float_t ptsum = 0;
 895   TString addpartlabels = "";
 896   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 897     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 898     if(!mcp)
 899       continue;
 900     Int_t firstd = mcp->GetFirstDaughter();
 901     Int_t lastd = mcp->GetLastDaughter();
 902     if((iTrack==6 || iTrack==7) && mcp->GetPdgCode()==22){
 903      for(Int_t id=firstd;id<=lastd;id++)
 904       addpartlabels += Form("%d,",id);
 905      continue;
 906     }
 907    if(mcp->Rho()>2.5)
 908       continue;
 909     else {
 910       if(fDebug)
 911         printf("      >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz());
 912     }
 913     Float_t dphi = mcp->Phi() - phiclus;
 914     Float_t deta = mcp->Eta() - etaclus;
 915     if(fDebug)
 916       printf("      >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta);
 917     if(deta>10)
 918       continue;
 919     Float_t dR = TMath::Sqrt(dphi*dphi +  deta*deta);
 920     if(dR>R)
 921       continue;
 922     if(addpartlabels.Contains(Form("%d",iTrack))){
 923       if(fDebug)
 924         printf("      >>>> list of particles included (and their daughters) %s\n",addpartlabels.Data());
 925       continue;
 926     }
 927     addpartlabels += Form("%d,",iTrack);
 928     for(Int_t id=firstd;id<=lastd;id++)
 929       addpartlabels += Form("%d,",id);
 930     ptsum += mcp->Pt();
 931   }
 932   return ptsum;
 933 }
 934 //________________________________________________________________________
 935 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 936 {
 937   // Called once at the end of the query.
 938 }