PWGGA/EMCALTasks/AliAnalysisTaskTrgContam.cxx

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