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 //
 579 //              //loop over tracks for Jet QA
 580 //              TList *l = event->GetList();
 581 //              TClonesArray *tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
 582 //
 583 //              for(int itrack = 0; itrack < NtrMatched; itrack++){
 584 //                      AliVTrack *trackcluster = static_cast<AliVTrack*>(tracks->At(itrack));
 585 //                      if (! trackcluster) {
 586 //                              AliError(Form("Couldn't get ESD track %d\n", itrack));
 587 //                              continue;
 588 //                      }
 589 //                      Double_t dphi = -999.0;
 590 //                      Double_t deta = -999.0;
 591 //                      AliPicoTrack::GetEtaPhiDiff(trackcluster, cluster, dphi, deta);
 592 //                      cout << "here 2" << endl;
 593 //                      Double_t dr = sqrt(dphi*dphi + deta+deta);
 594 //                      if(dr < fMinDistTrackToCluster){
 595 //                              if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2
 596 //                              return kFALSE;
 597 //                      }
 598 //
 599 //              }//loop over tracks
 600
 601
 602                 if(cluster->GetEmcCpvDistance() < fMinDistTrackToCluster){
 603                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //2
 604                         return kFALSE;
 605                 }
 606         }
 607         cutIndex++;//3, next cut
 608
 609         // exotic cell cut --IMPLEMENT LATER---
 610 //      if(!AcceptanceCuts(photon)){
 611 //              if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //3
 612 //              return kFALSE;
 613 //      }
 614         cutIndex++; //4, next cut
 615
 616         // minimum cell energy cut
 617         if (fUseMinEnergy){
 618                 if(cluster->E() < fMinEnergy){
 619                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //4
 620                         return kFALSE;
 621                 }
 622         }
 623         cutIndex++; //5, next cut
 624
 625         // minimum number of cells
 626         if (fUseNCells){
 627                 if(cluster->GetNCells() < fMinNCells) {
 628                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //5
 629                         return kFALSE;
 630                 }
 631         }
 632         cutIndex++; //6, next cut
 633
 634         // M02 cut
 635         if (fUseM02){
 636                 if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 ) {
 637                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //6
 638                         return kFALSE;
 639                 }
 640         }
 641         cutIndex++; //7, next cut
 642
 643         // M20 cut
 644         if (fUseM20){
 645                 if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 ) {
 646                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //7
 647                         return kFALSE;
 648                 }
 649         }
 650         cutIndex++; //8, next cut
 651
 652         // dispersion cut
 653         if (fUseDispersion){
 654                 if( cluster->GetDispersion()> fMaxDispersion) {
 655                         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //8
 656                         return kFALSE;
 657                 }
 658         }
 659         cutIndex++; //9, next cut
 660
 661         // NLM cut --IMPLEMENT LATER---
 662 //      if (fUseNLM){
 663 //              if( cluster->GetDispersion()> fMaxDispersion) {
 664 //                      if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //9
 665 //                      return kFALSE;
 666 //              }
 667 //      }
 668         cutIndex++; //9, next cut
 669
 670         // DONE with selecting photons
 671         if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex); //10
 672
 673         // Histos after Cuts
 674         Double_t vertex[3] = {0};
 675         event->GetPrimaryVertex()->GetXYZ(vertex);
 676         // TLorentzvector with cluster
 677         TLorentzVector clusterVector;
 678         cluster->GetMomentum(clusterVector,vertex);
 679         Double_t etaCluster = clusterVector.Eta();
 680         Double_t phiCluster = clusterVector.Phi();
 681
 682         if(fHistClusterEtavsPhiAfterQA) fHistClusterEtavsPhiAfterQA->Fill(phiCluster,etaCluster);
 683         if(fHistClusterTimevsEAfterQA) fHistClusterTimevsEAfterQA->Fill(cluster->GetTOF(), cluster->E());
 684 //      if(fHistExoticCellAfterQA) fHistExoticCellAfterQA->Fill(cluster->E(), );
 685         if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(minR);
 686         if(fHistDistanceTrackToClusterAfterQA) fHistDistanceTrackToClusterAfterQA->Fill(cluster->GetEmcCpvDistance());
 687         if(fHistEnergyOfClusterAfterQA) fHistEnergyOfClusterAfterQA->Fill(cluster->E());
 688         if(fHistNCellsAfterQA) fHistNCellsAfterQA->Fill(cluster->GetNCells());
 689         if(fHistM02AfterQA) fHistM02AfterQA->Fill(cluster->GetM02());
 690         if(fHistM20AfterQA) fHistM20AfterQA->Fill(cluster->GetM20());
 691         if(fHistDispersionAfterQA) fHistDispersionAfterQA->Fill(cluster->GetDispersion());
 692 //      if(fHistNLMBeforeQA) fHistNLMAfterQA->Fill(cluster->GetNExMax());
 693
 694         return kTRUE;
 695
 696 }
 697
 698
 699 ///________________________________________________________________________
 700 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, Bool_t isMC)
 701 {
 702         //Selection of Reconstructed photon clusters with Calorimeters
 703
 704         FillClusterCutIndex(kPhotonIn);
 705
 706         Double_t vertex[3] = {0};
 707         event->GetPrimaryVertex()->GetXYZ(vertex);
 708         // TLorentzvector with cluster
 709         TLorentzVector clusterVector;
 710         cluster->GetMomentum(clusterVector,vertex);
 711         Double_t etaCluster = clusterVector.Eta();
 712         Double_t phiCluster = clusterVector.Phi();
 713
 714         // Histos before cuts
 715         if(fHistClusterEtavsPhiBeforeAcc) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster);
 716
 717         // Cluster Selection - 0= accept any calo cluster
 718         if (fClusterType > 0){
 719                 //Select EMCAL cluster
 720                 if (fClusterType == 1 && !cluster->IsEMCAL()){
 721                         FillClusterCutIndex(kDetector);
 722                         return kFALSE;
 723                 }
 724                 //Select PHOS cluster
 725                 if (fClusterType == 2 && !cluster->IsPHOS()){
 726                         FillClusterCutIndex(kDetector);
 727                         return kFALSE;
 728                 }
 729         }
 730
 731         // Acceptance Cuts
 732         if(!AcceptanceCuts(cluster,event)){
 733                 FillClusterCutIndex(kAcceptance);
 734                 return kFALSE;
 735         }
 736         // Cluster Quality Cuts
 737         if(!ClusterQualityCuts(cluster,event,isMC)){
 738                 FillClusterCutIndex(kClusterQuality);
 739                 return kFALSE;
 740         }
 741
 742         // Photon passed cuts
 743         FillClusterCutIndex(kPhotonOut);
 744         return kTRUE;
 745 }
 746
 747
 748 ///________________________________________________________________________
 749 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event)
 750 {
 751    // Exclude certain areas for photon reconstruction
 752
 753         Int_t cutIndex=0;
 754         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
 755         cutIndex++;
 756
 757
 758         Double_t vertex[3] = {0};
 759         event->GetPrimaryVertex()->GetXYZ(vertex);
 760         // TLorentzvector with cluster
 761         TLorentzVector clusterVector;
 762         cluster->GetMomentum(clusterVector,vertex);
 763         Double_t etaCluster = clusterVector.Eta();
 764         Double_t phiCluster = clusterVector.Phi();
 765
 766         // check eta range
 767         if (fUseEtaCut){
 768                 if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut){
 769                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
 770                         return kFALSE;
 771                 }
 772         }
 773         cutIndex++;
 774
 775         // check phi range
 776         if (fUsePhiCut ){
 777                 if (phiCluster < fMinPhiCut || phiCluster > fMaxEtaCut){
 778                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
 779                         return kFALSE;
 780                 }
 781         }
 782         cutIndex++;
 783
 784         // check distance to bad channel
 785         if (fUseDistanceToBadChannel){
 786                 if (cluster->GetDistanceToBadChannel() < fMinDistanceToBadChannel){
 787                         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
 788                         return kFALSE;
 789                 }
 790         }
 791         cutIndex++;
 792         if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
 793
 794         // Histos after cuts
 795         if(fHistClusterEtavsPhiAfterAcc) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster);
 796
 797         return kTRUE;
 798 }
 799
 800 ///________________________________________________________________________
 801 Bool_t AliCaloPhotonCuts::UpdateCutString() {
 802    ///Update the cut string (if it has been created yet)
 803
 804    if(fCutString && fCutString->GetString().Length() == kNCuts) {
 805       fCutString->SetString(GetCutNumber());
 806    } else {
 807       return kFALSE;
 808    }
 809    return kTRUE;
 810 }
 811
 812 ///________________________________________________________________________
 813 Bool_t AliCaloPhotonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 814         // Initialize Cuts from a given Cut string
 815         AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
 816         if(analysisCutSelection.Length()!=kNCuts) {
 817                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 818                 return kFALSE;
 819         }
 820         if(!analysisCutSelection.IsDigit()){
 821                 AliError("Cut selection contains characters");
 822                 return kFALSE;
 823         }
 824
 825         const char *cutSelection = analysisCutSelection.Data();
 826         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 827         for(Int_t ii=0;ii<kNCuts;ii++){
 828                 ASSIGNARRAY(ii);
 829         }
 830
 831         // Set Individual Cuts
 832         for(Int_t ii=0;ii<kNCuts;ii++){
 833                 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 834         }
 835         PrintCutsWithValues();
 836         return kTRUE;
 837 }
 838
 839 ///________________________________________________________________________
 840 Bool_t AliCaloPhotonCuts::SetCut(cutIds cutID, const Int_t value) {
 841         ///Set individual cut ID
 842
 843         switch (cutID) {
 844
 845                 case kClusterType:
 846                         if( SetClusterTypeCut(value)) {
 847                                 fCuts[kClusterType] = value;
 848                                 UpdateCutString();
 849                                 return kTRUE;
 850                         } else return kFALSE;
 851
 852                 case kEtaMin:
 853                         if( SetMinEtaCut(value)) {
 854                                 fCuts[kEtaMin] = value;
 855                                 UpdateCutString();
 856                                 return kTRUE;
 857                         } else return kFALSE;
 858
 859                 case kEtaMax:
 860                         if( SetMaxEtaCut(value)) {
 861                                 fCuts[kEtaMax] = value;
 862                                 UpdateCutString();
 863                                 return kTRUE;
 864                         } else return kFALSE;
 865
 866                 case kPhiMin:
 867                         if( SetMinPhiCut(value)) {
 868                                 fCuts[kPhiMin] = value;
 869                                 UpdateCutString();
 870                                 return kTRUE;
 871                         } else return kFALSE;
 872
 873                 case kPhiMax:
 874                         if( SetMaxPhiCut(value)) {
 875                                 fCuts[kPhiMax] = value;
 876                                 UpdateCutString();
 877                                 return kTRUE;
 878                         } else return kFALSE;
 879
 880                 case kDistanceToBadChannel:
 881                         if( SetDistanceToBadChannelCut(value)) {
 882                                 fCuts[kDistanceToBadChannel] = value;
 883                                 UpdateCutString();
 884                                 return kTRUE;
 885                         } else return kFALSE;
 886
 887                 case kTiming:
 888                         if( SetTimingCut(value)) {
 889                                 fCuts[kTiming] = value;
 890                                 UpdateCutString();
 891                                 return kTRUE;
 892                         } else return kFALSE;
 893
 894                 case kTrackMatching:
 895                         if( SetTrackMatchingCut(value)) {
 896                                 fCuts[kTrackMatching] = value;
 897                                 UpdateCutString();
 898                                 return kTRUE;
 899                         } else return kFALSE;
 900
 901                 case kExoticCell:
 902                         if( SetExoticCellCut(value)) {
 903                                 fCuts[kExoticCell] = value;
 904                                 UpdateCutString();
 905                                 return kTRUE;
 906                         } else return kFALSE;
 907
 908                 case kMinEnery:
 909                         if( SetMinEnergyCut(value)) {
 910                                 fCuts[kMinEnery] = value;
 911                                 UpdateCutString();
 912                                 return kTRUE;
 913                         } else return kFALSE;
 914
 915                 case kNMinCells:
 916                         if( SetMinNCellsCut(value)) {
 917                                 fCuts[kNMinCells] = value;
 918                                 UpdateCutString();
 919                                 return kTRUE;
 920                         } else return kFALSE;
 921
 922                 case kMinM02:
 923                         if( SetMinM02(value)) {
 924                                 fCuts[kMinM02] = value;
 925                                 UpdateCutString();
 926                                 return kTRUE;
 927                         } else return kFALSE;
 928
 929                 case kMaxM02:
 930                         if( SetMaxM02(value)) {
 931                                 fCuts[kMaxM02] = value;
 932                                 UpdateCutString();
 933                                 return kTRUE;
 934                         } else return kFALSE;
 935
 936                 case kMinM20:
 937                         if( SetMinM20(value)) {
 938                                 fCuts[kMinM20] = value;
 939                                 UpdateCutString();
 940                                 return kTRUE;
 941                         } else return kFALSE;
 942
 943                 case kMaxM20:
 944                         if( SetMaxM20(value)) {
 945                                 fCuts[kMaxM20] = value;
 946                                 UpdateCutString();
 947                                 return kTRUE;
 948                         } else return kFALSE;
 949
 950                 case kDispersion:
 951                         if( SetDispersion(value)) {
 952                                 fCuts[kDispersion] = value;
 953                                 UpdateCutString();
 954                                 return kTRUE;
 955                         } else return kFALSE;
 956
 957                 case kNLM:
 958                         if( SetNLM(value)) {
 959                                 fCuts[kNLM] = value;
 960                                 UpdateCutString();
 961                                 return kTRUE;
 962                         } else return kFALSE;
 963
 964                 case kNCuts:
 965                         AliError("Cut id out of range");
 966                         return kFALSE;
 967         }
 968
 969         AliError("Cut id %d not recognized");
 970         return kFALSE;
 971
 972
 973 }
 974 ///________________________________________________________________________
 975 void AliCaloPhotonCuts::PrintCuts() {
 976    // Print out current Cut Selection
 977    for(Int_t ic = 0; ic < kNCuts; ic++) {
 978       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 979    }
 980 }
 981
 982 void AliCaloPhotonCuts::PrintCutsWithValues() {
 983         // Print out current Cut Selection with value
 984         printf("\nCluster cutnumber \n");
 985         for(Int_t ic = 0; ic < kNCuts; ic++) {
 986                 printf("%d",fCuts[ic]);
 987         }
 988         printf("\n\n");
 989
 990         printf("Acceptance cuts: \n");
 991         if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
 992         if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
 993         if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
 994         if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
 995         if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
 996         if (fUseDistanceToBadChannel) printf("\tcut on exotics applied \n");
 997
 998         printf("Cluster Quality cuts: \n");
 999         if (fUseTimeDiff) printf("\t time difference < %3.2f\n", fMaxTimeDiff );
1000         if (fUseDistTrackToCluster) printf("\tmin distance to track > %3.2f\n", fMinDistTrackToCluster );
1001         if (fUseExoticCell)printf("\t min distance to track > %3.2f\n", fMinDistTrackToCluster );
1002     if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
1003         if (fUseNCells) printf("\t number of cells per cluster > %d\n", fMinNCells );
1004         if (fUseM02) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
1005         if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
1006         if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
1007         if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
1008
1009 }
1010
1011 // EMCAL acceptance 2011
1012 // 1.39626, 3.125 (phi)
1013 // -0.66687,,0.66465
1014
1015
1016 ///________________________________________________________________________
1017 Bool_t AliCaloPhotonCuts::SetClusterTypeCut(Int_t clusterType)
1018 {   // Set Cut
1019         switch(clusterType){
1020         case 0: // all clusters
1021                 fClusterType=0;
1022                 break;
1023         case 1: // EMCAL clusters
1024                 fClusterType=1;
1025                 break;
1026         case 2: // PHOS clusters
1027                 fClusterType=2;
1028                 break;
1029         default:
1030                 AliError(Form("ClusterTypeCut not defined %d",clusterType));
1031                 return kFALSE;
1032         }
1033         return kTRUE;
1034 }
1035
1036 //___________________________________________________________________
1037 Bool_t AliCaloPhotonCuts::SetMinEtaCut(Int_t minEta)
1038 {
1039         switch(minEta){
1040         case 0:
1041                 if (!fUseEtaCut) fUseEtaCut=0;
1042                 fMinEtaCut=-10.;
1043                 break;
1044         case 1:
1045                 if (!fUseEtaCut) fUseEtaCut=1;
1046                 fMinEtaCut=-0.6687;
1047                 break;
1048         case 2:
1049                 if (!fUseEtaCut) fUseEtaCut=1;
1050                 fMinEtaCut=-0.5;
1051                 break;
1052         case 3:
1053                 if (!fUseEtaCut) fUseEtaCut=1;
1054                 fMinEtaCut=-2;
1055                 break;
1056         default:
1057                 AliError(Form("MinEta Cut not defined %d",minEta));
1058                 return kFALSE;
1059         }
1060         return kTRUE;
1061 }
1062
1063
1064 //___________________________________________________________________
1065 Bool_t AliCaloPhotonCuts::SetMaxEtaCut(Int_t maxEta)
1066 {
1067         switch(maxEta){
1068         case 0:
1069                 if (!fUseEtaCut) fUseEtaCut=0;
1070                 fMaxEtaCut=10;
1071                 break;
1072         case 1:
1073                 if (!fUseEtaCut) fUseEtaCut=1;
1074                 fMaxEtaCut=0.66465;
1075                 break;
1076         case 2:
1077                 if (!fUseEtaCut) fUseEtaCut=1;
1078                 fMaxEtaCut=0.5;
1079                 break;
1080         case 3:
1081                 if (!fUseEtaCut) fUseEtaCut=1;
1082                 fMaxEtaCut=2;
1083                 break;
1084         default:
1085                 AliError(Form("MaxEta Cut not defined %d",maxEta));
1086                 return kFALSE;
1087         }
1088         return kTRUE;
1089 }
1090
1091 //___________________________________________________________________
1092 Bool_t AliCaloPhotonCuts::SetMinPhiCut(Int_t minPhi)
1093 {
1094         switch(minPhi){
1095         case 0:
1096                 if (!fUsePhiCut) fUsePhiCut=0;
1097                 fMinPhiCut=-10000;
1098                 break;
1099         case 1:
1100                 if (!fUsePhiCut) fUsePhiCut=1;
1101                 fMinPhiCut=1.39626;
1102                 break;
1103         default:
1104                 AliError(Form("MinPhi Cut not defined %d",minPhi));
1105                 return kFALSE;
1106         }
1107         return kTRUE;
1108 }
1109
1110 //___________________________________________________________________
1111 Bool_t AliCaloPhotonCuts::SetMaxPhiCut(Int_t maxPhi)
1112 {
1113         switch(maxPhi){
1114         case 0:
1115                 if (!fUsePhiCut) fUsePhiCut=0;
1116                 fMaxPhiCut=10000;
1117                 break;
1118         case 1:
1119                 if (!fUsePhiCut) fUsePhiCut=1;
1120                 fMaxPhiCut=3.125;
1121                 break;
1122         default:
1123                 AliError(Form("Max Phi Cut not defined %d",maxPhi));
1124                 return kFALSE;
1125         }
1126         return kTRUE;
1127 }
1128
1129 //___________________________________________________________________
1130 Bool_t AliCaloPhotonCuts::SetDistanceToBadChannelCut(Int_t distanceToBadChannel)
1131 {
1132         switch(distanceToBadChannel){
1133         case 0:
1134                 fUseDistanceToBadChannel=0;
1135                 fMinDistanceToBadChannel=0;
1136                 break;
1137         case 1:
1138                 if (!fUseDistanceToBadChannel) fUseDistanceToBadChannel=1;
1139                 fMinDistanceToBadChannel=5;
1140                 break;
1141         default:
1142                 AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
1143                 return kFALSE;
1144         }
1145         return kTRUE;
1146 }
1147
1148 //___________________________________________________________________
1149 Bool_t AliCaloPhotonCuts::SetTimingCut(Int_t timing)
1150 {
1151         switch(timing){
1152         case 0:
1153                 fUseTimeDiff=0;
1154                 fMaxTimeDiff=500;
1155                 break;
1156         case 1:
1157                 if (!fUseTimeDiff) fUseTimeDiff=1;
1158                 fMaxTimeDiff=10e-7; //1000ns
1159                 break;
1160         case 2:
1161                 if (!fUseTimeDiff) fUseTimeDiff=1;
1162                 fMaxTimeDiff=50e-8; //500ns
1163                 break;
1164         case 3:
1165                 if (!fUseTimeDiff) fUseTimeDiff=1;
1166                 fMaxTimeDiff=20e-8; //200ns
1167                 break;
1168         case 4:
1169                 if (!fUseTimeDiff) fUseTimeDiff=1;
1170                 fMaxTimeDiff=10e-8; //100ns
1171                 break;
1172         case 5:
1173                 if (!fUseTimeDiff) fUseTimeDiff=1;
1174                 fMaxTimeDiff=50e-9; //50ns
1175                 break;
1176
1177         default:
1178                 AliError(Form("Timing Cut not defined %d",timing));
1179                 return kFALSE;
1180         }
1181         return kTRUE;
1182 }
1183
1184 //___________________________________________________________________
1185 Bool_t AliCaloPhotonCuts::SetTrackMatchingCut(Int_t trackMatching)
1186 {
1187         switch(trackMatching){
1188         case 0:
1189                 fUseDistTrackToCluster=0;
1190                 fMinDistTrackToCluster=0;
1191                 break;
1192         case 1:
1193                 if (!fUseDistTrackToCluster) fUseDistTrackToCluster=1;
1194                 fMinDistTrackToCluster=5;
1195                 break;
1196         default:
1197                 AliError(Form("Track Matching Cut not defined %d",trackMatching));
1198                 return kFALSE;
1199         }
1200         return kTRUE;
1201 }
1202
1203 //___________________________________________________________________
1204 Bool_t AliCaloPhotonCuts::SetExoticCellCut(Int_t exoticCell)
1205 {
1206         switch(exoticCell){
1207         case 0:
1208                 fUseExoticCell=0;
1209                 fExoticCell=0;
1210                 break;
1211         case 1:
1212                 if (!fUseExoticCell) fUseExoticCell=1;
1213                 fExoticCell=5;
1214                 break;
1215         default:
1216                 AliError(Form("Exotic cell Cut not defined %d",exoticCell));
1217                 return kFALSE;
1218         }
1219         return kTRUE;
1220 }
1221
1222 //___________________________________________________________________
1223 Bool_t AliCaloPhotonCuts::SetMinEnergyCut(Int_t minEnergy)
1224 {
1225         switch(minEnergy){
1226         case 0:
1227                 if (!fUseMinEnergy) fUseMinEnergy=0;
1228                 fMinEnergy=0;
1229                 break;
1230         case 1:
1231                 if (!fUseMinEnergy) fUseMinEnergy=1;
1232                 fMinEnergy=0.05;
1233                 break;
1234         case 2:
1235                 if (!fUseMinEnergy) fUseMinEnergy=1;
1236                 fMinEnergy=0.1;
1237                 break;
1238         case 3:
1239                 if (!fUseMinEnergy) fUseMinEnergy=1;
1240                 fMinEnergy=0.15;
1241                 break;
1242         case 4:
1243                 if (!fUseMinEnergy) fUseMinEnergy=1;
1244                 fMinEnergy=0.2;
1245                 break;
1246         case 5:
1247                 if (!fUseMinEnergy) fUseMinEnergy=1;
1248                 fMinEnergy=0.3;
1249                 break;
1250         case 6:
1251                 if (!fUseMinEnergy) fUseMinEnergy=1;
1252                 fMinEnergy=0.5;
1253                 break;
1254         case 7:
1255                 if (!fUseMinEnergy) fUseMinEnergy=1;
1256                 fMinEnergy=0.75;
1257                 break;
1258         case 8:
1259                 if (!fUseMinEnergy) fUseMinEnergy=1;
1260                 fMinEnergy=1.;
1261                 break;
1262         case 9:
1263                 if (!fUseMinEnergy) fUseMinEnergy=1;
1264                 fMinEnergy=1.25;
1265                 break;
1266         default:
1267                 AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
1268                 return kFALSE;
1269         }
1270         return kTRUE;
1271 }
1272
1273 //___________________________________________________________________
1274 Bool_t AliCaloPhotonCuts::SetMinNCellsCut(Int_t minNCells)
1275 {
1276         switch(minNCells){
1277         case 0:
1278                 if (!fUseNCells) fUseNCells=0;
1279                 fMinNCells=0;
1280                 break;
1281         case 1:
1282                 if (!fUseNCells) fUseNCells=1;
1283                 fMinNCells=1;
1284                 break;
1285         case 2:
1286                 if (!fUseNCells) fUseNCells=1;
1287                 fMinNCells=2;
1288                 break;
1289         case 3:
1290                 if (!fUseNCells) fUseNCells=1;
1291                 fMinNCells=3;
1292                 break;
1293         case 4:
1294                 if (!fUseNCells) fUseNCells=1;
1295                 fMinNCells=4;
1296                 break;
1297         case 5:
1298                 if (!fUseNCells) fUseNCells=1;
1299                 fMinNCells=5;
1300                 break;
1301         case 6:
1302                 if (!fUseNCells) fUseNCells=1;
1303                 fMinNCells=6;
1304                 break;
1305
1306         default:
1307                 AliError(Form("Min N cells Cut not defined %d",minNCells));
1308                 return kFALSE;
1309         }
1310         return kTRUE;
1311 }
1312
1313 //___________________________________________________________________
1314 Bool_t AliCaloPhotonCuts::SetMaxM02(Int_t maxM02)
1315 {
1316         switch(maxM02){
1317         case 0:
1318                 if (!fUseM02) fUseM02=0;
1319                 fMaxM02=100;
1320                 break;
1321         case 1:
1322                 if (!fUseM02) fUseM02=1;
1323                 fMaxM02=1.;
1324                 break;
1325         case 2:
1326                 if (!fUseM02) fUseM02=1;
1327                 fMaxM02=0.7;
1328                 break;
1329         case 3:
1330                 if (!fUseM02) fUseM02=1;
1331                 fMaxM02=0.5;
1332                 break;
1333         case 4:
1334                 if (!fUseM02) fUseM02=1;
1335                 fMaxM02=0.4;
1336                 break;
1337         default:
1338                 AliError(Form("Max M02 Cut not defined %d",maxM02));
1339                 return kFALSE;
1340         }
1341         return kTRUE;
1342 }
1343
1344 //___________________________________________________________________
1345 Bool_t AliCaloPhotonCuts::SetMinM02(Int_t minM02)
1346 {
1347         switch(minM02){
1348         case 0:
1349                 if (!fUseM02) fUseM02=0;
1350                 fMinM02=0;
1351                 break;
1352         case 1:
1353                 if (!fUseM02) fUseM02=1;
1354                 fMinM02=0.002;
1355                 break;
1356         default:
1357                 AliError(Form("Min M02 not defined %d",minM02));
1358                 return kFALSE;
1359         }
1360         return kTRUE;
1361 }
1362
1363 //___________________________________________________________________
1364 Bool_t AliCaloPhotonCuts::SetMaxM20(Int_t maxM20)
1365 {
1366         switch(maxM20){
1367         case 0:
1368                 if (!fUseM20) fUseM20=0;
1369                 fMaxM20=100;
1370                 break;
1371         case 1:
1372                 if (!fUseM20) fUseM20=1;
1373                 fMaxM20=0.5;
1374                 break;
1375         default:
1376                 AliError(Form("Max M20 Cut not defined %d",maxM20));
1377                 return kFALSE;
1378         }
1379         return kTRUE;
1380 }
1381
1382 //___________________________________________________________________
1383 Bool_t AliCaloPhotonCuts::SetMinM20(Int_t minM20)
1384 {
1385         switch(minM20){
1386         case 0:
1387                 if (!fUseM20) fUseM20=0;
1388                 fMinM20=0;
1389                 break;
1390         case 1:
1391                 if (!fUseM20) fUseM20=1;
1392                 fMinM20=0.002;
1393                 break;
1394         default:
1395                 AliError(Form("Min M20 Cut not defined %d",minM20));
1396                 return kFALSE;
1397         }
1398         return kTRUE;
1399 }
1400
1401 //___________________________________________________________________
1402 Bool_t AliCaloPhotonCuts::SetDispersion(Int_t dispersion)
1403 {
1404         switch(dispersion){
1405         case 0:
1406                 if (!fUseDispersion) fUseDispersion=0;
1407                 fMaxDispersion =100;
1408                 break;
1409         case 1:
1410                 if (!fUseDispersion) fUseDispersion=1;
1411                 fMaxDispersion=2.;
1412                 break;
1413         default:
1414                 AliError(Form("Maximum Dispersion Cut not defined %d",dispersion));
1415                 return kFALSE;
1416         }
1417         return kTRUE;
1418 }
1419
1420 //___________________________________________________________________
1421 Bool_t AliCaloPhotonCuts::SetNLM(Int_t nlm)
1422 {
1423         switch(nlm){
1424         case 0:
1425                 if (!fUseNLM) fUseNLM=0;
1426                 fMinNLM =0;
1427                 fMaxNLM =100;
1428                 break;
1429         case 1:
1430                 if (!fUseNLM) fUseNLM=1;
1431                 fMinNLM =0;
1432                 fMaxNLM =1;
1433                 break;
1434         default:
1435                 AliError(Form("NLM Cut not defined %d",nlm));
1436                 return kFALSE;
1437         }
1438         return kTRUE;
1439 }
1440
1441 ///________________________________________________________________________
1442 TString AliCaloPhotonCuts::GetCutNumber(){
1443    // returns TString with current cut number
1444    TString a(kNCuts);
1445    for(Int_t ii=0;ii<kNCuts;ii++){
1446       a.Append(Form("%d",fCuts[ii]));
1447    }
1448    return a;
1449 }
1450
1451