1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Authors: Friederike Bock, Baldo Sahlmueller *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class handling all kinds of selection cuts for
19 // Photon from EMCAL clusters
20 //---------------------------------------------
21 ////////////////////////////////////////////////
23 #include "AliCaloPhotonCuts.h"
24 #include "AliAnalysisManager.h"
25 #include "AliInputEventHandler.h"
26 #include "AliMCEventHandler.h"
27 #include "AliAODHandler.h"
32 #include "AliAODConversionMother.h"
33 #include "TObjString.h"
34 #include "AliAODEvent.h"
35 #include "AliESDEvent.h"
36 #include "AliCentrality.h"
40 #include "AliV0ReaderV1.h"
41 #include "AliAODMCParticle.h"
42 #include "AliAODMCHeader.h"
48 ClassImp(AliCaloPhotonCuts)
51 const char* AliCaloPhotonCuts::fgkCutNames[AliCaloPhotonCuts::kNCuts] = {
57 "DistanceToBadChannel", //5
67 "MaximumDispersion", //15
72 //________________________________________________________________________
73 AliCaloPhotonCuts::AliCaloPhotonCuts(const char *name,const char *title) :
74 AliAnalysisCuts(name,title),
81 fMinDistanceToBadChannel(0),
83 fMinDistTrackToCluster(0),
96 fHistAcceptanceCuts(NULL),
97 fHistClusterIdentificationCuts(NULL),
98 fHistClusterEtavsPhiBeforeAcc(NULL),
99 fHistClusterEtavsPhiAfterAcc(NULL),
100 fHistClusterEtavsPhiAfterQA(NULL),
101 fHistDistanceToBadChannelBeforeAcc(NULL),
102 fHistDistanceToBadChannelAfterAcc(NULL),
103 fHistClusterTimevsEBeforeQA(NULL),
104 fHistClusterTimevsEAfterQA(NULL),
105 fHistExoticCellBeforeQA(NULL),
106 fHistExoticCellAfterQA(NULL),
107 fHistDistanceTrackToClusterBeforeQA(NULL),
108 fHistDistanceTrackToClusterAfterQA(NULL),
109 fHistEnergyOfClusterBeforeQA(NULL),
110 fHistEnergyOfClusterAfterQA(NULL),
111 fHistNCellsBeforeQA(NULL),
112 fHistNCellsAfterQA(NULL),
113 fHistM02BeforeQA(NULL),
114 fHistM02AfterQA(NULL),
115 fHistM20BeforeQA(NULL),
116 fHistM20AfterQA(NULL),
117 fHistDispersionBeforeQA(NULL),
118 fHistDispersionAfterQA(NULL),
119 fHistNLMBeforeQA(NULL),
120 fHistNLMAfterQA(NULL)
122 for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
123 fCutString=new TObjString((GetCutNumber()).Data());
126 //________________________________________________________________________
127 AliCaloPhotonCuts::AliCaloPhotonCuts(const AliCaloPhotonCuts &ref) :
128 AliAnalysisCuts(ref),
130 fClusterType(ref.fClusterType),
131 fMinEtaCut(ref.fMinEtaCut),
132 fMaxEtaCut(ref.fMaxEtaCut),
133 fMinPhiCut(ref.fMinPhiCut),
134 fMaxPhiCut(ref.fMaxPhiCut),
135 fMinDistanceToBadChannel(ref.fMinDistanceToBadChannel),
136 fMaxTimeDiff(ref.fMaxTimeDiff),
137 fMinDistTrackToCluster(ref.fMinDistTrackToCluster),
138 fExoticCell(ref.fExoticCell),
139 fMinEnergy(ref.fMinEnergy),
140 fMinNCells(ref.fMinNCells),
141 fMaxM02(ref.fMaxM02),
142 fMinM02(ref.fMinM02),
143 fMaxM20(ref.fMaxM20),
144 fMinM20(ref.fMinM20),
145 fMaxDispersion(ref.fMaxDispersion),
146 fMinNLM(ref.fMinNLM),
147 fMaxNLM(ref.fMaxNLM),
150 fHistAcceptanceCuts(NULL),
151 fHistClusterIdentificationCuts(NULL),
152 fHistClusterEtavsPhiBeforeAcc(NULL),
153 fHistClusterEtavsPhiAfterAcc(NULL),
154 fHistClusterEtavsPhiAfterQA(NULL),
155 fHistDistanceToBadChannelBeforeAcc(NULL),
156 fHistDistanceToBadChannelAfterAcc(NULL),
157 fHistClusterTimevsEBeforeQA(NULL),
158 fHistClusterTimevsEAfterQA(NULL),
159 fHistExoticCellBeforeQA(NULL),
160 fHistExoticCellAfterQA(NULL),
161 fHistDistanceTrackToClusterBeforeQA(NULL),
162 fHistDistanceTrackToClusterAfterQA(NULL),
163 fHistEnergyOfClusterBeforeQA(NULL),
164 fHistEnergyOfClusterAfterQA(NULL),
165 fHistNCellsBeforeQA(NULL),
166 fHistNCellsAfterQA(NULL),
167 fHistM02BeforeQA(NULL),
168 fHistM02AfterQA(NULL),
169 fHistM20BeforeQA(NULL),
170 fHistM20AfterQA(NULL),
171 fHistDispersionBeforeQA(NULL),
172 fHistDispersionAfterQA(NULL),
173 fHistNLMBeforeQA(NULL),
174 fHistNLMAfterQA(NULL)
177 for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
178 fCutString=new TObjString((GetCutNumber()).Data());
183 //________________________________________________________________________
184 AliCaloPhotonCuts::~AliCaloPhotonCuts() {
186 //Deleting fHistograms leads to seg fault it it's added to output collection of a task
188 // delete fHistograms;
189 // fHistograms = NULL;
190 if(fCutString != NULL){
196 //________________________________________________________________________
197 void AliCaloPhotonCuts::InitCutHistograms(TString name){
199 // Initialize Cut Histograms for QA (only initialized and filled if function is called)
200 TH1::AddDirectory(kFALSE);
202 if(fHistograms != NULL){
206 if(fHistograms==NULL){
207 fHistograms=new TList();
208 fHistograms->SetOwner(kTRUE);
209 if(name=="")fHistograms->SetName(Form("CaloCuts_%s",GetCutNumber().Data()));
210 else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
214 fHistCutIndex=new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",5,-0.5,4.5);
215 fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in");
216 fHistCutIndex->GetXaxis()->SetBinLabel(kDetector+1,"detector");
217 fHistCutIndex->GetXaxis()->SetBinLabel(kAcceptance+1,"acceptance");
218 fHistCutIndex->GetXaxis()->SetBinLabel(kClusterQuality+1,"cluster QA");
219 fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonOut+1,"out");
220 fHistograms->Add(fHistCutIndex);
223 fHistAcceptanceCuts=new TH1F(Form("AcceptanceCuts %s",GetCutNumber().Data()),"AcceptanceCuts",5,-0.5,4.5);
224 fHistAcceptanceCuts->GetXaxis()->SetBinLabel(1,"in");
225 fHistAcceptanceCuts->GetXaxis()->SetBinLabel(2,"eta");
226 fHistAcceptanceCuts->GetXaxis()->SetBinLabel(3,"phi");
227 fHistAcceptanceCuts->GetXaxis()->SetBinLabel(4,"distance to bad channel");
228 fHistAcceptanceCuts->GetXaxis()->SetBinLabel(5,"out");
229 fHistograms->Add(fHistAcceptanceCuts);
232 fHistClusterIdentificationCuts =new TH1F(Form("ClusterQualityCuts %s",GetCutNumber().Data()),"ClusterQualityCuts",10,-0.5,9.5);
233 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(1,"in");
234 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(2,"timing");
235 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(3,"track matching");
236 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(4,"Exotics");
237 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(5,"minimum energy");
238 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(6,"M02");
239 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(7,"M20");
240 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(8,"dispersion");
241 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(9,"NLM");
242 fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(10,"out");
243 fHistograms->Add(fHistClusterIdentificationCuts);
245 // Acceptance related histogramms
246 fHistClusterEtavsPhiBeforeAcc=new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
247 fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
248 fHistClusterEtavsPhiAfterAcc=new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
249 fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
250 fHistClusterEtavsPhiAfterQA=new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
251 fHistograms->Add(fHistClusterEtavsPhiAfterQA);
252 fHistDistanceToBadChannelBeforeAcc = new TH1F(Form("DistanceToBadChannel_beforeAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_beforeAcceptance",200,0,40);
253 fHistograms->Add(fHistDistanceToBadChannelBeforeAcc);
254 fHistDistanceToBadChannelAfterAcc = new TH1F(Form("DistanceToBadChannel_afterAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_afterAcceptance",200,0,40);
255 fHistograms->Add(fHistDistanceToBadChannelAfterAcc);
257 // Cluster quality related histograms
258 fHistClusterTimevsEBeforeQA=new TH2F(Form("ClusterTimeVsE_beforeClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_beforeClusterQA",400,-10e-6,10e-6,100,0.,40);
259 fHistograms->Add(fHistClusterTimevsEBeforeQA);
260 fHistClusterTimevsEAfterQA=new TH2F(Form("ClusterTimeVsE_afterClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_afterClusterQA",400,-10e-6,10e-6,100,0.,40);
261 fHistograms->Add(fHistClusterTimevsEAfterQA);
262 fHistExoticCellBeforeQA=new TH2F(Form("ExoticCell_beforeClusterQA %s",GetCutNumber().Data()),"ExoticCell_beforeClusterQA",400,0,40,50,0.75,1);
263 fHistograms->Add(fHistExoticCellBeforeQA);
264 fHistExoticCellAfterQA=new TH2F(Form("ExoticCell_afterClusterQA %s",GetCutNumber().Data()),"ExoticCell_afterClusterQA",400,0,40,50,0.75,1);
265 fHistograms->Add(fHistExoticCellAfterQA);
266 fHistDistanceTrackToClusterBeforeQA = new TH1F(Form("DistanceToTrack_beforeClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_beforeClusterQA",200,0,40);
267 fHistograms->Add(fHistDistanceTrackToClusterBeforeQA);
268 fHistDistanceTrackToClusterAfterQA = new TH1F(Form("DistanceToTrack_afterClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_afterClusterQA",200,0,40);
269 fHistograms->Add(fHistDistanceTrackToClusterAfterQA);
270 fHistEnergyOfClusterBeforeQA = new TH1F(Form("EnergyOfCluster_beforeClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_beforeClusterQA",300,0,30);
271 fHistograms->Add(fHistEnergyOfClusterBeforeQA);
272 fHistEnergyOfClusterAfterQA = new TH1F(Form("EnergyOfCluster_afterClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_afterClusterQA",300,0,30);
273 fHistograms->Add(fHistEnergyOfClusterAfterQA);
274 fHistNCellsBeforeQA = new TH1F(Form("NCellPerCluster_beforeClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_beforeClusterQA",50,0,50);
275 fHistograms->Add(fHistNCellsBeforeQA);
276 fHistNCellsAfterQA = new TH1F(Form("NCellPerCluster_afterClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_afterClusterQA",50,0,50);
277 fHistograms->Add(fHistNCellsAfterQA);
278 fHistM02BeforeQA = new TH1F(Form("M02_beforeClusterQA %s",GetCutNumber().Data()),"M02_beforeClusterQA",100,0,5);
279 fHistograms->Add(fHistM02BeforeQA);
280 fHistM02AfterQA = new TH1F(Form("M02_afterClusterQA %s",GetCutNumber().Data()),"M02_afterClusterQA",100,0,5);
281 fHistograms->Add(fHistM02AfterQA);
282 fHistM20BeforeQA = new TH1F(Form("M20_beforeClusterQA %s",GetCutNumber().Data()),"M20_beforeClusterQA",100,0,2.5);
283 fHistograms->Add(fHistM20BeforeQA);
284 fHistM20AfterQA = new TH1F(Form("M20_afterClusterQA %s",GetCutNumber().Data()),"M20_afterClusterQA",100,0,2.5);
285 fHistograms->Add(fHistM20AfterQA);
286 fHistDispersionBeforeQA = new TH1F(Form("Dispersion_beforeClusterQA %s",GetCutNumber().Data()),"Dispersion_beforeClusterQA",100,0,4);
287 fHistograms->Add(fHistDispersionBeforeQA);
288 fHistDispersionAfterQA = new TH1F(Form("Dispersion_afterClusterQA %s",GetCutNumber().Data()),"Dispersion_afterClusterQA",100,0,4);
289 fHistograms->Add(fHistDispersionAfterQA);
290 fHistNLMBeforeQA = new TH1F(Form("NLM_beforeClusterQA %s",GetCutNumber().Data()),"NLM_beforeClusterQA",10,0,10);
291 fHistograms->Add(fHistNLMBeforeQA);
292 fHistNLMAfterQA = new TH1F(Form("NLM_afterClusterQA %s",GetCutNumber().Data()),"NLM_afterClusterQA",10,0,10);
293 fHistograms->Add(fHistNLMAfterQA);
295 TH1::AddDirectory(kTRUE);
299 ///________________________________________________________________________
300 Bool_t AliCaloPhotonCuts::ClusterIsSelectedMC(TParticle *particle,AliStack *fMCStack,Bool_t checkForConvertedGamma){
301 // MonteCarlo Photon Selection
303 if(!fMCStack)return kFALSE;
305 if (particle->GetPdgCode() == 22){
308 if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) )
311 if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) )
315 if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
316 return kFALSE; // no photon as mothers!
319 if(particle->GetMother(0) >= fMCStack->GetNprimary()){
320 return kFALSE; // the gamma has a mother, and it is not a primary particle
323 if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
325 // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
326 TParticle* ePos = NULL;
327 TParticle* eNeg = NULL;
329 if(particle->GetNDaughters() >= 2){
330 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
331 TParticle *tmpDaughter = fMCStack->Particle(daughterIndex);
332 if(tmpDaughter->GetUniqueID() == 5){
333 if(tmpDaughter->GetPdgCode() == 11){
335 } else if(tmpDaughter->GetPdgCode() == -11){
342 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
346 if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
347 return kFALSE; // no reconstruction below the Pt cut
350 if( ePos->Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) ||
351 eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) )
354 if(fEtaCutMin > -0.1){
355 if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) ||
356 (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) )
361 return kFALSE; // cuts on distance from collision point
364 if(abs(ePos->Vz()) > fMaxZ){
365 return kFALSE; // outside material
367 if(abs(eNeg->Vz()) > fMaxZ){
368 return kFALSE; // outside material
371 if( ePos->R() <= ((abs(ePos->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
372 return kFALSE; // line cut to exclude regions where we do not reconstruct
373 } else if ( fEtaCutMin != -0.1 && ePos->R() >= ((abs(ePos->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
377 if( eNeg->R() <= ((abs(eNeg->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
378 return kFALSE; // line cut to exclude regions where we do not reconstruct
379 } else if ( fEtaCutMin != -0.1 && eNeg->R() >= ((abs(eNeg->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
384 //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
388 ///________________________________________________________________________
389 Bool_t AliCaloPhotonCuts::ClusterIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray,Bool_t checkForConvertedGamma){
390 // MonteCarlo Photon Selection
392 if(!aodmcArray)return kFALSE;
394 if (particle->GetPdgCode() == 22){
395 if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) )
398 if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) )
402 if(particle->GetMother() > -1){
403 if((static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){
404 return kFALSE; // no photon as mothers!
406 if(!(static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother()))->IsPrimary())){
407 return kFALSE; // the gamma has a mother, and it is not a primary particle
411 if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
413 // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
414 AliAODMCParticle* ePos = NULL;
415 AliAODMCParticle* eNeg = NULL;
417 if(particle->GetNDaughters() >= 2){
418 for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
419 AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(aodmcArray->At(daughterIndex));
420 if(!tmpDaughter) continue;
421 if(((tmpDaughter->GetMCProcessCode())) == 5){ // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
422 if(tmpDaughter->GetPdgCode() == 11){
424 } else if(tmpDaughter->GetPdgCode() == -11){
431 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
435 if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
436 return kFALSE; // no reconstruction below the Pt cut
439 if( ePos->Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) ||
440 eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) )
443 if(fEtaCutMin > -0.1){
444 if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) ||
445 (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) )
449 Double_t rPos = sqrt( (ePos->Xv()*ePos->Xv()) + (ePos->Yv()*ePos->Yv()) );
450 Double_t rNeg = sqrt( (eNeg->Xv()*eNeg->Xv()) + (eNeg->Yv()*eNeg->Yv()) );
453 return kFALSE; // cuts on distance from collision point
455 if(abs(ePos->Zv()) > fMaxZ){
456 return kFALSE; // outside material
458 if(abs(eNeg->Zv()) > fMaxZ){
459 return kFALSE; // outside material
462 if( rPos <= ((abs(ePos->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
463 return kFALSE; // line cut to exclude regions where we do not reconstruct
464 } else if ( fEtaCutMin != -0.1 && rPos >= ((abs(ePos->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
468 if( rNeg <= ((abs(eNeg->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
469 return kFALSE; // line cut to exclude regions where we do not reconstruct
470 } else if ( fEtaCutMin != -0.1 && rNeg >= ((abs(eNeg->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
475 //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
482 ///________________________________________________________________________
483 // This function selects the clusters based on their quality criteria
484 ///________________________________________________________________________
485 Bool_t AliCaloPhotonCuts::ClusterQualityCuts(AliVCluster* cluster, AliVEvent *event, Bool_t isMC)
486 { // Specific Photon Cuts
489 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex);
493 // Fill Histos before Cuts
494 if(fHistClusterTimevsEBeforeQA) fHistClusterTimevsEBeforeQA->Fill(cluster->GetTOF(), cluster->E());
495 // if(fHistExoticCellBeforeQA) fHistExoticCellBeforeQA->Fill(cluster->E(), );
496 if(fHistDistanceTrackToClusterBeforeQA) fHistDistanceTrackToClusterBeforeQA->Fill(cluster->GetEmcCpvDistance());
497 if(fHistEnergyOfClusterBeforeQA) fHistEnergyOfClusterBeforeQA->Fill(cluster->E());
498 if(fHistNCellsBeforeQA) fHistNCellsBeforeQA->Fill(cluster->GetNCells());
499 if(fHistM02BeforeQA) fHistM02BeforeQA->Fill(cluster->GetM02());
500 if(fHistM20BeforeQA) fHistM20BeforeQA->Fill(cluster->GetM20());
501 if(fHistDispersionBeforeQA) fHistDispersionBeforeQA->Fill(cluster->GetDispersion());
502 // if(fHistNLMBeforeQA) fHistNLMBeforeQA->Fill(cluster->GetNExMax());
504 // Check wether timing is ok
505 if(abs(cluster->GetTOF()) > fMaxTimeDiff && !isMC){
506 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //1
509 cutIndex++; //2, next cut
511 // Minimum distance to track
512 if(cluster->GetEmcCpvDistance() < fMinDistTrackToCluster){
513 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2
516 cutIndex++;//3, next cut
518 // exotic cell cut --IMPLEMENT LATER---
519 // if(!AcceptanceCuts(photon)){
520 // if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //3
523 cutIndex++; //4, next cut
525 // minimum cell energy cut
526 if(cluster->E() < fMinEnergy){
527 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //4
530 cutIndex++; //5, next cut
532 // minimum number of cells
533 if(cluster->GetNCells() < fMinNCells) {
534 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //5
537 cutIndex++; //6, next cut
540 if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 ) {
541 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //6
544 cutIndex++; //7, next cut
547 if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 ) {
548 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //7
551 cutIndex++; //8, next cut
554 if( cluster->GetDispersion()> fMaxDispersion) {
555 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //8
558 cutIndex++; //9, next cut
560 // NLM cut --IMPLEMENT LATER---
561 if( cluster->GetDispersion()> fMaxDispersion) {
562 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //9
565 cutIndex++; //9, next cut
567 // DONE with selecting photons
568 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //10
571 Double_t vertex[3] = {0};
572 event->GetPrimaryVertex()->GetXYZ(vertex);
573 // TLorentzvector with cluster
574 TLorentzVector clusterVector;
575 cluster->GetMomentum(clusterVector,vertex);
576 Double_t etaCluster = clusterVector.Eta();
577 Double_t phiCluster = clusterVector.Phi();
579 if(fHistClusterEtavsPhiAfterQA) fHistClusterEtavsPhiAfterQA->Fill(phiCluster,etaCluster);
580 if(fHistClusterTimevsEAfterQA) fHistClusterTimevsEAfterQA->Fill(cluster->GetTOF(), cluster->E());
581 // if(fHistExoticCellAfterQA) fHistExoticCellAfterQA->Fill(cluster->E(), );
582 if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(cluster->GetEmcCpvDistance());
583 if(fHistEnergyOfClusterAfterQA) fHistEnergyOfClusterAfterQA->Fill(cluster->E());
584 if(fHistNCellsAfterQA) fHistNCellsAfterQA->Fill(cluster->GetNCells());
585 if(fHistM02AfterQA) fHistM02AfterQA->Fill(cluster->GetM02());
586 if(fHistM20AfterQA) fHistM20AfterQA->Fill(cluster->GetM20());
587 if(fHistDispersionAfterQA) fHistDispersionAfterQA->Fill(cluster->GetDispersion());
588 // if(fHistNLMBeforeQA) fHistNLMAfterQA->Fill(cluster->GetNExMax());
595 ///________________________________________________________________________
596 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, Bool_t isMC)
598 //Selection of Reconstructed photon clusters with Calorimeters
600 FillClusterCutIndex(kPhotonIn);
602 Double_t vertex[3] = {0};
603 event->GetPrimaryVertex()->GetXYZ(vertex);
604 // TLorentzvector with cluster
605 TLorentzVector clusterVector;
606 cluster->GetMomentum(clusterVector,vertex);
607 Double_t etaCluster = clusterVector.Eta();
608 Double_t phiCluster = clusterVector.Phi();
610 // Histos before cuts
611 if(fHistClusterEtavsPhiBeforeAcc) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster);
613 // Cluster Selection - 0= accept any calo cluster
614 if (fClusterType > 0){
615 //Select EMCAL cluster
616 if (fClusterType == 1 && !cluster->IsEMCAL()){
617 FillClusterCutIndex(kDetector);
619 //Select PHOS cluster
620 if (fClusterType == 2 && !cluster->IsPHOS()){
621 FillClusterCutIndex(kDetector);
626 if(!AcceptanceCuts(cluster,event)){
627 FillClusterCutIndex(kAcceptance);
630 // Cluster Quality Cuts
631 if(!ClusterQualityCuts(cluster,event,isMC)){
632 FillClusterCutIndex(kClusterQuality);
636 // Photon passed cuts
637 FillClusterCutIndex(kPhotonOut);
642 ///________________________________________________________________________
643 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event)
645 // Exclude certain areas for photon reconstruction
648 if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
652 Double_t vertex[3] = {0};
653 event->GetPrimaryVertex()->GetXYZ(vertex);
654 // TLorentzvector with cluster
655 TLorentzVector clusterVector;
656 cluster->GetMomentum(clusterVector,vertex);
657 Double_t etaCluster = clusterVector.Eta();
658 Double_t phiCluster = clusterVector.Phi();
661 if (fMinEtaCut !=-10 && fMaxEtaCut !=10 ){
662 if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut){
663 if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
670 if (fMinPhiCut !=-10000 && fMaxPhiCut !=-10000 ){
671 if (phiCluster < fMinPhiCut || phiCluster > fMaxEtaCut){
672 if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
678 // check distance to bad channel
679 if (cluster->GetDistanceToBadChannel() < fMinDistanceToBadChannel){
680 if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
684 if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
687 if(fHistClusterEtavsPhiAfterAcc) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster);
692 ///________________________________________________________________________
693 Bool_t AliCaloPhotonCuts::UpdateCutString() {
694 ///Update the cut string (if it has been created yet)
696 if(fCutString && fCutString->GetString().Length() == kNCuts) {
697 fCutString->SetString(GetCutNumber());
704 ///________________________________________________________________________
705 Bool_t AliCaloPhotonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
706 // Initialize Cuts from a given Cut string
707 AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
708 if(analysisCutSelection.Length()!=kNCuts) {
709 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
712 if(!analysisCutSelection.IsDigit()){
713 AliError("Cut selection contains characters");
717 const char *cutSelection = analysisCutSelection.Data();
718 #define ASSIGNARRAY(i) fCuts[i] = cutSelection[i] - '0'
719 for(Int_t ii=0;ii<kNCuts;ii++){
723 // Set Individual Cuts
724 for(Int_t ii=0;ii<kNCuts;ii++){
725 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
727 PrintCutsWithValues();
731 ///________________________________________________________________________
732 Bool_t AliCaloPhotonCuts::SetCut(cutIds cutID, const Int_t value) {
733 ///Set individual cut ID
738 if( SetClusterTypeCut(value)) {
739 fCuts[kClusterType] = value;
742 } else return kFALSE;
745 if( SetMinEtaCut(value)) {
746 fCuts[kEtaMin] = value;
749 } else return kFALSE;
752 if( SetMaxEtaCut(value)) {
753 fCuts[kEtaMax] = value;
756 } else return kFALSE;
759 if( SetMinPhiCut(value)) {
760 fCuts[kPhiMin] = value;
763 } else return kFALSE;
766 if( SetMaxPhiCut(value)) {
767 fCuts[kPhiMax] = value;
770 } else return kFALSE;
772 case kDistanceToBadChannel:
773 if( SetDistanceToBadChannelCut(value)) {
774 fCuts[kDistanceToBadChannel] = value;
777 } else return kFALSE;
780 if( SetTimingCut(value)) {
781 fCuts[kTiming] = value;
784 } else return kFALSE;
787 if( SetTrackMatchingCut(value)) {
788 fCuts[kTrackMatching] = value;
791 } else return kFALSE;
794 if( SetExoticCellCut(value)) {
795 fCuts[kExoticCell] = value;
798 } else return kFALSE;
801 if( SetMinEnergyCut(value)) {
802 fCuts[kMinEnery] = value;
805 } else return kFALSE;
808 if( SetMinNCellsCut(value)) {
809 fCuts[kNMinCells] = value;
812 } else return kFALSE;
815 if( SetMinM02(value)) {
816 fCuts[kMinM02] = value;
819 } else return kFALSE;
822 if( SetMaxM02(value)) {
823 fCuts[kMaxM02] = value;
826 } else return kFALSE;
829 if( SetMinM20(value)) {
830 fCuts[kMinM20] = value;
833 } else return kFALSE;
836 if( SetMaxM20(value)) {
837 fCuts[kMaxM20] = value;
840 } else return kFALSE;
843 if( SetDispersion(value)) {
844 fCuts[kDispersion] = value;
847 } else return kFALSE;
854 } else return kFALSE;
857 AliError("Cut id out of range");
861 AliError("Cut id %d not recognized");
866 ///________________________________________________________________________
867 void AliCaloPhotonCuts::PrintCuts() {
868 // Print out current Cut Selection
869 for(Int_t ic = 0; ic < kNCuts; ic++) {
870 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
874 void AliCaloPhotonCuts::PrintCutsWithValues() {
875 // Print out current Cut Selection with value
876 printf("Acceptance cuts: \n");
877 printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
878 printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
879 printf("\tmin distance to bad channel > %3.2f\n", fMinDistanceToBadChannel );
881 printf("Cluster Quality cuts: \n");
882 printf("\t time difference < %3.2f\n", fMaxTimeDiff );
883 printf("\t min distance to track > %3.2f\n", fMinDistTrackToCluster );
884 printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
885 printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
889 // EMCAL acceptance 2011
890 // 1.39626, 3.125 (phi)
894 ///________________________________________________________________________
895 Bool_t AliCaloPhotonCuts::SetClusterTypeCut(Int_t clusterType)
898 case 0: // all clusters
901 case 1: // EMCAL clusters
904 case 2: // PHOS clusters
908 AliError(Form("ClusterTypeCut not defined %d",clusterType));
914 //___________________________________________________________________
915 Bool_t AliCaloPhotonCuts::SetMinEtaCut(Int_t minEta)
928 AliError(Form("MinEta Cut not defined %d",minEta));
935 //___________________________________________________________________
936 Bool_t AliCaloPhotonCuts::SetMaxEtaCut(Int_t maxEta)
949 AliError(Form("MaxEta Cut not defined %d",maxEta));
955 //___________________________________________________________________
956 Bool_t AliCaloPhotonCuts::SetMinPhiCut(Int_t minPhi)
966 AliError(Form("MinPhi Cut not defined %d",minPhi));
972 //___________________________________________________________________
973 Bool_t AliCaloPhotonCuts::SetMaxPhiCut(Int_t maxPhi)
983 AliError(Form("Max Phi Cut not defined %d",maxPhi));
989 //___________________________________________________________________
990 Bool_t AliCaloPhotonCuts::SetDistanceToBadChannelCut(Int_t distanceToBadChannel)
992 switch(distanceToBadChannel){
994 fMinDistanceToBadChannel=0;
997 fMinDistanceToBadChannel=5;
1000 AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
1006 //___________________________________________________________________
1007 Bool_t AliCaloPhotonCuts::SetTimingCut(Int_t timing)
1014 fMaxTimeDiff=10e-7; //1000ns
1017 fMaxTimeDiff=50e-8; //500ns
1020 fMaxTimeDiff=20e-8; //200ns
1023 fMaxTimeDiff=10e-8; //100ns
1026 fMaxTimeDiff=50e-9; //50ns
1030 AliError(Form("Timing Cut not defined %d",timing));
1036 //___________________________________________________________________
1037 Bool_t AliCaloPhotonCuts::SetTrackMatchingCut(Int_t trackMatching)
1039 switch(trackMatching){
1041 fMinDistTrackToCluster=0;
1044 fMinDistTrackToCluster=5;
1047 AliError(Form("Track Matching Cut not defined %d",trackMatching));
1053 //___________________________________________________________________
1054 Bool_t AliCaloPhotonCuts::SetExoticCellCut(Int_t exoticCell)
1064 AliError(Form("Exotic cell Cut not defined %d",exoticCell));
1070 //___________________________________________________________________
1071 Bool_t AliCaloPhotonCuts::SetMinEnergyCut(Int_t minEnergy)
1105 AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
1111 //___________________________________________________________________
1112 Bool_t AliCaloPhotonCuts::SetMinNCellsCut(Int_t minNCells)
1138 AliError(Form("Min N cells Cut not defined %d",minNCells));
1144 //___________________________________________________________________
1145 Bool_t AliCaloPhotonCuts::SetMaxM02(Int_t maxM02)
1165 AliError(Form("Max M02 Cut not defined %d",maxM02));
1171 //___________________________________________________________________
1172 Bool_t AliCaloPhotonCuts::SetMinM02(Int_t minM02)
1182 AliError(Form("Min M02 not defined %d",minM02));
1188 //___________________________________________________________________
1189 Bool_t AliCaloPhotonCuts::SetMaxM20(Int_t maxM20)
1199 AliError(Form("Max M20 Cut not defined %d",maxM20));
1205 //___________________________________________________________________
1206 Bool_t AliCaloPhotonCuts::SetMinM20(Int_t minM20)
1216 AliError(Form("Min M20 Cut not defined %d",minM20));
1222 //___________________________________________________________________
1223 Bool_t AliCaloPhotonCuts::SetDispersion(Int_t dispersion)
1227 fMaxDispersion =100;
1233 AliError(Form("Maximum Dispersion Cut not defined %d",dispersion));
1239 //___________________________________________________________________
1240 Bool_t AliCaloPhotonCuts::SetNLM(Int_t nlm)
1252 AliError(Form("NLM Cut not defined %d",nlm));
1258 ///________________________________________________________________________
1259 TString AliCaloPhotonCuts::GetCutNumber(){
1260 // returns TString with current cut number
1262 for(Int_t ii=0;ii<kNCuts;ii++){
1263 a.Append(Form("%d",fCuts[ii]));