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