]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/GammaConv/AliCaloPhotonCuts.cxx
- changes in AliAnalysisTaskGammaConvCalo for debugging
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliCaloPhotonCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                                                                *
4  * Authors: Friederike Bock, Baldo Sahlmueller                                *
5  * Version 1.0                                                                                    *
6  *                                                                                                                *
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  **************************************************************************/
15
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class handling all kinds of selection cuts for
19 // Photon from EMCAL clusters
20 //---------------------------------------------
21 ////////////////////////////////////////////////
22
23 #include "AliCaloPhotonCuts.h"
24 #include "AliAnalysisManager.h"
25 #include "AliInputEventHandler.h"
26 #include "AliMCEventHandler.h"
27 #include "AliAODHandler.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TF1.h"
31 #include "AliStack.h"
32 #include "AliAODConversionMother.h"
33 #include "TObjString.h"
34 #include "AliAODEvent.h"
35 #include "AliESDEvent.h"
36 #include "AliCentrality.h"
37 #include "TList.h"
38 #include "TFile.h"
39 #include "AliLog.h"
40 #include "AliV0ReaderV1.h"
41 #include "AliAODMCParticle.h"
42 #include "AliAODMCHeader.h"
43 #include "AliPicoTrack.h"
44
45 class iostream;
46
47 using namespace std;
48
49 ClassImp(AliCaloPhotonCuts)
50
51
52 const char* AliCaloPhotonCuts::fgkCutNames[AliCaloPhotonCuts::kNCuts] = {
53         "ClusterType",          //0   0: all,    1: EMCAL,   2: PHOS
54         "EtaMin",               //1   0: -10,    1: -0.6687, 2: -0,5, 3: -2
55         "EtaMax",               //2   0: 10,     1: 0.66465, 2: 0.5,  3: 2
56         "PhiMin",               //3   0: -10000, 1: 1.39626
57         "PhiMax",               //4   0: 10000, 1: 3.125
58         "DistanceToBadChannel", //5   0: 0,      1: 5
59         "Timing",               //6   0: no cut
60         "TrackMatching",        //7   0: 0,      1: 5
61         "ExoticCell",           //8   0: no cut
62         "MinEnergy",            //9   0: no cut, 1: 0.05,    2: 0.1,  3: 0.15, 4: 0.2, 5: 0.3, 6: 0.5, 7: 0.75, 8: 1, 9: 1.25 (all GeV)
63         "MinNCells",            //10  0: no cut, 1: 1,       2: 2,    3: 3,    4: 4,   5: 5,   6: 6
64         "MinM02",                               //11
65         "MaxM02",                               //12
66         "MinM20",                               //13
67         "MaxM20",                               //14
68         "MaximumDispersion",    //15
69         "NLM"                                   //16
70 };
71
72
73 //________________________________________________________________________
74 AliCaloPhotonCuts::AliCaloPhotonCuts(const char *name,const char *title) :
75         AliAnalysisCuts(name,title),
76         fHistograms(NULL),      
77         fClusterType(0),
78         fMinEtaCut(-10),
79         fMaxEtaCut(10),
80         fUseEtaCut(0),
81         fMinPhiCut(-10000),
82         fMaxPhiCut(-10000),
83         fUsePhiCut(0),
84         fMinDistanceToBadChannel(0),
85         fUseDistanceToBadChannel(0),
86         fMaxTimeDiff(10e10),
87         fUseTimeDiff(0),
88         fMinDistTrackToCluster(0),
89         fUseDistTrackToCluster(0),
90         fExoticCell(0),
91         fUseExoticCell(0),
92         fMinEnergy(0),
93         fUseMinEnergy(0),
94         fMinNCells(0),
95         fUseNCells(0),
96         fMaxM02(1000),
97         fMinM02(0),
98         fUseM02(0),
99         fMaxM20(1000),
100         fMinM20(0),
101         fUseM20(0),
102         fMaxDispersion(1000),
103         fUseDispersion(0),
104         fMinNLM(0),
105         fMaxNLM(1000),
106         fUseNLM(0),
107         fCutString(NULL),
108         fHistCutIndex(NULL),
109         fHistAcceptanceCuts(NULL),
110         fHistClusterIdentificationCuts(NULL),
111         fHistClusterEtavsPhiBeforeAcc(NULL),
112         fHistClusterEtavsPhiAfterAcc(NULL),
113         fHistClusterEtavsPhiAfterQA(NULL),
114         fHistDistanceToBadChannelBeforeAcc(NULL),
115         fHistDistanceToBadChannelAfterAcc(NULL),
116         fHistClusterTimevsEBeforeQA(NULL),
117         fHistClusterTimevsEAfterQA(NULL),
118         fHistExoticCellBeforeQA(NULL),
119         fHistExoticCellAfterQA(NULL),
120         fHistNMatchedTracks(NULL),
121         fHistDistanceTrackToClusterBeforeQA(NULL),
122         fHistDistanceTrackToClusterAfterQA(NULL),
123         fHistEnergyOfClusterBeforeQA(NULL),
124         fHistEnergyOfClusterAfterQA(NULL),
125         fHistNCellsBeforeQA(NULL),
126         fHistNCellsAfterQA(NULL),
127         fHistM02BeforeQA(NULL),
128         fHistM02AfterQA(NULL),
129         fHistM20BeforeQA(NULL),
130         fHistM20AfterQA(NULL),
131         fHistDispersionBeforeQA(NULL),
132         fHistDispersionAfterQA(NULL),
133         fHistNLMBeforeQA(NULL),
134         fHistNLMAfterQA(NULL)
135 {
136    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
137    fCutString=new TObjString((GetCutNumber()).Data());
138 }
139
140 //________________________________________________________________________
141 AliCaloPhotonCuts::AliCaloPhotonCuts(const AliCaloPhotonCuts &ref) :
142    AliAnalysisCuts(ref),
143         fHistograms(NULL),      
144         fClusterType(ref.fClusterType),
145         fMinEtaCut(ref.fMinEtaCut),
146         fMaxEtaCut(ref.fMaxEtaCut),
147         fUseEtaCut(ref.fUseEtaCut),
148         fMinPhiCut(ref.fMinPhiCut),
149         fMaxPhiCut(ref.fMaxPhiCut),
150         fUsePhiCut(ref.fUsePhiCut),
151         fMinDistanceToBadChannel(ref.fMinDistanceToBadChannel),
152         fUseDistanceToBadChannel(ref.fUseDistanceToBadChannel),
153         fMaxTimeDiff(ref.fMaxTimeDiff),
154         fUseTimeDiff(ref.fUseTimeDiff),
155         fMinDistTrackToCluster(ref.fMinDistTrackToCluster),
156         fUseDistTrackToCluster(ref.fUseDistTrackToCluster),
157         fExoticCell(ref.fExoticCell),
158         fUseExoticCell(ref.fUseExoticCell),
159         fMinEnergy(ref.fMinEnergy),
160         fUseMinEnergy(ref.fUseMinEnergy),
161         fMinNCells(ref.fMinNCells),
162         fUseNCells(ref.fUseNCells),
163         fMaxM02(ref.fMaxM02),
164         fMinM02(ref.fMinM02),
165         fUseM02(ref.fUseM02),
166         fMaxM20(ref.fMaxM20),
167         fMinM20(ref.fMinM20),
168         fUseM20(ref.fUseDispersion),
169         fMaxDispersion(ref.fMaxDispersion),
170         fUseDispersion(ref.fUseDispersion),
171         fMinNLM(ref.fMinNLM),
172         fMaxNLM(ref.fMaxNLM),
173         fUseNLM(ref.fUseNLM),
174         fCutString(NULL),
175         fHistCutIndex(NULL),
176         fHistAcceptanceCuts(NULL),
177         fHistClusterIdentificationCuts(NULL),
178         fHistClusterEtavsPhiBeforeAcc(NULL),
179         fHistClusterEtavsPhiAfterAcc(NULL),
180         fHistClusterEtavsPhiAfterQA(NULL),
181         fHistDistanceToBadChannelBeforeAcc(NULL),
182         fHistDistanceToBadChannelAfterAcc(NULL),
183         fHistClusterTimevsEBeforeQA(NULL),
184         fHistClusterTimevsEAfterQA(NULL),
185         fHistExoticCellBeforeQA(NULL),
186         fHistExoticCellAfterQA(NULL),
187         fHistNMatchedTracks(NULL),
188         fHistDistanceTrackToClusterBeforeQA(NULL),
189         fHistDistanceTrackToClusterAfterQA(NULL),
190         fHistEnergyOfClusterBeforeQA(NULL),
191         fHistEnergyOfClusterAfterQA(NULL),
192         fHistNCellsBeforeQA(NULL),
193         fHistNCellsAfterQA(NULL),
194         fHistM02BeforeQA(NULL),
195         fHistM02AfterQA(NULL),
196         fHistM20BeforeQA(NULL),
197         fHistM20AfterQA(NULL),
198         fHistDispersionBeforeQA(NULL),
199         fHistDispersionAfterQA(NULL),
200         fHistNLMBeforeQA(NULL),
201         fHistNLMAfterQA(NULL)
202 {
203    // Copy Constructor
204    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
205    fCutString=new TObjString((GetCutNumber()).Data());
206
207 }
208
209
210 //________________________________________________________________________
211 AliCaloPhotonCuts::~AliCaloPhotonCuts() {
212    // Destructor
213    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
214    // if(fHistograms)
215    //    delete fHistograms;
216    // fHistograms = NULL;
217    if(fCutString != NULL){
218       delete fCutString;
219       fCutString = NULL;
220    }
221 }
222
223 //________________________________________________________________________
224 void AliCaloPhotonCuts::InitCutHistograms(TString name){
225
226         // Initialize Cut Histograms for QA (only initialized and filled if function is called)
227         TH1::AddDirectory(kFALSE);
228
229         if(fHistograms != NULL){
230                 delete fHistograms;
231                 fHistograms=NULL;
232         }
233         if(fHistograms==NULL){
234                 fHistograms=new TList();
235                 fHistograms->SetOwner(kTRUE);
236                 if(name=="")fHistograms->SetName(Form("CaloCuts_%s",GetCutNumber().Data()));
237                 else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
238         }
239
240         // IsPhotonSelected
241         fHistCutIndex=new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",5,-0.5,4.5);
242         fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in");
243         fHistCutIndex->GetXaxis()->SetBinLabel(kDetector+1,"detector");
244         fHistCutIndex->GetXaxis()->SetBinLabel(kAcceptance+1,"acceptance");
245         fHistCutIndex->GetXaxis()->SetBinLabel(kClusterQuality+1,"cluster QA");
246         fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonOut+1,"out");
247         fHistograms->Add(fHistCutIndex);
248
249         // Acceptance Cuts
250         fHistAcceptanceCuts=new TH1F(Form("AcceptanceCuts %s",GetCutNumber().Data()),"AcceptanceCuts",5,-0.5,4.5);
251         fHistAcceptanceCuts->GetXaxis()->SetBinLabel(1,"in");
252         fHistAcceptanceCuts->GetXaxis()->SetBinLabel(2,"eta");
253         fHistAcceptanceCuts->GetXaxis()->SetBinLabel(3,"phi");
254         fHistAcceptanceCuts->GetXaxis()->SetBinLabel(4,"distance to bad channel");
255         fHistAcceptanceCuts->GetXaxis()->SetBinLabel(5,"out");
256         fHistograms->Add(fHistAcceptanceCuts);
257
258         // Cluster Cuts
259         fHistClusterIdentificationCuts =new TH1F(Form("ClusterQualityCuts %s",GetCutNumber().Data()),"ClusterQualityCuts",11,-0.5,10.5);
260         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(1,"in");
261         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(2,"timing");
262         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(3,"track matching");
263         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(4,"Exotics");
264         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(5,"minimum energy");
265         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(6,"minimum NCells");
266         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(7,"M02");
267         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(8,"M20");
268         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(9,"dispersion");
269         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(10,"NLM");
270         fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(11,"out");
271         fHistograms->Add(fHistClusterIdentificationCuts);
272
273         // Acceptance related histogramms
274         fHistClusterEtavsPhiBeforeAcc=new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
275         fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
276         fHistClusterEtavsPhiAfterAcc=new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
277         fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
278         fHistClusterEtavsPhiAfterQA=new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",462,-TMath::Pi(),TMath::Pi(),110,-0.7,0.7);
279         fHistograms->Add(fHistClusterEtavsPhiAfterQA);
280         fHistDistanceToBadChannelBeforeAcc = new TH1F(Form("DistanceToBadChannel_beforeAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_beforeAcceptance",200,0,40);
281         fHistograms->Add(fHistDistanceToBadChannelBeforeAcc);
282         fHistDistanceToBadChannelAfterAcc = new TH1F(Form("DistanceToBadChannel_afterAcceptance %s",GetCutNumber().Data()),"DistanceToBadChannel_afterAcceptance",200,0,40);
283         fHistograms->Add(fHistDistanceToBadChannelAfterAcc);
284         
285         // Cluster quality related histograms
286         fHistClusterTimevsEBeforeQA=new TH2F(Form("ClusterTimeVsE_beforeClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_beforeClusterQA",400,-10e-6,10e-6,100,0.,40);
287         fHistograms->Add(fHistClusterTimevsEBeforeQA);
288         fHistClusterTimevsEAfterQA=new TH2F(Form("ClusterTimeVsE_afterClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_afterClusterQA",400,-10e-6,10e-6,100,0.,40);
289         fHistograms->Add(fHistClusterTimevsEAfterQA);
290         fHistExoticCellBeforeQA=new TH2F(Form("ExoticCell_beforeClusterQA %s",GetCutNumber().Data()),"ExoticCell_beforeClusterQA",400,0,40,50,0.75,1);
291         fHistograms->Add(fHistExoticCellBeforeQA);
292         fHistExoticCellAfterQA=new TH2F(Form("ExoticCell_afterClusterQA %s",GetCutNumber().Data()),"ExoticCell_afterClusterQA",400,0,40,50,0.75,1);
293         fHistograms->Add(fHistExoticCellAfterQA);
294         fHistNMatchedTracks = new TH1F(Form("NMatchedTracks_%s",GetCutNumber().Data()),"NMatchedTracks",22,-1.5,20.5);
295         fHistograms->Add(fHistNMatchedTracks);
296         fHistDistanceTrackToClusterBeforeQA = new TH1F(Form("DistanceToTrack_beforeClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_beforeClusterQA",200,0,2);
297         fHistograms->Add(fHistDistanceTrackToClusterBeforeQA);
298         fHistDistanceTrackToClusterAfterQA = new TH1F(Form("DistanceToTrack_afterClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_afterClusterQA",200,0,2);
299         fHistograms->Add(fHistDistanceTrackToClusterAfterQA);
300         fHistEnergyOfClusterBeforeQA = new TH1F(Form("EnergyOfCluster_beforeClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_beforeClusterQA",300,0,30);
301         fHistograms->Add(fHistEnergyOfClusterBeforeQA);
302         fHistEnergyOfClusterAfterQA = new TH1F(Form("EnergyOfCluster_afterClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_afterClusterQA",300,0,30);
303         fHistograms->Add(fHistEnergyOfClusterAfterQA);
304         fHistNCellsBeforeQA = new TH1F(Form("NCellPerCluster_beforeClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_beforeClusterQA",50,0,50);
305         fHistograms->Add(fHistNCellsBeforeQA);
306         fHistNCellsAfterQA = new TH1F(Form("NCellPerCluster_afterClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_afterClusterQA",50,0,50);
307         fHistograms->Add(fHistNCellsAfterQA);
308         fHistM02BeforeQA = new TH1F(Form("M02_beforeClusterQA %s",GetCutNumber().Data()),"M02_beforeClusterQA",100,0,5);
309         fHistograms->Add(fHistM02BeforeQA);
310         fHistM02AfterQA = new TH1F(Form("M02_afterClusterQA %s",GetCutNumber().Data()),"M02_afterClusterQA",100,0,5);
311         fHistograms->Add(fHistM02AfterQA);
312         fHistM20BeforeQA = new TH1F(Form("M20_beforeClusterQA %s",GetCutNumber().Data()),"M20_beforeClusterQA",100,0,2.5);
313         fHistograms->Add(fHistM20BeforeQA);
314         fHistM20AfterQA = new TH1F(Form("M20_afterClusterQA %s",GetCutNumber().Data()),"M20_afterClusterQA",100,0,2.5);
315         fHistograms->Add(fHistM20AfterQA);
316         fHistDispersionBeforeQA = new TH1F(Form("Dispersion_beforeClusterQA %s",GetCutNumber().Data()),"Dispersion_beforeClusterQA",100,0,4);
317         fHistograms->Add(fHistDispersionBeforeQA);
318         fHistDispersionAfterQA = new TH1F(Form("Dispersion_afterClusterQA %s",GetCutNumber().Data()),"Dispersion_afterClusterQA",100,0,4);
319         fHistograms->Add(fHistDispersionAfterQA);
320         fHistNLMBeforeQA = new TH1F(Form("NLM_beforeClusterQA %s",GetCutNumber().Data()),"NLM_beforeClusterQA",10,0,10);
321         fHistograms->Add(fHistNLMBeforeQA);
322         fHistNLMAfterQA = new TH1F(Form("NLM_afterClusterQA %s",GetCutNumber().Data()),"NLM_afterClusterQA",10,0,10);
323         fHistograms->Add(fHistNLMAfterQA);
324         
325         TH1::AddDirectory(kTRUE);
326 }
327
328 /*
329 ///________________________________________________________________________
330 Bool_t AliCaloPhotonCuts::ClusterIsSelectedMC(TParticle *particle,AliStack *fMCStack,Bool_t checkForConvertedGamma){
331    // MonteCarlo Photon Selection
332
333    if(!fMCStack)return kFALSE;
334
335    if (particle->GetPdgCode() == 22){
336
337
338      if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) )
339          return kFALSE;
340       if(fEtaCutMin>-0.1){
341          if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) )
342             return kFALSE;
343       }
344
345       if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
346          return kFALSE; // no photon as mothers!
347       }
348
349       if(particle->GetMother(0) >= fMCStack->GetNprimary()){
350          return kFALSE; // the gamma has a mother, and it is not a primary particle
351       }
352
353       if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
354
355       // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
356       TParticle* ePos = NULL;
357       TParticle* eNeg = NULL;
358
359       if(particle->GetNDaughters() >= 2){
360          for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
361             TParticle *tmpDaughter = fMCStack->Particle(daughterIndex);
362             if(tmpDaughter->GetUniqueID() == 5){
363                if(tmpDaughter->GetPdgCode() == 11){
364                   eNeg = tmpDaughter;
365                } else if(tmpDaughter->GetPdgCode() == -11){
366                   ePos = tmpDaughter;
367                }
368             }
369          }
370       }
371
372       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
373          return kFALSE;
374       }
375
376       if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
377          return kFALSE; // no reconstruction below the Pt cut
378       }
379
380       if( ePos->Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) ||
381           eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) )
382          return kFALSE;
383
384       if(fEtaCutMin > -0.1){
385          if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) ||
386              (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) )
387             return kFALSE;
388       }
389
390       if(ePos->R()>fMaxR){
391          return kFALSE; // cuts on distance from collision point
392       }
393
394       if(abs(ePos->Vz()) > fMaxZ){
395          return kFALSE;  // outside material
396       }
397       if(abs(eNeg->Vz()) > fMaxZ){
398          return kFALSE;  // outside material
399       }
400
401       if( ePos->R() <= ((abs(ePos->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
402          return kFALSE;  // line cut to exclude regions where we do not reconstruct
403       } else if ( fEtaCutMin != -0.1 &&   ePos->R() >= ((abs(ePos->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
404          return kFALSE;
405       }
406
407       if( eNeg->R() <= ((abs(eNeg->Vz()) * fLineCutZRSlope) - fLineCutZValue)){
408          return kFALSE; // line cut to exclude regions where we do not reconstruct
409       } else if ( fEtaCutMin != -0.1 &&   eNeg->R() >= ((abs(eNeg->Vz()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
410          return kFALSE;
411       }
412
413       return kTRUE;
414       //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
415    }
416    return kFALSE;
417 }
418 ///________________________________________________________________________
419 Bool_t AliCaloPhotonCuts::ClusterIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray,Bool_t checkForConvertedGamma){
420    // MonteCarlo Photon Selection
421
422    if(!aodmcArray)return kFALSE;
423
424    if (particle->GetPdgCode() == 22){
425       if( particle->Eta() > (fEtaCut) || particle->Eta() < (-fEtaCut) )
426          return kFALSE;
427       if(fEtaCutMin>-0.1){
428          if( particle->Eta() < (fEtaCutMin) && particle->Eta() > (-fEtaCutMin) )
429             return kFALSE;
430       }
431
432       if(particle->GetMother() > -1){
433          if((static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){
434             return kFALSE; // no photon as mothers!
435          }
436          if(!(static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother()))->IsPrimary())){
437             return kFALSE; // the gamma has a mother, and it is not a primary particle
438          }
439       }
440
441       if(!checkForConvertedGamma) return kTRUE; // return in case of accepted gamma
442
443       // looking for conversion gammas (electron + positron from pairbuilding (= 5) )
444       AliAODMCParticle* ePos = NULL;
445       AliAODMCParticle* eNeg = NULL;
446
447       if(particle->GetNDaughters() >= 2){
448          for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
449             AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(aodmcArray->At(daughterIndex));
450             if(!tmpDaughter) continue;
451             if(((tmpDaughter->GetMCProcessCode())) == 5){    // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
452                if(tmpDaughter->GetPdgCode() == 11){
453                   eNeg = tmpDaughter;
454                } else if(tmpDaughter->GetPdgCode() == -11){
455                   ePos = tmpDaughter;
456                }
457             }
458          }
459       }
460
461       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
462          return kFALSE;
463       }
464
465       if(ePos->Pt()<fSinglePtCut || eNeg->Pt()<fSinglePtCut){
466          return kFALSE; // no reconstruction below the Pt cut
467       }
468
469       if( ePos->Eta() > (fEtaCut) || ePos->Eta() < (-fEtaCut) ||
470           eNeg->Eta() > (fEtaCut) || eNeg->Eta() < (-fEtaCut) )
471          return kFALSE;
472
473       if(fEtaCutMin > -0.1){
474          if( (ePos->Eta() < (fEtaCutMin) && ePos->Eta() > (-fEtaCutMin)) ||
475              (eNeg->Eta() < (fEtaCutMin) && eNeg->Eta() > (-fEtaCutMin)) )
476             return kFALSE;
477       }
478
479       Double_t rPos = sqrt( (ePos->Xv()*ePos->Xv()) + (ePos->Yv()*ePos->Yv()) );
480       Double_t rNeg = sqrt( (eNeg->Xv()*eNeg->Xv()) + (eNeg->Yv()*eNeg->Yv()) );
481
482       if(rPos>fMaxR){
483          return kFALSE; // cuts on distance from collision point
484       }
485       if(abs(ePos->Zv()) > fMaxZ){
486          return kFALSE;  // outside material
487       }
488       if(abs(eNeg->Zv()) > fMaxZ){
489          return kFALSE;  // outside material
490       }
491
492       if( rPos <= ((abs(ePos->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
493          return kFALSE;  // line cut to exclude regions where we do not reconstruct
494       } else if ( fEtaCutMin != -0.1 &&   rPos >= ((abs(ePos->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
495          return kFALSE;
496       }
497
498       if( rNeg <= ((abs(eNeg->Zv()) * fLineCutZRSlope) - fLineCutZValue)){
499          return kFALSE; // line cut to exclude regions where we do not reconstruct
500       } else if ( fEtaCutMin != -0.1 &&   rNeg >= ((abs(eNeg->Zv()) * fLineCutZRSlopeMin) - fLineCutZValueMin)){
501          return kFALSE;
502       }
503
504       return kTRUE;
505       //if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
506    }
507    return kFALSE;
508 }*/
509
510
511
512 ///________________________________________________________________________
513 // This function selects the clusters based on their quality criteria
514 ///________________________________________________________________________
515 Bool_t AliCaloPhotonCuts::ClusterQualityCuts(AliVCluster* cluster, AliVEvent *event, Bool_t isMC)
516 {   // Specific Photon Cuts
517
518         Int_t cutIndex = 0;
519         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex);
520         cutIndex++;
521
522         Double_t minR = 999.0;
523         // get the minimum radius of tracks to cluster
524         if(fHistDistanceTrackToClusterBeforeQA || fHistDistanceTrackToClusterAfterQA){
525 //              Float_t pos[3];
526 //              cluster->GetPosition(pos);  // Get cluster position
527 //              
528 //              TVector3 cp(pos);
529 //              int NtrMatched = 0;
530 //              NtrMatched = cluster->GetNTracksMatched();
531 //              fHistNMatchedTracks->Fill(NtrMatched);
532 //              //loop over tracks for Jet QA
533 //              TList *l = event->GetList();
534 //              TClonesArray *tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
535 //              for(int itrack = 0; itrack < NtrMatched; itrack++){
536 //                      AliVTrack *trackcluster = static_cast<AliVTrack*>(tracks->At(itrack));
537 //                      if (! trackcluster) {
538 //                              AliError(Form("Couldn't get ESD track %d\n", itrack));
539 //                              continue;
540 //                      }
541 //                      Double_t dphi = -999.0;
542 //                      Double_t deta = -999.0;
543 //                      AliPicoTrack::GetEtaPhiDiff(trackcluster, cluster, dphi, deta);
544 //                      cout << "here" << endl;
545 //                      Double_t dr = sqrt(dphi*dphi + deta+deta);
546 //                      if(dr < minR)
547 //                              minR = dr;
548 //              }//loop over tracks
549         }
550         
551         // Fill Histos before Cuts
552         if(fHistClusterTimevsEBeforeQA) fHistClusterTimevsEBeforeQA->Fill(cluster->GetTOF(), cluster->E());
553 //      if(fHistExoticCellBeforeQA) fHistExoticCellBeforeQA->Fill(cluster->E(), );
554         if(fHistDistanceTrackToClusterBeforeQA) fHistDistanceTrackToClusterBeforeQA->Fill(minR);
555         if(fHistEnergyOfClusterBeforeQA) fHistEnergyOfClusterBeforeQA->Fill(cluster->E());
556         if(fHistNCellsBeforeQA) fHistNCellsBeforeQA->Fill(cluster->GetNCells());
557         if(fHistM02BeforeQA) fHistM02BeforeQA->Fill(cluster->GetM02());
558         if(fHistM20BeforeQA) fHistM20BeforeQA->Fill(cluster->GetM20());
559         if(fHistDispersionBeforeQA) fHistDispersionBeforeQA->Fill(cluster->GetDispersion());
560 //      if(fHistNLMBeforeQA) fHistNLMBeforeQA->Fill(cluster->GetNExMax());
561         
562         // Check wether timing is ok
563         if (fUseTimeDiff){
564                 if(abs(cluster->GetTOF()) > fMaxTimeDiff && !isMC){
565                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //1
566                         return kFALSE;
567                 }
568         }       
569         cutIndex++; //2, next cut
570
571         // Minimum distance to track
572         if (fUseDistTrackToCluster){
573                 Float_t pos[3];
574                 cluster->GetPosition(pos);  // Get cluster position
575                 TVector3 cp(pos);
576                 int NtrMatched = 0;
577                 NtrMatched = cluster->GetNTracksMatched();
578                 fHistNMatchedTracks->Fill(NtrMatched);
579                 if(NtrMatched > 0){
580                         //loop over tracks for QA
581                         TList *l = event->GetList();
582                         TClonesArray *tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
583                         
584                         Double_t dphi = 999.0;
585                         Double_t deta = 999.0;
586                         Double_t dr2 = 999.0;
587
588                         for(int itrack = 0; itrack < NtrMatched; itrack++){
589                                 AliVTrack *trackcluster = NULL;
590                                 trackcluster = static_cast<AliVTrack*>(tracks->At(itrack));
591                                 if (! trackcluster) {
592                                 AliError(Form("Couldn't get ESD track %d\n", itrack));
593                                 continue;
594                                 }
595                                 AliPicoTrack::GetEtaPhiDiff(trackcluster, cluster, dphi, deta);
596                                 dr2 = dphi*dphi + deta+deta;
597                                 //cout << dr << endl;
598                                 if(dr2 < fMinDistTrackToCluster*fMinDistTrackToCluster){
599                 //        if(dphi < fMinDistTrackToCluster || deta < fMinDistTrackToCluster){
600                                 if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2
601                                 return kFALSE;
602                                 }
603                                 
604                         }//loop over tracks
605                 }
606 //              if(cluster->GetEmcCpvDistance() < fMinDistTrackToCluster){
607 //                      if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2
608 //                      return kFALSE;
609 //              }
610         }
611         cutIndex++;//3, next cut
612
613         // exotic cell cut --IMPLEMENT LATER---
614 //      if(!AcceptanceCuts(photon)){
615 //              if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //3
616 //              return kFALSE;
617 //      }
618         cutIndex++; //4, next cut
619         
620         // minimum cell energy cut
621         if (fUseMinEnergy){
622                 if(cluster->E() < fMinEnergy){
623                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //4
624                         return kFALSE;
625                 }
626         }       
627         cutIndex++; //5, next cut
628         
629         // minimum number of cells
630         if (fUseNCells){
631                 if(cluster->GetNCells() < fMinNCells) {
632                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //5
633                         return kFALSE;
634                 }
635         }       
636         cutIndex++; //6, next cut
637         
638         // M02 cut
639         if (fUseM02){
640                 if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 ) {
641                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //6
642                         return kFALSE;
643                 }
644         }       
645         cutIndex++; //7, next cut
646         
647         // M20 cut
648         if (fUseM20){
649                 if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 ) {
650                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //7
651                         return kFALSE;
652                 }
653         }       
654         cutIndex++; //8, next cut
655         
656         // dispersion cut
657         if (fUseDispersion){
658                 if( cluster->GetDispersion()> fMaxDispersion) {
659                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //8
660                         return kFALSE;
661                 }
662         }       
663         cutIndex++; //9, next cut
664         
665         // NLM cut --IMPLEMENT LATER---
666 //      if (fUseNLM){
667 //              if( cluster->GetDispersion()> fMaxDispersion) {
668 //                      if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //9
669 //                      return kFALSE;
670 //              }
671 //      }       
672         cutIndex++; //9, next cut
673         
674         // DONE with selecting photons
675         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //10
676
677         // Histos after Cuts
678         Double_t vertex[3] = {0};
679         event->GetPrimaryVertex()->GetXYZ(vertex);
680         // TLorentzvector with cluster
681         TLorentzVector clusterVector;
682         cluster->GetMomentum(clusterVector,vertex);
683         Double_t etaCluster = clusterVector.Eta();
684         Double_t phiCluster = clusterVector.Phi();
685         
686         if(fHistClusterEtavsPhiAfterQA) fHistClusterEtavsPhiAfterQA->Fill(phiCluster,etaCluster);
687         if(fHistClusterTimevsEAfterQA) fHistClusterTimevsEAfterQA->Fill(cluster->GetTOF(), cluster->E());
688 //      if(fHistExoticCellAfterQA) fHistExoticCellAfterQA->Fill(cluster->E(), );
689         if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(minR);
690         if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(cluster->GetEmcCpvDistance());
691         if(fHistEnergyOfClusterAfterQA) fHistEnergyOfClusterAfterQA->Fill(cluster->E());
692         if(fHistNCellsAfterQA) fHistNCellsAfterQA->Fill(cluster->GetNCells());
693         if(fHistM02AfterQA) fHistM02AfterQA->Fill(cluster->GetM02());
694         if(fHistM20AfterQA) fHistM20AfterQA->Fill(cluster->GetM20());
695         if(fHistDispersionAfterQA) fHistDispersionAfterQA->Fill(cluster->GetDispersion());
696 //      if(fHistNLMBeforeQA) fHistNLMAfterQA->Fill(cluster->GetNExMax());
697
698         return kTRUE;
699
700 }
701
702
703 ///________________________________________________________________________
704 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, Bool_t isMC)
705 {
706         //Selection of Reconstructed photon clusters with Calorimeters
707
708         FillClusterCutIndex(kPhotonIn);
709
710         Double_t vertex[3] = {0};
711         event->GetPrimaryVertex()->GetXYZ(vertex);
712         // TLorentzvector with cluster
713         TLorentzVector clusterVector;
714         cluster->GetMomentum(clusterVector,vertex);
715         Double_t etaCluster = clusterVector.Eta();
716         Double_t phiCluster = clusterVector.Phi();
717
718         // Histos before cuts
719         if(fHistClusterEtavsPhiBeforeAcc) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster);
720         
721         // Cluster Selection - 0= accept any calo cluster
722         if (fClusterType > 0){
723                 //Select EMCAL cluster
724                 if (fClusterType == 1 && !cluster->IsEMCAL()){
725                         FillClusterCutIndex(kDetector);
726                         return kFALSE;
727                 }
728                 //Select PHOS cluster
729                 if (fClusterType == 2 && !cluster->IsPHOS()){
730                         FillClusterCutIndex(kDetector);
731                         return kFALSE;
732                 }
733         }
734         
735         // Acceptance Cuts
736         if(!AcceptanceCuts(cluster,event)){
737                 FillClusterCutIndex(kAcceptance);
738                 return kFALSE;
739         }
740         // Cluster Quality Cuts
741         if(!ClusterQualityCuts(cluster,event,isMC)){
742                 FillClusterCutIndex(kClusterQuality);
743                 return kFALSE;
744         }
745
746         // Photon passed cuts
747         FillClusterCutIndex(kPhotonOut);
748         return kTRUE;
749 }
750
751
752 ///________________________________________________________________________
753 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event) 
754 {
755    // Exclude certain areas for photon reconstruction
756
757         Int_t cutIndex=0;
758         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
759         cutIndex++;
760
761         
762         Double_t vertex[3] = {0};
763         event->GetPrimaryVertex()->GetXYZ(vertex);
764         // TLorentzvector with cluster
765         TLorentzVector clusterVector;
766         cluster->GetMomentum(clusterVector,vertex);
767         Double_t etaCluster = clusterVector.Eta();
768         Double_t phiCluster = clusterVector.Phi();
769         
770         // check eta range
771         if (fUseEtaCut){
772                 if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut){
773                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
774                         return kFALSE;
775                 }
776         }
777         cutIndex++;
778         
779         // check phi range
780         if (fUsePhiCut ){
781                 if (phiCluster < fMinPhiCut || phiCluster > fMaxEtaCut){
782                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
783                         return kFALSE;
784                 }
785         }
786         cutIndex++;
787         
788         // check distance to bad channel
789         if (fUseDistanceToBadChannel){
790                 if (cluster->GetDistanceToBadChannel() < fMinDistanceToBadChannel){
791                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
792                         return kFALSE;
793                 }       
794         }
795         cutIndex++;
796         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
797
798         // Histos after cuts
799         if(fHistClusterEtavsPhiAfterAcc) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster);
800         
801         return kTRUE;
802 }
803
804 ///________________________________________________________________________
805 Bool_t AliCaloPhotonCuts::UpdateCutString() {
806    ///Update the cut string (if it has been created yet)
807
808    if(fCutString && fCutString->GetString().Length() == kNCuts) {
809       fCutString->SetString(GetCutNumber());
810    } else {
811       return kFALSE;
812    }
813    return kTRUE;
814 }
815
816 ///________________________________________________________________________
817 Bool_t AliCaloPhotonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
818         // Initialize Cuts from a given Cut string
819         AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
820         if(analysisCutSelection.Length()!=kNCuts) {
821                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
822                 return kFALSE;
823         }
824         if(!analysisCutSelection.IsDigit()){
825                 AliError("Cut selection contains characters");
826                 return kFALSE;
827         }
828
829         const char *cutSelection = analysisCutSelection.Data();
830         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
831         for(Int_t ii=0;ii<kNCuts;ii++){
832                 ASSIGNARRAY(ii);
833         }
834
835         // Set Individual Cuts
836         for(Int_t ii=0;ii<kNCuts;ii++){
837                 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
838         }
839         PrintCutsWithValues();
840         return kTRUE;
841 }
842
843 ///________________________________________________________________________
844 Bool_t AliCaloPhotonCuts::SetCut(cutIds cutID, const Int_t value) {
845         ///Set individual cut ID
846
847         switch (cutID) {                
848                 
849                 case kClusterType:
850                         if( SetClusterTypeCut(value)) {
851                                 fCuts[kClusterType] = value;
852                                 UpdateCutString();
853                                 return kTRUE;
854                         } else return kFALSE;
855                 
856                 case kEtaMin:
857                         if( SetMinEtaCut(value)) {
858                                 fCuts[kEtaMin] = value;
859                                 UpdateCutString();
860                                 return kTRUE;
861                         } else return kFALSE;
862
863                 case kEtaMax:
864                         if( SetMaxEtaCut(value)) {
865                                 fCuts[kEtaMax] = value;
866                                 UpdateCutString();
867                                 return kTRUE;
868                         } else return kFALSE;
869
870                 case kPhiMin:
871                         if( SetMinPhiCut(value)) {
872                                 fCuts[kPhiMin] = value;
873                                 UpdateCutString();
874                                 return kTRUE;
875                         } else return kFALSE;
876
877                 case kPhiMax:
878                         if( SetMaxPhiCut(value)) {
879                                 fCuts[kPhiMax] = value;
880                                 UpdateCutString();
881                                 return kTRUE;
882                         } else return kFALSE;
883
884                 case kDistanceToBadChannel:
885                         if( SetDistanceToBadChannelCut(value)) {
886                                 fCuts[kDistanceToBadChannel] = value;
887                                 UpdateCutString();
888                                 return kTRUE;
889                         } else return kFALSE;
890
891                 case kTiming:
892                         if( SetTimingCut(value)) {
893                                 fCuts[kTiming] = value;
894                                 UpdateCutString();
895                                 return kTRUE;
896                         } else return kFALSE;
897
898                 case kTrackMatching:
899                         if( SetTrackMatchingCut(value)) {
900                                 fCuts[kTrackMatching] = value;
901                                 UpdateCutString();
902                                 return kTRUE;
903                         } else return kFALSE;
904
905                 case kExoticCell:
906                         if( SetExoticCellCut(value)) {
907                                 fCuts[kExoticCell] = value;
908                                 UpdateCutString();
909                                 return kTRUE;
910                         } else return kFALSE;
911
912                 case kMinEnery:
913                         if( SetMinEnergyCut(value)) {
914                                 fCuts[kMinEnery] = value;
915                                 UpdateCutString();
916                                 return kTRUE;
917                         } else return kFALSE;
918
919                 case kNMinCells:
920                         if( SetMinNCellsCut(value)) {
921                                 fCuts[kNMinCells] = value;
922                                 UpdateCutString();
923                                 return kTRUE;
924                         } else return kFALSE;
925                         
926                 case kMinM02:
927                         if( SetMinM02(value)) {
928                                 fCuts[kMinM02] = value;
929                                 UpdateCutString();
930                                 return kTRUE;
931                         } else return kFALSE;
932
933                 case kMaxM02:
934                         if( SetMaxM02(value)) {
935                                 fCuts[kMaxM02] = value;
936                                 UpdateCutString();
937                                 return kTRUE;
938                         } else return kFALSE;
939                 
940                 case kMinM20:
941                         if( SetMinM20(value)) {
942                                 fCuts[kMinM20] = value;
943                                 UpdateCutString();
944                                 return kTRUE;
945                         } else return kFALSE;
946
947                 case kMaxM20:
948                         if( SetMaxM20(value)) {
949                                 fCuts[kMaxM20] = value;
950                                 UpdateCutString();
951                                 return kTRUE;
952                         } else return kFALSE;
953
954                 case kDispersion:
955                         if( SetDispersion(value)) {
956                                 fCuts[kDispersion] = value;
957                                 UpdateCutString();
958                                 return kTRUE;
959                         } else return kFALSE;
960
961                 case kNLM:
962                         if( SetNLM(value)) {
963                                 fCuts[kNLM] = value;
964                                 UpdateCutString();
965                                 return kTRUE;
966                         } else return kFALSE;
967
968                 case kNCuts:
969                         AliError("Cut id out of range");
970                         return kFALSE;
971         }
972
973         AliError("Cut id %d not recognized");
974         return kFALSE;
975
976
977 }
978 ///________________________________________________________________________
979 void AliCaloPhotonCuts::PrintCuts() {
980    // Print out current Cut Selection
981    for(Int_t ic = 0; ic < kNCuts; ic++) {
982       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
983    }
984 }
985
986 void AliCaloPhotonCuts::PrintCutsWithValues() {
987         // Print out current Cut Selection with value
988         printf("\nCluster cutnumber \n");
989         for(Int_t ic = 0; ic < kNCuts; ic++) {
990                 printf("%d",fCuts[ic]);
991         }
992         printf("\n\n");
993
994         printf("Acceptance cuts: \n");
995         if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
996         if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
997         if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
998         if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
999         if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
1000         if (fUseDistanceToBadChannel) printf("\tcut on exotics applied \n");
1001         
1002         printf("Cluster Quality cuts: \n");
1003         if (fUseTimeDiff) printf("\t time difference < %3.2f\n", fMaxTimeDiff );
1004         if (fUseDistTrackToCluster) printf("\tmin distance to track > %3.2f\n", fMinDistTrackToCluster );
1005         if (fUseExoticCell)printf("\t min distance to track > %3.2f\n", fMinDistTrackToCluster );
1006     if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
1007         if (fUseNCells) printf("\t number of cells per cluster > %d\n", fMinNCells );
1008         if (fUseM02) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
1009         if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
1010         if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
1011         if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
1012         
1013 }
1014
1015 // EMCAL acceptance 2011
1016 // 1.39626, 3.125 (phi)
1017 // -0.66687,,0.66465
1018
1019
1020 ///________________________________________________________________________
1021 Bool_t AliCaloPhotonCuts::SetClusterTypeCut(Int_t clusterType)
1022 {   // Set Cut
1023         switch(clusterType){
1024         case 0: // all clusters
1025                 fClusterType=0;
1026                 break;
1027         case 1: // EMCAL clusters
1028                 fClusterType=1;
1029                 break;
1030         case 2: // PHOS clusters
1031                 fClusterType=2;
1032                 break;
1033         default:
1034                 AliError(Form("ClusterTypeCut not defined %d",clusterType));
1035                 return kFALSE;
1036         }
1037         return kTRUE;
1038 }
1039
1040 //___________________________________________________________________
1041 Bool_t AliCaloPhotonCuts::SetMinEtaCut(Int_t minEta)
1042 {
1043         switch(minEta){
1044         case 0:
1045                 if (!fUseEtaCut) fUseEtaCut=0;
1046                 fMinEtaCut=-10.;
1047                 break;
1048         case 1:
1049                 if (!fUseEtaCut) fUseEtaCut=1;
1050                 fMinEtaCut=-0.6687;
1051                 break;
1052         case 2: 
1053                 if (!fUseEtaCut) fUseEtaCut=1;
1054                 fMinEtaCut=-0.5;
1055                 break;
1056         case 3: 
1057                 if (!fUseEtaCut) fUseEtaCut=1;
1058                 fMinEtaCut=-2;
1059                 break;
1060         default:
1061                 AliError(Form("MinEta Cut not defined %d",minEta));
1062                 return kFALSE;
1063         }
1064         return kTRUE;
1065 }
1066
1067
1068 //___________________________________________________________________
1069 Bool_t AliCaloPhotonCuts::SetMaxEtaCut(Int_t maxEta)
1070 {
1071         switch(maxEta){
1072         case 0: 
1073                 if (!fUseEtaCut) fUseEtaCut=0;
1074                 fMaxEtaCut=10;
1075                 break;          
1076         case 1:
1077                 if (!fUseEtaCut) fUseEtaCut=1;
1078                 fMaxEtaCut=0.66465;
1079                 break;
1080         case 2: 
1081                 if (!fUseEtaCut) fUseEtaCut=1;
1082                 fMaxEtaCut=0.5;
1083                 break;
1084         case 3: 
1085                 if (!fUseEtaCut) fUseEtaCut=1;
1086                 fMaxEtaCut=2;
1087                 break;
1088         default:
1089                 AliError(Form("MaxEta Cut not defined %d",maxEta));
1090                 return kFALSE;
1091         }
1092         return kTRUE;
1093 }
1094
1095 //___________________________________________________________________
1096 Bool_t AliCaloPhotonCuts::SetMinPhiCut(Int_t minPhi)
1097 {
1098         switch(minPhi){
1099         case 0: 
1100                 if (!fUsePhiCut) fUsePhiCut=0;
1101                 fMinPhiCut=-10000;
1102                 break;
1103         case 1: 
1104                 if (!fUsePhiCut) fUsePhiCut=1;
1105                 fMinPhiCut=1.39626;
1106                 break;
1107         default:
1108                 AliError(Form("MinPhi Cut not defined %d",minPhi));
1109                 return kFALSE;
1110         }
1111         return kTRUE;
1112 }
1113
1114 //___________________________________________________________________
1115 Bool_t AliCaloPhotonCuts::SetMaxPhiCut(Int_t maxPhi)
1116 {
1117         switch(maxPhi){
1118         case 0: 
1119                 if (!fUsePhiCut) fUsePhiCut=0;
1120                 fMaxPhiCut=10000;
1121                 break;
1122         case 1: 
1123                 if (!fUsePhiCut) fUsePhiCut=1;
1124                 fMaxPhiCut=3.125;
1125                 break;
1126         default:
1127                 AliError(Form("Max Phi Cut not defined %d",maxPhi));
1128                 return kFALSE;
1129         }
1130         return kTRUE;
1131 }
1132
1133 //___________________________________________________________________
1134 Bool_t AliCaloPhotonCuts::SetDistanceToBadChannelCut(Int_t distanceToBadChannel)
1135 {
1136         switch(distanceToBadChannel){
1137         case 0: 
1138                 fUseDistanceToBadChannel=0;
1139                 fMinDistanceToBadChannel=0;
1140                 break;
1141         case 1: 
1142                 if (!fUseDistanceToBadChannel) fUseDistanceToBadChannel=1;
1143                 fMinDistanceToBadChannel=5;
1144                 break;
1145         default:
1146                 AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
1147                 return kFALSE;
1148         }
1149         return kTRUE;
1150 }
1151
1152 //___________________________________________________________________
1153 Bool_t AliCaloPhotonCuts::SetTimingCut(Int_t timing)
1154 {
1155         switch(timing){
1156         case 0: 
1157                 fUseTimeDiff=0;
1158                 fMaxTimeDiff=500;
1159                 break;
1160         case 1: 
1161                 if (!fUseTimeDiff) fUseTimeDiff=1;
1162                 fMaxTimeDiff=10e-7; //1000ns
1163                 break;
1164         case 2: 
1165                 if (!fUseTimeDiff) fUseTimeDiff=1;
1166                 fMaxTimeDiff=50e-8; //500ns
1167                 break;
1168         case 3: 
1169                 if (!fUseTimeDiff) fUseTimeDiff=1;
1170                 fMaxTimeDiff=20e-8; //200ns
1171                 break;
1172         case 4: 
1173                 if (!fUseTimeDiff) fUseTimeDiff=1;
1174                 fMaxTimeDiff=10e-8; //100ns
1175                 break;
1176         case 5: 
1177                 if (!fUseTimeDiff) fUseTimeDiff=1;
1178                 fMaxTimeDiff=50e-9; //50ns
1179                 break;
1180
1181         default:
1182                 AliError(Form("Timing Cut not defined %d",timing));
1183                 return kFALSE;
1184         }
1185         return kTRUE;
1186 }
1187
1188 //___________________________________________________________________
1189 Bool_t AliCaloPhotonCuts::SetTrackMatchingCut(Int_t trackMatching)
1190 {
1191         switch(trackMatching){
1192         case 0: 
1193                 fUseDistTrackToCluster=0;
1194                 fMinDistTrackToCluster=0;
1195                 break;
1196         case 1: 
1197                 if (!fUseDistTrackToCluster) fUseDistTrackToCluster=1;
1198                 fMinDistTrackToCluster=5; 
1199                 break;
1200         default:
1201                 AliError(Form("Track Matching Cut not defined %d",trackMatching));
1202                 return kFALSE;
1203         }
1204         return kTRUE;
1205 }
1206
1207 //___________________________________________________________________
1208 Bool_t AliCaloPhotonCuts::SetExoticCellCut(Int_t exoticCell)
1209 {
1210         switch(exoticCell){
1211         case 0: 
1212                 fUseExoticCell=0;
1213                 fExoticCell=0;
1214                 break;
1215         case 1: 
1216                 if (!fUseExoticCell) fUseExoticCell=1;
1217                 fExoticCell=5; 
1218                 break;
1219         default:
1220                 AliError(Form("Exotic cell Cut not defined %d",exoticCell));
1221                 return kFALSE;
1222         }
1223         return kTRUE;
1224 }
1225                 
1226 //___________________________________________________________________
1227 Bool_t AliCaloPhotonCuts::SetMinEnergyCut(Int_t minEnergy)
1228 {
1229         switch(minEnergy){
1230         case 0: 
1231                 if (!fUseMinEnergy) fUseMinEnergy=0;
1232                 fMinEnergy=0;
1233                 break;
1234         case 1: 
1235                 if (!fUseMinEnergy) fUseMinEnergy=1;
1236                 fMinEnergy=0.05; 
1237                 break;
1238         case 2: 
1239                 if (!fUseMinEnergy) fUseMinEnergy=1;
1240                 fMinEnergy=0.1; 
1241                 break;
1242         case 3: 
1243                 if (!fUseMinEnergy) fUseMinEnergy=1;
1244                 fMinEnergy=0.15; 
1245                 break;
1246         case 4: 
1247                 if (!fUseMinEnergy) fUseMinEnergy=1;
1248                 fMinEnergy=0.2; 
1249                 break;
1250         case 5: 
1251                 if (!fUseMinEnergy) fUseMinEnergy=1;
1252                 fMinEnergy=0.3; 
1253                 break;
1254         case 6: 
1255                 if (!fUseMinEnergy) fUseMinEnergy=1;
1256                 fMinEnergy=0.5; 
1257                 break;
1258         case 7: 
1259                 if (!fUseMinEnergy) fUseMinEnergy=1;
1260                 fMinEnergy=0.75; 
1261                 break;
1262         case 8: 
1263                 if (!fUseMinEnergy) fUseMinEnergy=1;
1264                 fMinEnergy=1.; 
1265                 break;
1266         case 9: 
1267                 if (!fUseMinEnergy) fUseMinEnergy=1;
1268                 fMinEnergy=1.25; 
1269                 break;
1270         default:
1271                 AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
1272                 return kFALSE;
1273         }
1274         return kTRUE;
1275 }
1276                 
1277 //___________________________________________________________________
1278 Bool_t AliCaloPhotonCuts::SetMinNCellsCut(Int_t minNCells)
1279 {
1280         switch(minNCells){
1281         case 0:
1282                 if (!fUseNCells) fUseNCells=0;
1283                 fMinNCells=0;
1284                 break;
1285         case 1: 
1286                 if (!fUseNCells) fUseNCells=1;
1287                 fMinNCells=1; 
1288                 break;
1289         case 2: 
1290                 if (!fUseNCells) fUseNCells=1;
1291                 fMinNCells=2; 
1292                 break;
1293         case 3: 
1294                 if (!fUseNCells) fUseNCells=1;
1295                 fMinNCells=3; 
1296                 break;
1297         case 4: 
1298                 if (!fUseNCells) fUseNCells=1;
1299                 fMinNCells=4; 
1300                 break;
1301         case 5: 
1302                 if (!fUseNCells) fUseNCells=1;
1303                 fMinNCells=5; 
1304                 break;
1305         case 6: 
1306                 if (!fUseNCells) fUseNCells=1;
1307                 fMinNCells=6; 
1308                 break;
1309
1310         default:
1311                 AliError(Form("Min N cells Cut not defined %d",minNCells));
1312                 return kFALSE;
1313         }
1314         return kTRUE;
1315 }
1316
1317 //___________________________________________________________________
1318 Bool_t AliCaloPhotonCuts::SetMaxM02(Int_t maxM02)
1319 {
1320         switch(maxM02){
1321         case 0: 
1322                 if (!fUseM02) fUseM02=0;
1323                 fMaxM02=100;
1324                 break;
1325         case 1: 
1326                 if (!fUseM02) fUseM02=1;
1327                 fMaxM02=1.; 
1328                 break;
1329         case 2: 
1330                 if (!fUseM02) fUseM02=1;
1331                 fMaxM02=0.7; 
1332                 break;
1333         case 3: 
1334                 if (!fUseM02) fUseM02=1;
1335                 fMaxM02=0.5; 
1336                 break;
1337         case 4: 
1338                 if (!fUseM02) fUseM02=1;
1339                 fMaxM02=0.4; 
1340                 break;
1341         default:
1342                 AliError(Form("Max M02 Cut not defined %d",maxM02));
1343                 return kFALSE;
1344         }
1345         return kTRUE;
1346 }
1347
1348 //___________________________________________________________________
1349 Bool_t AliCaloPhotonCuts::SetMinM02(Int_t minM02)
1350 {
1351         switch(minM02){
1352         case 0: 
1353                 if (!fUseM02) fUseM02=0;
1354                 fMinM02=0;
1355                 break;
1356         case 1: 
1357                 if (!fUseM02) fUseM02=1;
1358                 fMinM02=0.002; 
1359                 break;
1360         default:
1361                 AliError(Form("Min M02 not defined %d",minM02));
1362                 return kFALSE;
1363         }
1364         return kTRUE;
1365 }
1366
1367 //___________________________________________________________________
1368 Bool_t AliCaloPhotonCuts::SetMaxM20(Int_t maxM20)
1369 {
1370         switch(maxM20){
1371         case 0: 
1372                 if (!fUseM20) fUseM20=0;
1373                 fMaxM20=100;
1374                 break;
1375         case 1: 
1376                 if (!fUseM20) fUseM20=1;
1377                 fMaxM20=0.5; 
1378                 break;
1379         default:
1380                 AliError(Form("Max M20 Cut not defined %d",maxM20));
1381                 return kFALSE;
1382         }
1383         return kTRUE;
1384 }
1385
1386 //___________________________________________________________________
1387 Bool_t AliCaloPhotonCuts::SetMinM20(Int_t minM20)
1388 {
1389         switch(minM20){
1390         case 0: 
1391                 if (!fUseM20) fUseM20=0;
1392                 fMinM20=0;
1393                 break;
1394         case 1: 
1395                 if (!fUseM20) fUseM20=1;
1396                 fMinM20=0.002; 
1397                 break;
1398         default:
1399                 AliError(Form("Min M20 Cut not defined %d",minM20));
1400                 return kFALSE;
1401         }
1402         return kTRUE;
1403 }
1404
1405 //___________________________________________________________________
1406 Bool_t AliCaloPhotonCuts::SetDispersion(Int_t dispersion)
1407 {
1408         switch(dispersion){
1409         case 0: 
1410                 if (!fUseDispersion) fUseDispersion=0;
1411                 fMaxDispersion =100;
1412                 break;
1413         case 1: 
1414                 if (!fUseDispersion) fUseDispersion=1;
1415                 fMaxDispersion=2.; 
1416                 break;
1417         default:
1418                 AliError(Form("Maximum Dispersion Cut not defined %d",dispersion));
1419                 return kFALSE;
1420         }
1421         return kTRUE;
1422 }
1423
1424 //___________________________________________________________________
1425 Bool_t AliCaloPhotonCuts::SetNLM(Int_t nlm)
1426 {
1427         switch(nlm){
1428         case 0: 
1429                 if (!fUseNLM) fUseNLM=0;
1430                 fMinNLM =0;
1431                 fMaxNLM =100;
1432                 break;
1433         case 1: 
1434                 if (!fUseNLM) fUseNLM=1;
1435                 fMinNLM =0;
1436                 fMaxNLM =1;
1437                 break;
1438         default:
1439                 AliError(Form("NLM Cut not defined %d",nlm));
1440                 return kFALSE;
1441         }
1442         return kTRUE;
1443 }
1444         
1445 ///________________________________________________________________________
1446 TString AliCaloPhotonCuts::GetCutNumber(){
1447    // returns TString with current cut number
1448    TString a(kNCuts);
1449    for(Int_t ii=0;ii<kNCuts;ii++){
1450       a.Append(Form("%d",fCuts[ii]));
1451    }
1452    return a;
1453 }
1454         
1455