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