PWGGA/EMCALTasks/AliAnalysisTaskTrgContam.cxx

   1 // $Id$
   2 //
   3 // Task to compute the trigger contamination by exotics.
   4 //
   5 // Authors: M.Cosentino
   6
   7 #include "TChain.h"
   8 #include "TTree.h"
   9 #include "TH1F.h"
  10 #include "TH2F.h"
  11 #include "TCanvas.h"
  12
  13 #include "AliAnalysisTask.h"
  14 #include "AliAnalysisManager.h"
  15 #include "AliESDEvent.h"
  16 #include "AliESDHeader.h"
  17 #include "AliESDUtils.h"
  18 #include "AliESDInputHandler.h"
  19 #include "AliESDpid.h"
  20 #include "AliKFParticle.h"
  21 #include "AliMCEventHandler.h"
  22 #include "AliMCEvent.h"
  23 #include "AliStack.h"
  24 #include "AliESDtrackCuts.h"
  25 #include "AliESDv0.h"
  26 #include "AliV0vertexer.h"
  27 #include "AliESDCaloCluster.h"
  28 #include "AliESDCaloCells.h"
  29 #include "AliEMCALGeometry.h"
  30 #include "AliEMCALRecoUtils.h"
  31 #include "TLorentzVector.h"
  32 #include "AliVCluster.h"
  33 #include "AliAnalysisTaskTrgContam.h"
  34 #include "TFile.h"
  35
  36 ClassImp(AliAnalysisTaskTrgContam)
  37
  38 AliAnalysisTaskTrgContam::AliAnalysisTaskTrgContam() :
  39   AliAnalysisTaskSE(),
  40   fCaloClusters(0),
  41   fEMCalCells(0),
  42   fGeom(0x0),
  43   fGeoName("EMCAL_COMPLETEV1"),
  44   fPeriod("LHC11c"),
  45   fIsTrain(0),
  46   fTrigThresh(4.8),
  47   fExoticCut(0.97),
  48   fESD(0),
  49   fOutputList(0),
  50   fEvtSel(0),
  51   fClusEt(0),
  52   fClusEtTM(0),
  53   fClusEtLead(0),
  54   fClusEtSubLead(0),
  55   fClusEtLeadTM(0),
  56   fClusEtSubLeadTM(0),
  57   fClusEtExotic(0),
  58   fClusEtExoticTM(0),
  59   fClusEtSingleExotic(0),
  60   fCellEnergy(0),
  61   fM02Et(0),
  62   fM02EtTM(0),
  63   fM02EtExot(0),
  64   fM02EtExotTM(0)
  65 {
  66   // Default constructor.
  67 }
  68
  69 //________________________________________________________________________
  70 AliAnalysisTaskTrgContam::AliAnalysisTaskTrgContam(const char *name) :
  71   AliAnalysisTaskSE(name),
  72   fCaloClusters(0),
  73   fEMCalCells(0),
  74   fGeom(0x0),
  75   fGeoName("EMCAL_COMPLETEV1"),
  76   fPeriod("LHC11c"),
  77   fIsTrain(0),
  78   fTrigThresh(4.8),
  79   fExoticCut(0.97),
  80   fESD(0),
  81   fOutputList(0),
  82   fEvtSel(0),
  83   fClusEt(0),
  84   fClusEtTM(0),
  85   fClusEtLead(0),
  86   fClusEtSubLead(0),
  87   fClusEtLeadTM(0),
  88   fClusEtSubLeadTM(0),
  89   fClusEtExotic(0),
  90   fClusEtExoticTM(0),
  91   fClusEtSingleExotic(0),
  92   fCellEnergy(0),
  93   fM02Et(0),
  94   fM02EtTM(0),
  95   fM02EtExot(0),
  96   fM02EtExotTM(0)
  97 {
  98   // Constructor
  99
 100   // Define input and output slots here
 101   // Input slot #0 works with a TChain
 102   DefineInput(0, TChain::Class());
 103   // Output slot #0 id reserved by the base class for AOD
 104   // Output slot #1 writes into a TH1 container
 105   DefineOutput(1, TList::Class());
 106 }
 107
 108 //________________________________________________________________________
 109 void AliAnalysisTaskTrgContam::UserCreateOutputObjects()
 110 {
 111   // Create histograms, called once.
 112
 113   fCaloClusters = new TRefArray();
 114
 115   fOutputList = new TList();
 116   fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
 117
 118   fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data());
 119
 120   fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
 121   fOutputList->Add(fEvtSel);
 122
 123   fClusEt = new TH1F("hClusEt","Clusters E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 124   fOutputList->Add(fClusEt);
 125
 126   fClusEtTM = new TH1F("hClusEtTM","Clusters (track-matched) E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 127   fOutputList->Add(fClusEtTM);
 128
 129   fClusEtLead = new TH1F("hClusEtLead","Clusters (leading-trig) E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 130   fOutputList->Add(fClusEtLead);
 131
 132   fClusEtSubLead = new TH1F("hClusEtSubLead","Clusters (subleading-trig) E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 133   fOutputList->Add(fClusEtSubLead);
 134
 135   fClusEtLeadTM = new TH1F("hClusEtLeadTM","Clusters (leading-trig, TM) E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 136   fOutputList->Add(fClusEtLeadTM);
 137
 138   fClusEtSubLeadTM = new TH1F("hClusEtSubLeadTM","Clusters (subleading-trig, TM) E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 139   fOutputList->Add(fClusEtSubLeadTM);
 140
 141   fClusEtExotic = new TH1F("hClusEtExotic","Exotic trigger clusters  E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 142   fOutputList->Add(fClusEtExotic);
 143
 144   fClusEtExoticTM = new TH1F("hClusEtExoticTM","Exotic trigger clusters (TM)  E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 145   fOutputList->Add(fClusEtExoticTM);
 146
 147   fClusEtSingleExotic = new TH1F("hClusEtSingleExotic","Exotic trigger clusters (only this above thrs.)  E_{T} ;E_{T} ;dN/dE_{T}",400,0,200);
 148   fOutputList->Add(fClusEtSingleExotic);
 149
 150   fCellEnergy = new TH1F("hCellE","cell energy spectrum;E_{cell} (GeV);entries",200,0,20);
 151   fOutputList->Add(fCellEnergy);
 152
 153   fM02Et = new TH2F("hM02Et","#lambda_{0}^{2} vs. E_{T} for trigger clusters ;E_{T} ;#lambda_{0}^{2}",400,0,200, 400,0,4);
 154   fOutputList->Add(fM02Et);
 155
 156   fM02EtTM = new TH2F("hM02EtTM","#lambda_{0}^{2} vs. E_{T} for trigger clusters(TM) ;E_{T} ;#lambda_{0}^{2}",400,0,200, 400,0,4);
 157   fOutputList->Add(fM02EtTM);
 158
 159   fM02EtExot = new TH2F("hM02EtExot","#lambda_{0}^{2} vs. E_{T} for trigger clusters(Exotic) ;E_{T} ;#lambda_{0}^{2}",400,0,200, 400,0,4);
 160   fOutputList->Add(fM02EtExot);
 161
 162   fM02EtExotTM = new TH2F("hM02EtExotTM","#lambda_{0}^{2} vs. E_{T} for trigger clusters(TM+Exotic) ;E_{T} ;#lambda_{0}^{2}",400,0,200, 400,0,4);
 163   fOutputList->Add(fM02EtExotTM);
 164
 165   PostData(1, fOutputList);
 166 }
 167
 168 //________________________________________________________________________
 169 void AliAnalysisTaskTrgContam::UserExec(Option_t *)
 170 {
 171   // User exec. Called once per event.
 172
 173   Bool_t isSelected = 0;
 174   if(fPeriod.Contains("11a"))
 175     isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
 176   else
 177     isSelected =  ((((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ||
 178                    (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral));
 179   if(!isSelected )
 180     return;
 181
 182   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 183   if (!fESD) {
 184     printf("ERROR: fESD not available\n");
 185     return;
 186   }
 187
 188   fEvtSel->Fill(0);
 189   if (0) {
 190     AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
 191     if(!pv)
 192       return;
 193     if(TMath::Abs(pv->GetZ())>15)
 194       return;
 195   }
 196   fEvtSel->Fill(1);
 197   if (!fIsTrain) {
 198     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 199       if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
 200         break;
 201       fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
 202     }
 203   }
 204   fESD->GetEMCALClusters(fCaloClusters);
 205   fEMCalCells = fESD->GetEMCALCells();
 206   for(int i=0;i<fEMCalCells->GetNumberOfCells();i++){
 207     Double_t e = fEMCalCells->GetCellAmplitude(TMath::Abs(fEMCalCells->GetAmplitude(i)));
 208     fCellEnergy->Fill(e);
 209   }
 210   FillClusHists();
 211
 212   fCaloClusters->Clear();
 213   PostData(1, fOutputList);
 214 }
 215
 216 //________________________________________________________________________
 217 void AliAnalysisTaskTrgContam::FillClusHists()
 218 {
 219   // Fill cluster histograms.
 220
 221   if(!fCaloClusters)
 222     return;
 223   const Int_t nclus = fCaloClusters->GetEntries();
 224   if(nclus==0)
 225     return;
 226   Double_t EtArray[nclus];
 227   Bool_t isTM[nclus];
 228   Bool_t isEx[nclus];
 229   Int_t index[nclus];
 230   Int_t nthresholds = 0;
 231   for(Int_t ic=0;ic<nclus;ic++){
 232     EtArray[ic]=0;
 233     isTM[ic] = 0;
 234     isEx[ic] = 0;
 235     index[ic]=0;
 236     AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
 237     if(!c)
 238       continue;
 239     if(!c->IsEMCAL())
 240       continue;
 241     if(c->E()<fTrigThresh)
 242       continue;
 243     nthresholds++;
 244     Short_t id;
 245     Double_t Emax = GetMaxCellEnergy( c, id);
 246     Double_t Ecross = GetCrossEnergy( c, id);
 247     if((1-Ecross/Emax)>fExoticCut)
 248       isEx[ic] = 1;
 249     Float_t clsPos[3] = {0,0,0};
 250     c->GetPosition(clsPos);
 251     TVector3 clsVec(clsPos);
 252     Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
 253     EtArray[ic] = Et;
 254     fClusEt->Fill(Et);
 255     fM02Et->Fill(Et, c->GetM02());
 256     if(isEx[ic]){
 257       fClusEtExotic->Fill(Et);
 258       fM02EtExot->Fill(Et,c->GetM02());
 259      }
 260     Double_t dR = TMath::Sqrt(pow(c->GetTrackDx(),2)+pow(c->GetTrackDz(),2));
 261     if(dR<0.025){
 262       isTM[ic]=1;
 263       fClusEtTM->Fill(Et);
 264       fM02EtTM->Fill(Et, c->GetM02());
 265       if(isEx[ic]){
 266         fClusEtExoticTM->Fill(Et);
 267         fM02EtExotTM->Fill(Et,c->GetM02());
 268       }
 269     }
 270   }
 271   TMath::Sort(nclus,EtArray,index, kTRUE);
 272   if(EtArray[index[0]]>0){
 273     fClusEtLead->Fill(EtArray[index[0]]);
 274     if(nthresholds==1 && isEx[index[0]])
 275        fClusEtSingleExotic->Fill(EtArray[index[0]]);
 276   }
 277   if(nclus>1)if(EtArray[index[1]]>0)
 278     fClusEtSubLead->Fill(EtArray[index[1]]);
 279   if(isTM[index[0]] && EtArray[index[0]]>0)
 280     fClusEtLeadTM->Fill(EtArray[index[0]]);
 281   if(nclus>1)if(isTM[index[1]] && EtArray[index[1]]>0)
 282     fClusEtSubLeadTM->Fill(EtArray[index[1]]);
 283 }
 284
 285 //________________________________________________________________________
 286 Double_t AliAnalysisTaskTrgContam::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
 287 {
 288   // Calculate the energy of cross cells around the leading cell.
 289
 290   AliVCaloCells *cells = 0;
 291   cells = fESD->GetEMCALCells();
 292   if (!cells)
 293     return 0;
 294
 295   if (!fGeom)
 296     return 0;
 297
 298   Int_t iSupMod = -1;
 299   Int_t iTower  = -1;
 300   Int_t iIphi   = -1;
 301   Int_t iIeta   = -1;
 302   Int_t iphi    = -1;
 303   Int_t ieta    = -1;
 304   Int_t iphis   = -1;
 305   Int_t ietas   = -1;
 306
 307   Double_t crossEnergy = 0;
 308
 309   fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
 310   fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
 311
 312   Int_t ncells = cluster->GetNCells();
 313   for (Int_t i=0; i<ncells; i++) {
 314     Int_t cellAbsId = cluster->GetCellAbsId(i);
 315     fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
 316     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
 317     Int_t aphidiff = TMath::Abs(iphi-iphis);
 318     if (aphidiff>1)
 319       continue;
 320     Int_t aetadiff = TMath::Abs(ieta-ietas);
 321     if (aetadiff>1)
 322       continue;
 323     if ( (aphidiff==1 && aetadiff==0) ||
 324         (aphidiff==0 && aetadiff==1) ) {
 325       crossEnergy += cells->GetCellAmplitude(cellAbsId);
 326     }
 327   }
 328
 329   return crossEnergy;
 330 }
 331
 332 //________________________________________________________________________
 333 Double_t AliAnalysisTaskTrgContam ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
 334 {
 335   // Get maximum energy of attached cell.
 336
 337   id = -1;
 338
 339   AliVCaloCells *cells = 0;
 340   cells = fESD->GetEMCALCells();
 341   if (!cells)
 342     return 0;
 343
 344   Double_t maxe = 0;
 345   Int_t ncells = cluster->GetNCells();
 346   for (Int_t i=0; i<ncells; i++) {
 347     Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
 348     if (e>maxe) {
 349       maxe = e;
 350       id   = cluster->GetCellAbsId(i);
 351     }
 352   }
 353   return maxe;
 354 }
 355
 356 //________________________________________________________________________
 357 void AliAnalysisTaskTrgContam::Terminate(Option_t *)
 358 {
 359   // Called once at the end of the query
 360 }