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 "AliESDpid.h"
  25 #include "AliESDtrack.h"
  26 #include "AliESDtrackCuts.h"
  27 #include "AliESDv0.h"
  28 #include "AliKFParticle.h"
  29 #include "AliMCEvent.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliStack.h"
  32 #include "AliV0vertexer.h"
  33 #include "AliVCluster.h"
  34
  35 ClassImp(AliAnalysisTaskEMCALIsoPhoton)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
  39   AliAnalysisTaskSE(),
  40   fCaloClusters(0),
  41   fSelPrimTracks(0),
  42   fTracks(0),
  43   fEMCalCells(0),
  44   fPrTrCuts(0),
  45   fGeom(0x0),
  46   fGeoName("EMCAL_COMPLETEV1"),
  47   fPeriod("LHC11c"),
  48   fTrigBit("kEMC7"),
  49   fIsTrain(0),
  50   fExoticCut(0.97),
  51   fIsoConeR(0.4),
  52   fNDimensions(7),
  53   fECut(3.),
  54   fESD(0),
  55   fOutputList(0),
  56   fEvtSel(0),
  57   fPVtxZ(0),
  58   fCellAbsIdVsAmpl(0),
  59   fNClusHighClusE(0),
  60   fHnOutput(0)
  61 {
  62   // Default constructor.
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
  67   AliAnalysisTaskSE(name),
  68   fCaloClusters(0),
  69   fSelPrimTracks(0),
  70   fTracks(0),
  71   fEMCalCells(0),
  72   fPrTrCuts(0),
  73   fGeom(0x0),
  74   fGeoName("EMCAL_COMPLETEV1"),
  75   fPeriod("LHC11c"),
  76   fTrigBit("kEMC7"),
  77   fIsTrain(0),
  78   fExoticCut(0.97),
  79   fIsoConeR(0.4),
  80   fNDimensions(7),
  81   fECut(3.),
  82   fESD(0),
  83   fOutputList(0),
  84   fEvtSel(0),
  85   fPVtxZ(0),
  86   fCellAbsIdVsAmpl(0),
  87   fNClusHighClusE(0),
  88   fHnOutput(0)
  89 {
  90   // Constructor
  91
  92   // Define input and output slots here
  93   // Input slot #0 works with a TChain
  94   DefineInput(0, TChain::Class());
  95   // Output slot #0 id reserved by the base class for AOD
  96   // Output slot #1 writes into a TH1 container
  97   DefineOutput(1, TList::Class());
  98 }
  99
 100 //________________________________________________________________________
 101 void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
 102 {
 103   // Create histograms, called once.
 104
 105   fCaloClusters = new TRefArray();
 106   fSelPrimTracks = new TObjArray();
 107
 108
 109   fOutputList = new TList();
 110   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 111
 112   fGeom = AliEMCALGeometry::GetInstance(fGeoName);
 113
 114   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 115   fOutputList->Add(fEvtSel);
 116
 117
 118   fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
 119   fOutputList->Add(fPVtxZ);
 120
 121   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);
 122   fOutputList->Add(fCellAbsIdVsAmpl);
 123
 124   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);
 125   fOutputList->Add(fNClusHighClusE);
 126
 127   const Int_t ndims =   fNDimensions;
 128   Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000, nTrClDphi=200, nTrClDeta=100;
 129   Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nTrClDphi, nTrClDeta};
 130   Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -0.1,-0.05};
 131   Double_t xmax[] = { 100., 4., 190., 190., 190., 0.1, 0.05};
 132   fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, d#phi_{trk},d#eta_{trk}", ndims, bins, xmin, xmax);
 133   fOutputList->Add(fHnOutput);
 134
 135
 136
 137   PostData(1, fOutputList);
 138 }
 139
 140 //________________________________________________________________________
 141 void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
 142 {
 143   // Main loop, called for each event.
 144
 145   // event trigger selection
 146   Bool_t isSelected = 0;
 147   if(fPeriod.Contains("11a")){
 148     if(fTrigBit.Contains("kEMC"))
 149       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 150     else
 151       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 152   }
 153   else{
 154     if(fTrigBit.Contains("kEMC"))
 155       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
 156     else
 157       isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
 158   }
 159   if(!isSelected )
 160         return;
 161
 162   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 163   if (!fESD) {
 164     printf("ERROR: fESD not available\n");
 165     return;
 166   }
 167
 168   fEvtSel->Fill(0);
 169
 170   AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 171   if(!pv)
 172     return;
 173   fPVtxZ->Fill(pv->GetZ());
 174   if(TMath::Abs(pv->GetZ())>15)
 175     return;
 176
 177   fEvtSel->Fill(1);
 178
 179   if (!fTracks)
 180     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
 181   // Track loop to fill a pT spectrum
 182   const Int_t Ntracks = fTracks->GetEntriesFast();
 183   for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
 184     //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 185     //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
 186     AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks));
 187     if (!track)
 188       continue;
 189     if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
 190       fSelPrimTracks->Add(track);
 191       //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
 192     }
 193   }
 194
 195   if(!fIsTrain){
 196     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 197       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 198         break;
 199       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 200     }
 201   }
 202   fESD->GetEMCALClusters(fCaloClusters);
 203   fEMCalCells = fESD->GetEMCALCells();
 204
 205
 206   FillClusHists();
 207
 208   fCaloClusters->Clear();
 209   fSelPrimTracks->Clear();
 210   PostData(1, fOutputList);
 211 }
 212
 213 //________________________________________________________________________
 214 void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
 215 {
 216   // Fill cluster histograms.
 217
 218   if(!fCaloClusters)
 219     return;
 220   const Int_t nclus = fCaloClusters->GetEntries();
 221   if(nclus==0)
 222     return;
 223   Double_t maxE = 0;
 224   for(Int_t ic=0;ic<nclus;ic++){
 225     maxE=0;
 226     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 227     if(!c)
 228       continue;
 229     if(!c->IsEMCAL())
 230       continue;
 231     if(c->E()<fECut)
 232       continue;
 233     Short_t id;
 234     Double_t Emax = GetMaxCellEnergy( c, id);
 235     Double_t Ecross = GetCrossEnergy( c, id);
 236     if((1-Ecross/Emax)>fExoticCut)
 237       continue;
 238     Float_t clsPos[3] = {0,0,0};
 239     c->GetPosition(clsPos);
 240     TVector3 clsVec(clsPos);
 241     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 242     Float_t ceiso, cephiband, cecore;
 243     Float_t triso, trphiband, trcore;
 244     Float_t alliso, allphiband, allcore;
 245     Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
 246     Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
 247     GetCeIso(clsVec, ceiso, cephiband, cecore);
 248     GetTrIso(clsVec, triso, trphiband, trcore);
 249     alliso = ceiso + triso;
 250     allphiband = cephiband + trphiband;
 251     allcore = cecore + trcore;
 252     Float_t ceisoue =  cephiband/phibandArea*netConeArea;
 253     Float_t trisoue =  trphiband/phibandArea*netConeArea;
 254     Float_t allisoue =  allphiband/phibandArea*netConeArea;
 255     const Int_t ndims =   fNDimensions;
 256     Double_t outputValues[ndims];
 257     outputValues[0] = Et;
 258     outputValues[1] = c->GetM02();
 259     outputValues[2] = ceiso-cecore-ceisoue;
 260     outputValues[3] = triso-trisoue;
 261     outputValues[4] = alliso-cecore-allisoue;
 262     outputValues[5] = c->GetTrackDx();
 263     outputValues[6] = c->GetTrackDz();
 264     fHnOutput->Fill(outputValues);
 265     if(c->E()>maxE)
 266       maxE = c->E();
 267   }
 268   fNClusHighClusE->Fill(maxE,nclus);
 269 }
 270
 271 //________________________________________________________________________
 272 void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 273 {
 274   // Get cell isolation.
 275
 276   if(!fEMCalCells)
 277     return;
 278   const Int_t ncells = fEMCalCells->GetNumberOfCells();
 279   Float_t totiso=0;
 280   Float_t totband=0;
 281   Float_t totcore=0;
 282   Float_t etacl = vec.Eta();
 283   Float_t phicl = vec.Phi();
 284   Float_t thetacl = vec.Theta();
 285   if(phicl<0)
 286     phicl+=TMath::TwoPi();
 287   Int_t absid = -1;
 288   Float_t eta=-1, phi=-1;
 289   for(int icell=0;icell<ncells;icell++){
 290     absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
 291     if(!fGeom)
 292       return;
 293     fGeom->EtaPhiFromIndex(absid,eta,phi);
 294     Float_t dphi = TMath::Abs(phi-phicl);
 295     Float_t deta = TMath::Abs(eta-etacl);
 296     Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 297     Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 298     if(R<fIsoConeR){
 299       totiso += etcell;
 300       if(R<0.04)
 301         totcore += etcell;
 302     }
 303     else{
 304       if(dphi>fIsoConeR)
 305         continue;
 306       if(deta<fIsoConeR)
 307         continue;
 308       totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
 309     }
 310   }
 311   iso = totiso;
 312   phiband = totband;
 313   core = totcore;
 314 }
 315
 316 //________________________________________________________________________
 317 void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
 318 {
 319   // Get track isolation.
 320
 321   if(!fSelPrimTracks)
 322     return;
 323   const Int_t ntracks = fSelPrimTracks->GetEntries();
 324   Double_t totiso=0;
 325   Double_t totband=0;
 326   Double_t totcore=0;
 327   Double_t etacl = vec.Eta();
 328   Double_t phicl = vec.Phi();
 329   if(phicl<0)
 330     phicl+=TMath::TwoPi();
 331   for(int itrack=0;itrack<ntracks;itrack++){
 332     AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
 333     if(!track)
 334       continue;
 335     Double_t dphi = TMath::Abs(track->Phi()-phicl);
 336     Double_t deta = TMath::Abs(track->Eta()-etacl);
 337     Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
 338     Double_t pt = track->Pt();
 339     if(R<fIsoConeR){
 340       totiso += track->Pt();
 341       if(R<0.04)
 342         totcore += pt;
 343     }
 344     else{
 345       if(dphi>fIsoConeR)
 346         continue;
 347       if(deta<fIsoConeR)
 348         continue;
 349       totband += track->Pt();
 350     }
 351   }
 352   iso = totiso;
 353   phiband = totband;
 354   core = totcore;
 355 }
 356
 357 //________________________________________________________________________
 358 Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 359 {
 360   // Calculate the energy of cross cells around the leading cell.
 361
 362   AliVCaloCells *cells = 0;
 363   cells = fESD->GetEMCALCells();
 364   if (!cells)
 365     return 0;
 366
 367   if (!fGeom)
 368     return 0;
 369
 370   Int_t iSupMod = -1;
 371   Int_t iTower  = -1;
 372   Int_t iIphi   = -1;
 373   Int_t iIeta   = -1;
 374   Int_t iphi    = -1;
 375   Int_t ieta    = -1;
 376   Int_t iphis   = -1;
 377   Int_t ietas   = -1;
 378
 379   Double_t crossEnergy = 0;
 380
 381   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 382   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 383
 384   Int_t ncells = cluster->GetNCells();
 385   for (Int_t i=0; i<ncells; i++) {
 386     Int_t cellAbsId = cluster->GetCellAbsId(i);
 387     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 388     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 389     Int_t aphidiff = TMath::Abs(iphi-iphis);
 390     if (aphidiff>1)
 391       continue;
 392     Int_t aetadiff = TMath::Abs(ieta-ietas);
 393     if (aetadiff>1)
 394       continue;
 395     if ( (aphidiff==1 && aetadiff==0) ||
 396         (aphidiff==0 && aetadiff==1) ) {
 397       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 398     }
 399   }
 400
 401   return crossEnergy;
 402 }
 403
 404 //________________________________________________________________________
 405 Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 406 {
 407   // Get maximum energy of attached cell.
 408
 409   id = -1;
 410
 411   AliVCaloCells *cells = 0;
 412   cells = fESD->GetEMCALCells();
 413   if (!cells)
 414     return 0;
 415
 416   Double_t maxe = 0;
 417   Int_t ncells = cluster->GetNCells();
 418   for (Int_t i=0; i<ncells; i++) {
 419     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 420     if (e>maxe) {
 421       maxe = e;
 422       id   = cluster->GetCellAbsId(i);
 423     }
 424   }
 425   return maxe;
 426 }
 427
 428 //________________________________________________________________________
 429 void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
 430 {
 431   // Called once at the end of the query.
 432 }