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