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 <THnSparse.h>
  11 #include <TLorentzVector.h>
  12 #include <TList.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliAnalysisTask.h"
  16 #include "AliEMCALGeometry.h"
  17 #include "AliEMCALRecoUtils.h"
  18 #include "AliESDCaloCells.h"
  19 #include "AliESDCaloCluster.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDHeader.h"
  22 #include "AliESDInputHandler.h"
  23 #include "AliESDUtils.h"
  24 #include "AliESDtrack.h"
  25 #include "AliESDtrackCuts.h"
  26 #include "AliMCEvent.h"
  27 #include "AliMCEventHandler.h"
  28 #include "AliStack.h"
  29 #include "AliV0vertexer.h"
  30 #include "AliVCluster.h"
  31
  32 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  36   AliAnalysisTaskSE(),
  37   fCaloClusters(0),
  38   fSelPrimTracks(0),
  39   fTracks(0),
  40   fEMCalCells(0),
  41   fPrTrCuts(0),
  42   fGeom(0x0),
  43   fGeoName("EMCAL_COMPLETEV1"),
  44   fPeriod("LHC11c"),
  45   fTrigBit("kEMC7"),
  46   fIsTrain(0),
  47   fIsMc(0),
  48   fDebug(0),
  49   fExoticCut(0.97),
  50   fIsoConeR(0.4),
  51   fNDimensions(7),
  52   fECut(3.),
  53   fTrackMult(0),
  54   fESD(0),
  55   fMCEvent(0),
  56   fStack(0),
  57   fOutputList(0),
  58   fEvtSel(0),
  59   fNClusEt10(0),
  60   fPVtxZ(0),
  61   fDirPhotonPtMC(0),
  62   fDecayPhotonPtMC(0),
  63   fCellAbsIdVsAmpl(0),
  64   fNClusHighClusE(0),
  65   fHnOutput(0)
  66 {
  67   // Default constructor.
  68 }
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
  72   AliAnalysisTaskSE(name),
  73   fCaloClusters(0),
  74   fSelPrimTracks(0),
  75   fTracks(0),
  76   fEMCalCells(0),
  77   fPrTrCuts(0),
  78   fGeom(0x0),
  79   fGeoName("EMCAL_COMPLETEV1"),
  80   fPeriod("LHC11c"),
  81   fTrigBit("kEMC7"),
  82   fIsTrain(0),
  83   fIsMc(0),
  84   fDebug(0),
  85   fExoticCut(0.97),
  86   fIsoConeR(0.4),
  87   fNDimensions(7),
  88   fECut(3.),
  89   fTrackMult(0),
  90   fESD(0),
  91   fMCEvent(0),
  92   fStack(0),
  93   fOutputList(0),
  94   fEvtSel(0),
  95   fNClusEt10(0),
  96   fPVtxZ(0),
  97   fDirPhotonPtMC(0),
  98   fDecayPhotonPtMC(0),
  99   fCellAbsIdVsAmpl(0),
 100   fNClusHighClusE(0),
 101   fHnOutput(0)
 102 {
 103   // Constructor
 104
 105   // Define input and output slots here
 106   // Input slot #0 works with a TChain
 107   DefineInput(0, TChain::Class());
 108   // Output slot #0 id reserved by the base class for AOD
 109   // Output slot #1 writes into a TH1 container
 110   DefineOutput(1, TList::Class());
 111 }
 112
 113 //________________________________________________________________________
 114 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 115 {
 116   // Create histograms, called once.
 117
 118   fCaloClusters = new TRefArray();
 119   fSelPrimTracks = new TObjArray();
 120
 121
 122   fOutputList = new TList();
 123   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 124
 125   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 126
 127   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 128   fOutputList->Add(fEvtSel);
 129
 130   fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
 131   fOutputList->Add(fNClusEt10);
 132
 133   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 134   fOutputList->Add(fPVtxZ);
 135
 136   fDirPhotonPtMC = new TH1F("hDirPhotonPtMC","photon (gq->#gammaq) p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
 137   fDirPhotonPtMC->Sumw2();
 138   fOutputList->Add(fDirPhotonPtMC);
 139
 140   fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
 141   fDecayPhotonPtMC->Sumw2();
 142   fOutputList->Add(fDecayPhotonPtMC);
 143
 144   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);
 145   fOutputList->Add(fCellAbsIdVsAmpl);
 146
 147   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);
 148   fOutputList->Add(fNClusHighClusE);
 149
 150   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=20;
 151   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult};
 152   fNDimensions = sizeof(bins)/sizeof(Int_t);
 153   const Int_t ndims =   fNDimensions;
 154   Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5};
 155   Double_t xmax[] = { 100., 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5};
 156   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", ndims, bins, xmin, xmax);
 157   fOutputList->Add(fHnOutput);
 158
 159
 160
 161   PostData(1, fOutputList);
 162 }
 163
 164 //________________________________________________________________________
 165 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 166 {
 167   // Main loop, called for each event.
 168
 169   // event trigger selection
 170   Bool_t isSelected = 0;
 171   if(fPeriod.Contains("11a")){
 172     if(fTrigBit.Contains("kEMC"))
 173       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 174     else
 175       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 176   }
 177   else{
 178     if(fTrigBit.Contains("kEMC"))
 179       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 180     else
 181       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 182   }
 183   if(fIsMc)
 184     isSelected = kTRUE;
 185   if(fDebug)
 186     printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
 187   if(!isSelected )
 188         return;
 189
 190   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 191   if (!fESD) {
 192     printf("ERROR: fESD not available\n");
 193     return;
 194   }
 195
 196   fEvtSel->Fill(0);
 197   if(fDebug)
 198     printf("fESD is ok\n");
 199
 200   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 201   if(!pv)
 202     return;
 203   fPVtxZ->Fill(pv->GetZ());
 204   if(TMath::Abs(pv->GetZ())>15)
 205     return;
 206   if(fDebug)
 207     printf("passed vertex cut\n");
 208
 209   fEvtSel->Fill(1);
 210
 211   if (!fTracks)
 212     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 213   // Track loop to fill a pT spectrum
 214   const Int_t Ntracks = fTracks->GetEntriesFast();
 215   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 216     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 217     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 218     AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks));
 219     if (!track)
 220       continue;
 221     if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 222       fSelPrimTracks->Add(track);
 223       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 224     }
 225   }
 226
 227   if(!fIsTrain){
 228     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 229       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 230         break;
 231       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 232     }
 233   }
 234   AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
 235   fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
 236
 237   fESD->GetEMCALClusters(fCaloClusters);
 238   fEMCalCells = fESD->GetEMCALCells();
 239
 240
 241   FillClusHists();
 242   if(fDebug)
 243     printf("passed calling of FillClusHists\n");
 244
 245   fCaloClusters->Clear();
 246   fSelPrimTracks->Clear();
 247
 248   fMCEvent = MCEvent();
 249   if(fMCEvent){
 250     if(fDebug)
 251       cout<<"MCevent exists\n";
 252     fStack = (AliStack*)fMCEvent->Stack();
 253   }
 254   else{
 255     if(fDebug)
 256       cout<<"ERROR: NO MC EVENT!!!!!!\n";
 257   }
 258   FillMcHists();
 259   if(fDebug)
 260     printf("passed calling of FillMcHists\n");
 261
 262
 263   PostData(1, fOutputList);
 264 }
 265
 266 //________________________________________________________________________
 267 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 268 {
 269   // Fill cluster histograms.
 270
 271   if(!fCaloClusters)
 272     return;
 273   const Int_t nclus = fCaloClusters->GetEntries();
 274   if(nclus==0)
 275     return;
 276   if(fDebug)
 277     printf("Inside FillClusHists and there are %d clusters\n",nclus);
 278   Double_t maxE = 0;
 279   Int_t nclus10 = 0;
 280   for(Int_t ic=0;ic<nclus;ic++){
 281     maxE=0;
 282     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 283     if(!c)
 284       continue;
 285     if(!c->IsEMCAL())
 286       continue;
 287     if(c->E()<fECut)
 288       continue;
 289     Short_t id;
 290     Double_t Emax = GetMaxCellEnergy( c, id);
 291     Double_t Ecross = GetCrossEnergy( c, id);
 292     if((1-Ecross/Emax)>fExoticCut)
 293       continue;
 294     Float_t clsPos[3] = {0,0,0};
 295     c->GetPosition(clsPos);
 296     TVector3 clsVec(clsPos);
 297     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 298     if(Et>10)
 299       nclus10++;
 300     Float_t ceiso, cephiband, cecore;
 301     Float_t triso, trphiband, trcore;
 302     Float_t alliso, allphiband, allcore;
 303     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 304     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 305     GetCeIso(clsVec, ceiso, cephiband, cecore);
 306     GetTrIso(clsVec, triso, trphiband, trcore);
 307     alliso = ceiso + triso;
 308     allphiband = cephiband + trphiband;
 309     allcore = cecore + trcore;
 310     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 311     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 312     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 313     const Int_t ndims =   fNDimensions;
 314     Double_t outputValues[ndims];
 315     outputValues[0] = Et;
 316     outputValues[1] = c->GetM02();
 317     outputValues[2] = ceiso-cecore-ceisoue;
 318     outputValues[3] = triso-trisoue;
 319     outputValues[4] = alliso-cecore-allisoue;
 320     outputValues[5] = ceiso-Et;
 321     outputValues[6] = alliso-Et;
 322     outputValues[7] = c->GetTrackDx();
 323     outputValues[8] = c->GetTrackDz();
 324     outputValues[9] = clsVec.Eta();
 325     outputValues[10] = clsVec.Phi();
 326     outputValues[11] = fEMCalCells->GetCellTime(id);
 327     outputValues[12] = fTrackMult;
 328     fHnOutput->Fill(outputValues);
 329     if(c->E()>maxE)
 330       maxE = c->E();
 331   }
 332   fNClusHighClusE->Fill(maxE,nclus);
 333   fNClusEt10->Fill(nclus10);
 334 }
 335
 336 //________________________________________________________________________
 337 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 338 {
 339   // Get cell isolation.
 340
 341   if(!fEMCalCells)
 342     return;
 343   const Int_t ncells = fEMCalCells->GetNumberOfCells();
 344   Float_t totiso=0;
 345   Float_t totband=0;
 346   Float_t totcore=0;
 347   Float_t etacl = vec.Eta();
 348   Float_t phicl = vec.Phi();
 349   Float_t thetacl = vec.Theta();
 350   if(phicl<0)
 351     phicl+=TMath::TwoPi();
 352   Int_t absid = -1;
 353   Float_t eta=-1, phi=-1;
 354   for(int icell=0;icell<ncells;icell++){
 355     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 356     if(!fGeom)
 357       return;
 358     fGeom->EtaPhiFromIndex(absid,eta,phi);
 359     Float_t dphi = TMath::Abs(phi-phicl);
 360     Float_t deta = TMath::Abs(eta-etacl);
 361     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 362     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 363     if(R<fIsoConeR){
 364       totiso += etcell;
 365       if(R<0.04)
 366         totcore += etcell;
 367     }
 368     else{
 369       if(dphi>fIsoConeR)
 370         continue;
 371       if(deta<fIsoConeR)
 372         continue;
 373       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 374     }
 375   }
 376   iso = totiso;
 377   phiband = totband;
 378   core = totcore;
 379 }
 380
 381 //________________________________________________________________________
 382 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 383 {
 384   // Get track isolation.
 385
 386   if(!fSelPrimTracks)
 387     return;
 388   const Int_t ntracks = fSelPrimTracks->GetEntries();
 389   Double_t totiso=0;
 390   Double_t totband=0;
 391   Double_t totcore=0;
 392   Double_t etacl = vec.Eta();
 393   Double_t phicl = vec.Phi();
 394   if(phicl<0)
 395     phicl+=TMath::TwoPi();
 396   for(int itrack=0;itrack<ntracks;itrack++){
 397     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 398     if(!track)
 399       continue;
 400     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 401     Double_t deta = TMath::Abs(track->Eta()-etacl);
 402     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 403     Double_t pt = track->Pt();
 404     if(R<fIsoConeR){
 405       totiso += track->Pt();
 406       if(R<0.04)
 407         totcore += pt;
 408     }
 409     else{
 410       if(dphi>fIsoConeR)
 411         continue;
 412       if(deta<fIsoConeR)
 413         continue;
 414       totband += track->Pt();
 415     }
 416   }
 417   iso = totiso;
 418   phiband = totband;
 419   core = totcore;
 420 }
 421
 422 //________________________________________________________________________
 423 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 424 {
 425   // Calculate the energy of cross cells around the leading cell.
 426
 427   AliVCaloCells *cells = 0;
 428   cells = fESD->GetEMCALCells();
 429   if (!cells)
 430     return 0;
 431
 432   if (!fGeom)
 433     return 0;
 434
 435   Int_t iSupMod = -1;
 436   Int_t iTower  = -1;
 437   Int_t iIphi   = -1;
 438   Int_t iIeta   = -1;
 439   Int_t iphi    = -1;
 440   Int_t ieta    = -1;
 441   Int_t iphis   = -1;
 442   Int_t ietas   = -1;
 443
 444   Double_t crossEnergy = 0;
 445
 446   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 447   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 448
 449   Int_t ncells = cluster->GetNCells();
 450   for (Int_t i=0; i<ncells; i++) {
 451     Int_t cellAbsId = cluster->GetCellAbsId(i);
 452     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 453     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 454     Int_t aphidiff = TMath::Abs(iphi-iphis);
 455     if (aphidiff>1)
 456       continue;
 457     Int_t aetadiff = TMath::Abs(ieta-ietas);
 458     if (aetadiff>1)
 459       continue;
 460     if ( (aphidiff==1 && aetadiff==0) ||
 461         (aphidiff==0 && aetadiff==1) ) {
 462       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 463     }
 464   }
 465
 466   return crossEnergy;
 467 }
 468
 469 //________________________________________________________________________
 470 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 471 {
 472   // Get maximum energy of attached cell.
 473
 474   id = -1;
 475
 476   AliVCaloCells *cells = 0;
 477   cells = fESD->GetEMCALCells();
 478   if (!cells)
 479     return 0;
 480
 481   Double_t maxe = 0;
 482   Int_t ncells = cluster->GetNCells();
 483   for (Int_t i=0; i<ncells; i++) {
 484     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 485     if (e>maxe) {
 486       maxe = e;
 487       id   = cluster->GetCellAbsId(i);
 488     }
 489   }
 490   return maxe;
 491 }
 492
 493 //________________________________________________________________________
 494 void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
 495 {
 496   if(!fStack)
 497     return;
 498   Int_t nTracks = fStack->GetNtrack();
 499   if(fDebug)
 500     printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
 501   for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
 502     TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));
 503     if(!mcp)
 504       continue;
 505     Int_t pdg = mcp->GetPdgCode();
 506     if(pdg!=22)
 507       continue;
 508     if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
 509       continue;
 510     Int_t imom = mcp->GetMother(0);
 511     if(imom<0 || imom>nTracks)
 512       continue;
 513     TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));
 514     if(!mcmom)
 515       continue;
 516     Int_t pdgMom = mcmom->GetPdgCode();
 517     if(imom==6 || imom==7 && pdgMom==22) {
 518       fDirPhotonPtMC->Fill(mcp->Pt());
 519       if(fDebug)
 520         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());
 521       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());
 522     }
 523     else{
 524       if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
 525         fDecayPhotonPtMC->Fill(mcp->Pt());
 526     }
 527   }
 528 }
 529 //________________________________________________________________________
 530 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 531 {
 532   // Called once at the end of the query.
 533 }