PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

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