PWGGA/EMCALTasks/AliAnalysisTaskTrgContam.cxx

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