PWGGA/GammaConv/AliConversionMesonCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Authors: Svein Lindal, Daniel Lohner                                   *
   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 // Gamma Conversion analysis
  20 //---------------------------------------------
  21 ////////////////////////////////////////////////
  22
  23 #include "AliConversionMesonCuts.h"
  24
  25 #include "AliKFVertex.h"
  26 #include "AliAODTrack.h"
  27 #include "AliESDtrack.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliInputEventHandler.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliAODHandler.h"
  32 #include "AliPIDResponse.h"
  33 #include "TH1.h"
  34 #include "TH2.h"
  35 #include "AliStack.h"
  36 #include "AliAODConversionMother.h"
  37 #include "TObjString.h"
  38 #include "AliAODEvent.h"
  39 #include "AliESDEvent.h"
  40 #include "AliCentrality.h"
  41 #include "TList.h"
  42 #include "TPDGCode.h"
  43 #include "TDatabasePDG.h"
  44 #include "AliAODMCParticle.h"
  45
  46 class iostream;
  47
  48 using namespace std;
  49
  50 ClassImp(AliConversionMesonCuts)
  51
  52
  53 const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = {
  54    "MesonKind", //0
  55    "BackgroundScheme", //1
  56    "NumberOfBGEvents", //2
  57    "DegreesForRotationMethod", //3
  58    "RapidityMesonCut", //4
  59    "RCut", //5
  60    "AlphaMesonCut", //6
  61    "Chi2MesonCut", //7
  62    "SharedElectronCuts", //8
  63    "RejectToCloseV0s", //9
  64    "UseMCPSmearing", //10
  65    "DcaGammaGamma", //11
  66    "DcaRPrimVtx", //12
  67    "DcaZPrimVtx" //13
  68 };
  69
  70
  71 //________________________________________________________________________
  72 AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) :
  73    AliAnalysisCuts(name,title),
  74    fHistograms(NULL),
  75    fMesonKind(0),
  76    fMaxR(200),
  77    fChi2CutMeson(1000),
  78    fAlphaMinCutMeson(0),
  79    fAlphaCutMeson(1),
  80    fRapidityCutMeson(1),
  81    fUseRotationMethodInBG(kFALSE),
  82    fDoBG(kTRUE),
  83    fdoBGProbability(kFALSE),
  84    fUseTrackMultiplicityForBG(kFALSE),
  85    fnDegreeRotationPMForBG(0),
  86    fNumberOfBGEvents(0),
  87    fOpeningAngle(0.005),
  88    fDoToCloseV0sCut(kFALSE),
  89    fminV0Dist(200.),
  90    fDoSharedElecCut(kFALSE),
  91    fUseMCPSmearing(kFALSE),
  92    fPBremSmearing(0),
  93    fPSigSmearing(0),
  94    fPSigSmearingCte(0),
  95    fBrem(NULL),
  96    fRandom(0),
  97    fFAlphaCut(0),
  98    fAlphaPtDepCut(kFALSE),
  99    fElectronLabelArraySize(500),
 100    fElectronLabelArray(NULL),
 101    fDCAGammaGammaCut(1000),
 102    fDCAZMesonPrimVtxCut(1000),
 103    fDCARMesonPrimVtxCut(1000),
 104    fBackgroundHandler(0),
 105    fCutString(NULL),
 106    hMesonCuts(NULL),
 107    hMesonBGCuts(NULL),
 108    hDCAGammaGammaMesonBefore(NULL),
 109    hDCAZMesonPrimVtxBefore(NULL),
 110    hDCARMesonPrimVtxBefore(NULL),
 111    hDCAGammaGammaMesonAfter(NULL),
 112    hDCAZMesonPrimVtxAfter(NULL),
 113    hDCARMesonPrimVtxAfter(NULL)
 114
 115 {
 116    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
 117    fCutString=new TObjString((GetCutNumber()).Data());
 118    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
 119    if (fBrem == NULL){
 120       fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
 121       // tests done with 1.0e-14
 122       fBrem->SetParameter(0,fPBremSmearing);
 123       fBrem->SetNpx(100000);
 124    }
 125
 126 }
 127
 128 //________________________________________________________________________
 129 AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref) :
 130    AliAnalysisCuts(ref),
 131    fHistograms(NULL),
 132    fMesonKind(ref.fMesonKind),
 133    fMaxR(ref.fMaxR),
 134    fChi2CutMeson(ref.fChi2CutMeson),
 135    fAlphaMinCutMeson(ref.fAlphaMinCutMeson),
 136    fAlphaCutMeson(ref.fAlphaCutMeson),
 137    fRapidityCutMeson(ref.fRapidityCutMeson),
 138    fUseRotationMethodInBG(ref.fUseRotationMethodInBG),
 139    fDoBG(ref.fDoBG),
 140    fdoBGProbability(ref.fdoBGProbability),
 141    fUseTrackMultiplicityForBG(ref.fUseTrackMultiplicityForBG),
 142    fnDegreeRotationPMForBG(ref.fnDegreeRotationPMForBG),
 143    fNumberOfBGEvents(ref. fNumberOfBGEvents),
 144    fOpeningAngle(ref.fOpeningAngle),
 145    fDoToCloseV0sCut(ref.fDoToCloseV0sCut),
 146    fminV0Dist(ref.fminV0Dist),
 147    fDoSharedElecCut(ref.fDoSharedElecCut),
 148    fUseMCPSmearing(ref.fUseMCPSmearing),
 149    fPBremSmearing(ref.fPBremSmearing),
 150    fPSigSmearing(ref.fPSigSmearing),
 151    fPSigSmearingCte(ref.fPSigSmearingCte),
 152    fBrem(NULL),
 153    fRandom(ref.fRandom),
 154    fFAlphaCut(NULL),
 155    fAlphaPtDepCut(ref.fAlphaPtDepCut),
 156    fElectronLabelArraySize(ref.fElectronLabelArraySize),
 157    fElectronLabelArray(NULL),
 158    fDCAGammaGammaCut(ref.fDCAGammaGammaCut),
 159    fDCAZMesonPrimVtxCut(ref.fDCAZMesonPrimVtxCut),
 160    fDCARMesonPrimVtxCut(ref.fDCARMesonPrimVtxCut),
 161    fBackgroundHandler(ref.fBackgroundHandler),
 162    fCutString(NULL),
 163    hMesonCuts(NULL),
 164    hMesonBGCuts(NULL),
 165    hDCAGammaGammaMesonBefore(NULL),
 166    hDCAZMesonPrimVtxBefore(NULL),
 167    hDCARMesonPrimVtxBefore(NULL),
 168    hDCAGammaGammaMesonAfter(NULL),
 169    hDCAZMesonPrimVtxAfter(NULL),
 170    hDCARMesonPrimVtxAfter(NULL)
 171 {
 172    // Copy Constructor
 173    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
 174    fCutString=new TObjString((GetCutNumber()).Data());
 175    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
 176    if (fBrem == NULL)fBrem = (TF1*)ref.fBrem->Clone("fBrem");
 177    // Histograms are not copied, if you need them, call InitCutHistograms
 178 }
 179
 180
 181 //________________________________________________________________________
 182 AliConversionMesonCuts::~AliConversionMesonCuts() {
 183    // Destructor
 184    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
 185    // if(fHistograms)
 186    //   delete fHistograms;
 187    // fHistograms = NULL;
 188    if(fCutString != NULL){
 189       delete fCutString;
 190       fCutString = NULL;
 191    }
 192    if(fElectronLabelArray){
 193       delete fElectronLabelArray;
 194       fElectronLabelArray = NULL;
 195    }
 196
 197 }
 198
 199 //________________________________________________________________________
 200 void AliConversionMesonCuts::InitCutHistograms(TString name, Bool_t additionalHists){
 201
 202    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
 203    TH1::AddDirectory(kFALSE);
 204
 205    if(fHistograms != NULL){
 206       delete fHistograms;
 207       fHistograms=NULL;
 208    }
 209    if(fHistograms==NULL){
 210       fHistograms=new TList();
 211       fHistograms->SetOwner(kTRUE);
 212       if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
 213       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
 214    }
 215
 216    // Meson Cuts
 217    hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",13,-0.5,12.5);
 218    hMesonCuts->GetXaxis()->SetBinLabel(1,"in");
 219    hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 220    hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 221    hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 222    hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 223    hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 224    hMesonCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
 225    hMesonCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
 226    hMesonCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
 227    hMesonCuts->GetXaxis()->SetBinLabel(10,"out");
 228    fHistograms->Add(hMesonCuts);
 229
 230    hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",13,-0.5,12.5);
 231    hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in");
 232    hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 233    hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 234    hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 235    hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 236    hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 237    hMesonBGCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
 238    hMesonBGCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
 239    hMesonBGCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
 240    hMesonBGCuts->GetXaxis()->SetBinLabel(10,"out");
 241
 242    fHistograms->Add(hMesonBGCuts);
 243
 244    if (additionalHists){
 245       hDCAGammaGammaMesonBefore=new TH1F(Form("DCAGammaGammaMeson Before %s",GetCutNumber().Data()),"DCAGammaGammaMeson Before",200,0,10);
 246       fHistograms->Add(hDCAGammaGammaMesonBefore);
 247
 248       hDCARMesonPrimVtxBefore=new TH1F(Form("DCARMesonPrimVtx Before %s",GetCutNumber().Data()),"DCARMesonPrimVtx Before",200,0,10);
 249       fHistograms->Add(hDCARMesonPrimVtxBefore);
 250
 251       hDCAZMesonPrimVtxBefore=new TH1F(Form("DCAZMesonPrimVtx Before %s",GetCutNumber().Data()),"DCAZMesonPrimVtx Before",401,-10,10);
 252       fHistograms->Add(hDCAZMesonPrimVtxBefore);
 253
 254    }
 255
 256    hDCAGammaGammaMesonAfter=new TH1F(Form("DCAGammaGammaMeson After %s",GetCutNumber().Data()),"DCAGammaGammaMeson After",200,0,10);
 257    fHistograms->Add(hDCAGammaGammaMesonAfter);
 258
 259    hDCAZMesonPrimVtxAfter=new TH2F(Form("InvMassDCAZMesonPrimVtx After %s",GetCutNumber().Data()),"InvMassDCAZMesonPrimVtx After",800,0,0.8,401,-10,10);
 260    fHistograms->Add(hDCAZMesonPrimVtxAfter);
 261
 262    hDCARMesonPrimVtxAfter=new TH1F(Form("DCARMesonPrimVtx After %s",GetCutNumber().Data()),"DCARMesonPrimVtx After",200,0,10);
 263    fHistograms->Add(hDCARMesonPrimVtxAfter);
 264
 265    TH1::AddDirectory(kTRUE);
 266 }
 267
 268 //________________________________________________________________________
 269 Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack, Double_t fRapidityShift){
 270    // Returns true for all pions within acceptance cuts for decay into 2 photons
 271    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 272
 273    if(!fMCStack)return kFALSE;
 274
 275    if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
 276       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 277
 278       Double_t rapidity = 10.;
 279       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 280          rapidity=8.-fRapidityShift;
 281       } else{
 282          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 283       }
 284
 285       // Rapidity Cut
 286       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 287
 288       // Select only -> 2y decay channel
 289       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 290
 291       for(Int_t i=0;i<2;i++){
 292          TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
 293          // Is Daughter a Photon?
 294          if(MDaughter->GetPdgCode()!=22)return kFALSE;
 295          // Is Photon in Acceptance?
 296          //   if(bMCDaughtersInAcceptance){
 297          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
 298          //   }
 299       }
 300       return kTRUE;
 301    }
 302    return kFALSE;
 303 }
 304 //________________________________________________________________________
 305 Bool_t AliConversionMesonCuts::MesonIsSelectedAODMC(AliAODMCParticle *MCMother,TClonesArray *AODMCArray, Double_t fRapidityShift){
 306    // Returns true for all pions within acceptance cuts for decay into 2 photons
 307    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 308
 309    if(!AODMCArray)return kFALSE;
 310
 311    if(MCMother->GetPdgCode()==111 || MCMother->GetPdgCode()==221){
 312       Double_t rMeson = sqrt( (MCMother->Xv()*MCMother->Xv()) + (MCMother->Yv()*MCMother->Yv()) ) ;
 313       if(rMeson>fMaxR)  return kFALSE; // cuts on distance from collision point
 314
 315       Double_t rapidity = 10.;
 316       if(MCMother->E() - MCMother->Pz() == 0 || MCMother->E() + MCMother->Pz() == 0){
 317          rapidity=8.-fRapidityShift;
 318       } else{
 319          rapidity = 0.5*(TMath::Log((MCMother->E()+MCMother->Pz()) / (MCMother->E()-MCMother->Pz())))-fRapidityShift;
 320       }
 321
 322       // Rapidity Cut
 323       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 324
 325       // Select only -> 2y decay channel
 326       if(MCMother->GetNDaughters()!=2)return kFALSE;
 327
 328       for(Int_t i=0;i<2;i++){
 329          AliAODMCParticle *MDaughter=static_cast<AliAODMCParticle*>(AODMCArray->At(MCMother->GetDaughter(i)));
 330          // Is Daughter a Photon?
 331          if(MDaughter->GetPdgCode()!=22)return kFALSE;
 332          // Is Photon in Acceptance?
 333          //   if(bMCDaughtersInAcceptance){
 334          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
 335          //   }
 336       }
 337       return kTRUE;
 338    }
 339    return kFALSE;
 340 }
 341
 342 //________________________________________________________________________
 343 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelelectron, Int_t &labelpositron, Int_t &labelgamma, Double_t fRapidityShift){
 344
 345    // Returns true for all pions within acceptance cuts for decay into 2 photons
 346    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 347
 348    if( !fMCStack )return kFALSE;
 349
 350    if(  fMCMother->GetPdgCode() != 111 && fMCMother->GetPdgCode() != 221 ) return kFALSE;
 351
 352    if(  fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 353
 354    Double_t rapidity = 10.;
 355
 356    if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 357       rapidity=8.-fRapidityShift;
 358    }
 359    else{
 360       rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 361    }
 362
 363    // Rapidity Cut
 364    if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 365
 366    // Select only -> Dalitz decay channel
 367    if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 368
 369    TParticle *positron = 0x0;
 370    TParticle *electron = 0x0;
 371    TParticle    *gamma = 0x0;
 372
 373    for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 374
 375       TParticle* temp = (TParticle*)fMCStack->Particle( index );
 376
 377       switch( temp->GetPdgCode() ) {
 378       case ::kPositron:
 379          positron      =  temp;
 380          labelpositron = index;
 381          break;
 382       case ::kElectron:
 383          electron      =  temp;
 384          labelelectron = index;
 385          break;
 386       case ::kGamma:
 387          gamma         =  temp;
 388          labelgamma    = index;
 389          break;
 390       }
 391    }
 392
 393    if( positron && electron && gamma) return kTRUE;
 394    return kFALSE;
 395 }
 396
 397 //________________________________________________________________________
 398 Bool_t AliConversionMesonCuts::MesonIsSelectedMCEtaPiPlPiMiGamma(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelGamma, Double_t fRapidityShift){
 399
 400         // Returns true for all pions within acceptance cuts for decay into 2 photons
 401         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 402
 403         if( !fMCStack )return kFALSE;
 404
 405         if( fMCMother->GetPdgCode() != 221 ) return kFALSE;
 406
 407         if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 408
 409         Double_t rapidity = 10.;
 410
 411         if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 412                 rapidity=8.-fRapidityShift;
 413         }
 414         else{
 415                 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 416         }
 417
 418         // Rapidity Cut
 419         if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 420
 421         // Select only -> Dalitz decay channel
 422         if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 423
 424         TParticle *posPion = 0x0;
 425         TParticle *negPion = 0x0;
 426         TParticle    *gamma = 0x0;
 427
 428         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 429
 430                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
 431
 432                 switch( temp->GetPdgCode() ) {
 433                 case 211:
 434                         posPion      =  temp;
 435                         labelPosPion = index;
 436                         break;
 437                 case -211:
 438                         negPion      =  temp;
 439                         labelNegPion = index;
 440                         break;
 441                 case ::kGamma:
 442                         gamma         =  temp;
 443                         labelGamma    = index;
 444                         break;
 445                 }
 446         }
 447
 448         if( posPion && negPion && gamma) return kTRUE;
 449         return kFALSE;
 450 }
 451
 452 //________________________________________________________________________
 453 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
 454    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
 455    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 456
 457    if(!fMCStack)return kFALSE;
 458    // if(fMCMother->GetPdgCode()==20443 ){
 459    //    return kFALSE;
 460    // }
 461    if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
 462       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 463
 464       Double_t rapidity = 10.;
 465       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 466          rapidity=8.-fRapidityShift;
 467       }
 468       else{
 469          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 470       }
 471
 472       // Rapidity Cut
 473       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 474
 475       // Select only -> ChiC radiative (JPsi+gamma) decay channel
 476       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 477
 478       TParticle *jpsi   = 0x0;
 479       TParticle *gamma  = 0x0;
 480       TParticle *positron = 0x0;
 481       TParticle *electron = 0x0;
 482
 483       Int_t labeljpsiChiC = -1;
 484
 485       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 486
 487          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 488
 489          switch( temp->GetPdgCode() ) {
 490          case 443:
 491             jpsi =  temp;
 492             labeljpsiChiC = index;
 493             break;
 494          case 22:
 495             gamma    =  temp;
 496             labelgammaChiC = index;
 497             break;
 498          }
 499       }
 500
 501       if ( !jpsi || ! gamma) return kFALSE;
 502       if(jpsi->GetNDaughters()!=2)return kFALSE;
 503
 504
 505       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
 506          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 507          switch( temp->GetPdgCode() ) {
 508          case -11:
 509             electron =  temp;
 510             labelelectronChiC = index;
 511             break;
 512          case 11:
 513             positron =  temp;
 514             labelpositronChiC = index;
 515             break;
 516          }
 517       }
 518       if( !electron || !positron) return kFALSE;
 519       if( positron && electron && gamma) return kTRUE;
 520    }
 521    return kFALSE;
 522 }
 523
 524
 525
 526 ///________________________________________________________________________
 527 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
 528 {
 529
 530    // Selection of reconstructed Meson candidates
 531    // Use flag IsSignal in order to fill Fill different
 532    // histograms for Signal and Background
 533    TH1 *hist=0x0;
 534
 535    if(IsSignal){hist=hMesonCuts;}
 536    else{hist=hMesonBGCuts;}
 537
 538    Int_t cutIndex=0;
 539
 540    if(hist)hist->Fill(cutIndex);
 541    cutIndex++;
 542
 543    // Undefined Rapidity -> Floating Point exception
 544    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
 545       if(hist)hist->Fill(cutIndex);
 546       cutIndex++;
 547 //        cout << "undefined rapidity" << endl;
 548       return kFALSE;
 549    }
 550    else{
 551       // PseudoRapidity Cut --> But we cut on Rapidity !!!
 552       cutIndex++;
 553       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
 554          if(hist)hist->Fill(cutIndex);
 555 //               cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
 556          return kFALSE;
 557       }
 558    }
 559    cutIndex++;
 560
 561    // Opening Angle Cut
 562    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
 563    if( pi0->GetOpeningAngle() < fOpeningAngle){
 564 //        cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl;
 565       if(hist)hist->Fill(cutIndex);
 566       return kFALSE;
 567    }
 568    cutIndex++;
 569
 570    if ( fAlphaPtDepCut == kTRUE ) {
 571
 572         fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
 573    }
 574
 575
 576    // Alpha Max Cut
 577    if(pi0->GetAlpha()>fAlphaCutMeson){
 578 //         cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl;
 579       if(hist)hist->Fill(cutIndex);
 580       return kFALSE;
 581    }
 582    cutIndex++;
 583
 584    // Alpha Min Cut
 585    if(pi0->GetAlpha()<fAlphaMinCutMeson){
 586 //        cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl;
 587       if(hist)hist->Fill(cutIndex);
 588       return kFALSE;
 589    }
 590    cutIndex++;
 591
 592    if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
 593    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
 594
 595    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
 596 //        cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl;
 597       if(hist)hist->Fill(cutIndex);
 598       return kFALSE;
 599    }
 600    cutIndex++;
 601
 602    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
 603 //         cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl;
 604       if(hist)hist->Fill(cutIndex);
 605       return kFALSE;
 606    }
 607    cutIndex++;
 608
 609
 610    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
 611
 612    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
 613 //        cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl;
 614       if(hist)hist->Fill(cutIndex);
 615       return kFALSE;
 616    }
 617    cutIndex++;
 618
 619
 620    if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
 621    if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
 622    if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
 623
 624    if(hist)hist->Fill(cutIndex);
 625    return kTRUE;
 626 }
 627
 628
 629
 630 ///________________________________________________________________________
 631 ///________________________________________________________________________
 632 Bool_t AliConversionMesonCuts::UpdateCutString() {
 633    ///Update the cut string (if it has been created yet)
 634
 635    if(fCutString && fCutString->GetString().Length() == kNCuts) {
 636       fCutString->SetString(GetCutNumber());
 637    } else {
 638       return kFALSE;
 639    }
 640    return kTRUE;
 641 }
 642
 643 ///________________________________________________________________________
 644 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 645    // Initialize Cuts from a given Cut string
 646    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
 647    if(analysisCutSelection.Length()!=kNCuts) {
 648       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 649       return kFALSE;
 650    }
 651    if(!analysisCutSelection.IsDigit()){
 652       AliError("Cut selection contains characters");
 653       return kFALSE;
 654    }
 655
 656    const char *cutSelection = analysisCutSelection.Data();
 657 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 658    for(Int_t ii=0;ii<kNCuts;ii++){
 659       ASSIGNARRAY(ii);
 660    }
 661
 662    // Set Individual Cuts
 663    for(Int_t ii=0;ii<kNCuts;ii++){
 664       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 665    }
 666
 667    //PrintCuts();
 668    return kTRUE;
 669 }
 670 ///________________________________________________________________________
 671 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
 672    ///Set individual cut ID
 673
 674    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 675    switch (cutID) {
 676    case kMesonKind:
 677       if( SetMesonKind(value)) {
 678          fCuts[kMesonKind] = value;
 679          UpdateCutString();
 680          return kTRUE;
 681       } else return kFALSE;
 682    case kchi2MesonCut:
 683       if( SetChi2MesonCut(value)) {
 684          fCuts[kchi2MesonCut] = value;
 685          UpdateCutString();
 686          return kTRUE;
 687       } else return kFALSE;
 688    case kalphaMesonCut:
 689       if( SetAlphaMesonCut(value)) {
 690          fCuts[kalphaMesonCut] = value;
 691          UpdateCutString();
 692          return kTRUE;
 693       } else return kFALSE;
 694    case kRCut:
 695       if( SetRCut(value)) {
 696          fCuts[kRCut] = value;
 697          UpdateCutString();
 698          return kTRUE;
 699       } else return kFALSE;
 700
 701    case kRapidityMesonCut:
 702       if( SetRapidityMesonCut(value)) {
 703          fCuts[kRapidityMesonCut] = value;
 704          UpdateCutString();
 705          return kTRUE;
 706       } else return kFALSE;
 707
 708    case kBackgroundScheme:
 709       if( SetBackgroundScheme(value)) {
 710          fCuts[kBackgroundScheme] = value;
 711          UpdateCutString();
 712          return kTRUE;
 713       } else return kFALSE;
 714
 715    case kDegreesForRotationMethod:
 716       if( SetNDegreesForRotationMethod(value)) {
 717          fCuts[kDegreesForRotationMethod] = value;
 718          UpdateCutString();
 719          return kTRUE;
 720       } else return kFALSE;
 721
 722    case kNumberOfBGEvents:
 723       if( SetNumberOfBGEvents(value)) {
 724          fCuts[kNumberOfBGEvents] = value;
 725          UpdateCutString();
 726          return kTRUE;
 727       } else return kFALSE;
 728
 729    case kuseMCPSmearing:
 730       if( SetMCPSmearing(value)) {
 731          fCuts[kuseMCPSmearing] = value;
 732          UpdateCutString();
 733          return kTRUE;
 734       } else return kFALSE;
 735    case kElecShare:
 736       if( SetSharedElectronCut(value)) {
 737          fCuts[kElecShare] = value;
 738          UpdateCutString();
 739          return kTRUE;
 740       } else return kFALSE;
 741    case kToCloseV0s:
 742       if( SetToCloseV0sCut(value)) {
 743          fCuts[kToCloseV0s] = value;
 744          UpdateCutString();
 745          return kTRUE;
 746       } else return kFALSE;
 747    case kDcaGammaGamma:
 748       if( SetDCAGammaGammaCut(value)) {
 749          fCuts[kDcaGammaGamma] = value;
 750          UpdateCutString();
 751          return kTRUE;
 752       } else return kFALSE;
 753    case kDcaZPrimVtx:
 754       if( SetDCAZMesonPrimVtxCut(value)) {
 755          fCuts[kDcaZPrimVtx] = value;
 756          UpdateCutString();
 757          return kTRUE;
 758       } else return kFALSE;
 759    case kDcaRPrimVtx:
 760       if( SetDCARMesonPrimVtxCut(value)) {
 761          fCuts[kDcaRPrimVtx] = value;
 762          UpdateCutString();
 763          return kTRUE;
 764       } else return kFALSE;
 765
 766    case kNCuts:
 767       cout << "Error:: Cut id out of range"<< endl;
 768       return kFALSE;
 769    }
 770
 771    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 772    return kFALSE;
 773
 774 }
 775
 776
 777 ///________________________________________________________________________
 778 void AliConversionMesonCuts::PrintCuts() {
 779    // Print out current Cut Selection
 780    for(Int_t ic = 0; ic < kNCuts; ic++) {
 781       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 782    }
 783 }
 784
 785 ///________________________________________________________________________
 786 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
 787    // Set Cut
 788    switch(mesonKind){
 789    case 0:
 790       fMesonKind=0;
 791       break;
 792    case 1:
 793       fMesonKind=1;;
 794       break;
 795    default:
 796       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
 797       return kFALSE;
 798    }
 799    return kTRUE;
 800 }
 801
 802 ///________________________________________________________________________
 803 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
 804    // Set Cut
 805    switch(rCut){
 806    case 0:
 807       fMaxR = 180.;
 808       break;
 809    case 1:
 810       fMaxR = 180.;
 811       break;
 812    case 2:
 813       fMaxR = 180.;
 814       break;
 815    case 3:
 816       fMaxR = 70.;
 817       break;
 818    case 4:
 819       fMaxR = 70.;
 820       break;
 821    case 5:
 822       fMaxR = 180.;
 823       break;
 824       // High purity cuts for PbPb
 825    case 9:
 826       fMaxR = 180.;
 827       break;
 828    default:
 829       cout<<"Warning: rCut not defined"<<rCut<<endl;
 830       return kFALSE;
 831    }
 832    return kTRUE;
 833 }
 834
 835 ///________________________________________________________________________
 836 Bool_t AliConversionMesonCuts::SetChi2MesonCut(Int_t chi2MesonCut){
 837    // Set Cut
 838    switch(chi2MesonCut){
 839    case 0:  // 100.
 840       fChi2CutMeson = 100.;
 841       break;
 842    case 1:  // 50.
 843       fChi2CutMeson = 50.;
 844       break;
 845    case 2:  // 30.
 846       fChi2CutMeson = 30.;
 847       break;
 848    case 3:
 849       fChi2CutMeson = 200.;
 850       break;
 851    case 4:
 852       fChi2CutMeson = 500.;
 853       break;
 854    case 5:
 855       fChi2CutMeson = 1000.;
 856       break;
 857    default:
 858       cout<<"Warning: Chi2MesonCut not defined "<<chi2MesonCut<<endl;
 859       return kFALSE;
 860    }
 861    return kTRUE;
 862 }
 863
 864
 865 ///________________________________________________________________________
 866 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 867 {   // Set Cut
 868    switch(alphaMesonCut){
 869    case 0:      // 0- 0.7
 870       fAlphaMinCutMeson  = 0.0;
 871       fAlphaCutMeson     = 0.7;
 872       fAlphaPtDepCut = kFALSE;
 873       break;
 874    case 1:      // Updated 31 October 2013 before 0.0 - 0.5
 875       if( fFAlphaCut ) delete fFAlphaCut;
 876       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 877       fFAlphaCut->SetParameter(0,0.7);
 878       fFAlphaCut->SetParameter(1,1.2);
 879       fAlphaMinCutMeson  =  0.0;
 880       fAlphaCutMeson     = -1.0;
 881       fAlphaPtDepCut = kTRUE;
 882       break;
 883    case 2:      // Updated 31 October 2013 before 0.5-1
 884       if( fFAlphaCut ) delete fFAlphaCut;
 885       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 886       fFAlphaCut->SetParameter(0,0.8);
 887       fFAlphaCut->SetParameter(1,1.2);
 888       fAlphaMinCutMeson  =  0.0;
 889       fAlphaCutMeson     = -1.0;
 890       fAlphaPtDepCut = kTRUE;
 891       break;
 892    case 3:      // 0.0-1
 893       fAlphaMinCutMeson  = 0.0;
 894       fAlphaCutMeson     = 1.;
 895       fAlphaPtDepCut = kFALSE;
 896       break;
 897    case 4:      // 0-0.65
 898       fAlphaMinCutMeson  = 0.0;
 899       fAlphaCutMeson     = 0.65;
 900       fAlphaPtDepCut = kFALSE;
 901       break;
 902    case 5:      // 0-0.75
 903       fAlphaMinCutMeson  = 0.0;
 904       fAlphaCutMeson     = 0.75;
 905       fAlphaPtDepCut = kFALSE;
 906       break;
 907    case 6:      // 0-0.8
 908       fAlphaMinCutMeson  = 0.0;
 909       fAlphaCutMeson     = 0.8;
 910       fAlphaPtDepCut = kFALSE;
 911       break;
 912    case 7:      // 0.0-0.85
 913       fAlphaMinCutMeson  = 0.0;
 914       fAlphaCutMeson     = 0.85;
 915       fAlphaPtDepCut = kFALSE;
 916       break;
 917    case 8:      // 0.0-0.6
 918       fAlphaMinCutMeson  = 0.0;
 919       fAlphaCutMeson     = 0.6;
 920       fAlphaPtDepCut = kFALSE;
 921       break;
 922    case 9: // Updated 11 November 2013 before 0.0 - 0.3
 923       if( fFAlphaCut ) delete fFAlphaCut;
 924       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 925       fFAlphaCut->SetParameter(0,0.65);
 926       fFAlphaCut->SetParameter(1,1.2);
 927       fAlphaMinCutMeson  =  0.0;
 928       fAlphaCutMeson     = -1.0;
 929       fAlphaPtDepCut = kTRUE;
 930       break;
 931    default:
 932       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
 933       return kFALSE;
 934    }
 935    return kTRUE;
 936 }
 937
 938 ///________________________________________________________________________
 939 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
 940    // Set Cut
 941    switch(RapidityMesonCut){
 942    case 0:  // changed from 0.9 to 1.35
 943       fRapidityCutMeson   = 1.35;
 944       break;
 945    case 1:  //
 946       fRapidityCutMeson   = 0.8;
 947       break;
 948    case 2:  //
 949       fRapidityCutMeson   = 0.7;
 950       break;
 951    case 3:  //
 952       fRapidityCutMeson   = 0.6;
 953       break;
 954    case 4:  //
 955       fRapidityCutMeson   = 0.5;
 956       break;
 957    case 5:  //
 958       fRapidityCutMeson   = 0.85;
 959       break;
 960    case 6:  //
 961       fRapidityCutMeson   = 0.75;
 962       break;
 963    case 7:  //
 964       fRapidityCutMeson   = 0.3;
 965       break;
 966    case 8:  //
 967       fRapidityCutMeson   = 0.35;
 968       break;
 969    case 9:  //
 970       fRapidityCutMeson   = 0.4;
 971       break;
 972    default:
 973       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
 974       return kFALSE;
 975    }
 976    return kTRUE;
 977 }
 978
 979
 980 ///________________________________________________________________________
 981 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
 982    // Set Cut
 983    switch(BackgroundScheme){
 984    case 0: //Rotation
 985       fUseRotationMethodInBG=kTRUE;
 986       fdoBGProbability=kFALSE;
 987       break;
 988    case 1: // mixed event with V0 multiplicity
 989       fUseRotationMethodInBG=kFALSE;
 990       fUseTrackMultiplicityForBG=kFALSE;
 991       fdoBGProbability=kFALSE;
 992       break;
 993    case 2: // mixed event with track multiplicity
 994       fUseRotationMethodInBG=kFALSE;
 995       fUseTrackMultiplicityForBG=kTRUE;
 996       fdoBGProbability=kFALSE;
 997       break;
 998    case 3: //Rotation
 999       fUseRotationMethodInBG=kTRUE;
1000       fdoBGProbability=kTRUE;
1001       break;
1002    case 4: //No BG calculation
1003       cout << "no BG calculation should be done" << endl;
1004       fUseRotationMethodInBG=kFALSE;
1005       fdoBGProbability=kFALSE;
1006       fDoBG=kFALSE;
1007       break;
1008    case 5: //Rotation
1009       fUseRotationMethodInBG=kTRUE;
1010       fdoBGProbability=kFALSE;
1011       fBackgroundHandler = 1;
1012       break;
1013    case 6: // mixed event with V0 multiplicity
1014       fUseRotationMethodInBG=kFALSE;
1015       fUseTrackMultiplicityForBG=kFALSE;
1016       fdoBGProbability=kFALSE;
1017       fBackgroundHandler = 1;
1018       break;
1019    case 7: // mixed event with track multiplicity
1020       fUseRotationMethodInBG=kFALSE;
1021       fUseTrackMultiplicityForBG=kTRUE;
1022       fdoBGProbability=kFALSE;
1023       fBackgroundHandler = 1;
1024       break;
1025    case 8: //Rotation
1026       fUseRotationMethodInBG=kTRUE;
1027       fdoBGProbability=kTRUE;
1028       fBackgroundHandler = 1;
1029       break;
1030    default:
1031       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1032       return kFALSE;
1033    }
1034    return kTRUE;
1035 }
1036
1037
1038 ///________________________________________________________________________
1039 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
1040    // Set Cut
1041    switch(DegreesForRotationMethod){
1042    case 0:
1043       fnDegreeRotationPMForBG = 5;
1044       break;
1045    case 1:
1046       fnDegreeRotationPMForBG = 10;
1047       break;
1048    case 2:
1049       fnDegreeRotationPMForBG = 15;
1050       break;
1051    case 3:
1052       fnDegreeRotationPMForBG = 20;
1053       break;
1054    default:
1055       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1056       return kFALSE;
1057    }
1058    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1059    return kTRUE;
1060 }
1061
1062 ///________________________________________________________________________
1063 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
1064    // Set Cut
1065    switch(NumberOfBGEvents){
1066    case 0:
1067       fNumberOfBGEvents = 5;
1068       break;
1069    case 1:
1070       fNumberOfBGEvents = 10;
1071       break;
1072    case 2:
1073       fNumberOfBGEvents = 15;
1074       break;
1075    case 3:
1076       fNumberOfBGEvents = 20;
1077       break;
1078    case 4:
1079       fNumberOfBGEvents = 2;
1080       break;
1081    case 5:
1082       fNumberOfBGEvents = 50;
1083       break;
1084    case 6:
1085       fNumberOfBGEvents = 80;
1086       break;
1087    case 7:
1088       fNumberOfBGEvents = 100;
1089       break;
1090    default:
1091       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1092       return kFALSE;
1093    }
1094    return kTRUE;
1095 }
1096 ///________________________________________________________________________
1097 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1098
1099    switch(sharedElec){
1100    case 0:
1101       fDoSharedElecCut = kFALSE;
1102       break;
1103    case 1:
1104       fDoSharedElecCut = kTRUE;
1105       break;
1106    default:
1107       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1108       return kFALSE;
1109    }
1110
1111    return kTRUE;
1112 }
1113
1114 ///________________________________________________________________________
1115 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1116
1117    switch(toClose){
1118    case 0:
1119       fDoToCloseV0sCut = kFALSE;
1120       fminV0Dist = 250;
1121       break;
1122    case 1:
1123       fDoToCloseV0sCut = kTRUE;
1124       fminV0Dist = 1;
1125       break;
1126    case 2:
1127       fDoToCloseV0sCut = kTRUE;
1128       fminV0Dist = 2;
1129       break;
1130    case 3:
1131       fDoToCloseV0sCut = kTRUE;
1132       fminV0Dist = 3;
1133       break;
1134    default:
1135       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1136       return kFALSE;
1137    }
1138    return kTRUE;
1139 }
1140
1141 ///________________________________________________________________________
1142 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1143 {// Set Cut
1144    switch(useMCPSmearing){
1145    case 0:
1146       fUseMCPSmearing=0;
1147       fPBremSmearing=1.;
1148       fPSigSmearing=0.;
1149       fPSigSmearingCte=0.;
1150       break;
1151    case 1:
1152       fUseMCPSmearing=1;
1153       fPBremSmearing=1.0e-14;
1154       fPSigSmearing=0.;
1155       fPSigSmearingCte=0.;
1156       break;
1157    case 2:
1158       fUseMCPSmearing=1;
1159       fPBremSmearing=1.0e-15;
1160       fPSigSmearing=0.0;
1161       fPSigSmearingCte=0.;
1162       break;
1163    case 3:
1164       fUseMCPSmearing=1;
1165       fPBremSmearing=1.;
1166       fPSigSmearing=0.003;
1167       fPSigSmearingCte=0.002;
1168       break;
1169    case 4:
1170       fUseMCPSmearing=1;
1171       fPBremSmearing=1.;
1172       fPSigSmearing=0.003;
1173       fPSigSmearingCte=0.007;
1174       break;
1175    case 5:
1176       fUseMCPSmearing=1;
1177       fPBremSmearing=1.;
1178       fPSigSmearing=0.003;
1179       fPSigSmearingCte=0.016;
1180       break;
1181    case 6:
1182       fUseMCPSmearing=1;
1183       fPBremSmearing=1.;
1184       fPSigSmearing=0.007;
1185       fPSigSmearingCte=0.016;
1186       break;
1187    case 7:
1188       fUseMCPSmearing=1;
1189       fPBremSmearing=1.0e-16;
1190       fPSigSmearing=0.0;
1191       fPSigSmearingCte=0.;
1192       break;
1193    case 8:
1194       fUseMCPSmearing=1;
1195       fPBremSmearing=1.;
1196       fPSigSmearing=0.007;
1197       fPSigSmearingCte=0.014;
1198       break;
1199    case 9:
1200       fUseMCPSmearing=1;
1201       fPBremSmearing=1.;
1202       fPSigSmearing=0.007;
1203       fPSigSmearingCte=0.011;
1204       break;
1205
1206    default:
1207       AliError("Warning: UseMCPSmearing not defined");
1208       return kFALSE;
1209    }
1210    return kTRUE;
1211 }
1212
1213
1214 ///________________________________________________________________________
1215 Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1216    // Set Cut
1217    switch(DCAGammaGamma){
1218    case 0:  //
1219       fDCAGammaGammaCut   = 1000;
1220       break;
1221    case 1:  //
1222       fDCAGammaGammaCut   = 10;
1223       break;
1224    case 2:  //
1225       fDCAGammaGammaCut   = 5;
1226       break;
1227    case 3:  //
1228       fDCAGammaGammaCut   = 4;
1229       break;
1230    case 4:  //
1231       fDCAGammaGammaCut   = 3;
1232       break;
1233    case 5:  //
1234       fDCAGammaGammaCut   = 2.5;
1235       break;
1236    case 6:  //
1237       fDCAGammaGammaCut   = 2;
1238       break;
1239    case 7:  //
1240       fDCAGammaGammaCut   = 1.5;
1241       break;
1242    case 8:  //
1243       fDCAGammaGammaCut   = 1;
1244       break;
1245    case 9:  //
1246       fDCAGammaGammaCut   = 0.5;
1247       break;
1248    default:
1249       cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1250       return kFALSE;
1251    }
1252    return kTRUE;
1253 }
1254
1255 ///________________________________________________________________________
1256 Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1257    // Set Cut
1258    switch(DCAZMesonPrimVtx){
1259    case 0:  //
1260       fDCAZMesonPrimVtxCut   = 1000;
1261       break;
1262    case 1:  //
1263       fDCAZMesonPrimVtxCut   = 10;
1264       break;
1265    case 2:  //
1266       fDCAZMesonPrimVtxCut   = 5;
1267       break;
1268    case 3:  //
1269       fDCAZMesonPrimVtxCut   = 4;
1270       break;
1271    case 4:  //
1272       fDCAZMesonPrimVtxCut   = 3;
1273       break;
1274    case 5:  //
1275       fDCAZMesonPrimVtxCut   = 2.5;
1276       break;
1277    case 6:  //
1278       fDCAZMesonPrimVtxCut   = 2;
1279       break;
1280    case 7:  //
1281       fDCAZMesonPrimVtxCut   = 1.5;
1282       break;
1283    case 8:  //
1284       fDCAZMesonPrimVtxCut   = 1;
1285       break;
1286    case 9:  //
1287       fDCAZMesonPrimVtxCut   = 0.5;
1288       break;
1289    default:
1290       cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1291       return kFALSE;
1292    }
1293    return kTRUE;
1294 }
1295
1296 ///________________________________________________________________________
1297 Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1298    // Set Cut
1299    switch(DCARMesonPrimVtx){
1300    case 0:  //
1301       fDCARMesonPrimVtxCut   = 1000;
1302       break;
1303    case 1:  //
1304       fDCARMesonPrimVtxCut   = 10;
1305       break;
1306    case 2:  //
1307       fDCARMesonPrimVtxCut   = 5;
1308       break;
1309    case 3:  //
1310       fDCARMesonPrimVtxCut   = 4;
1311       break;
1312    case 4:  //
1313       fDCARMesonPrimVtxCut   = 3;
1314       break;
1315    case 5:  //
1316       fDCARMesonPrimVtxCut   = 2.5;
1317       break;
1318    case 6:  //
1319       fDCARMesonPrimVtxCut   = 2;
1320       break;
1321    case 7:  //
1322       fDCARMesonPrimVtxCut   = 1.5;
1323       break;
1324    case 8:  //
1325       fDCARMesonPrimVtxCut   = 1;
1326       break;
1327    case 9:  //
1328       fDCARMesonPrimVtxCut   = 0.5;
1329       break;
1330    default:
1331       cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1332       return kFALSE;
1333    }
1334    return kTRUE;
1335 }
1336
1337
1338 ///________________________________________________________________________
1339 TString AliConversionMesonCuts::GetCutNumber(){
1340    // returns TString with current cut number
1341    TString a(kNCuts);
1342    for(Int_t ii=0;ii<kNCuts;ii++){
1343       a.Append(Form("%d",fCuts[ii]));
1344    }
1345    return a;
1346 }
1347
1348 ///________________________________________________________________________
1349 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
1350
1351    Int_t posLabel = photon->GetTrackLabelPositive();
1352    Int_t negLabel = photon->GetTrackLabelNegative();
1353
1354    fElectronLabelArray[nV0*2] = posLabel;
1355    fElectronLabelArray[(nV0*2)+1] = negLabel;
1356 }
1357
1358 ///________________________________________________________________________
1359 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1360
1361    Int_t posLabel = photon->GetTrackLabelPositive();
1362    Int_t negLabel = photon->GetTrackLabelNegative();
1363
1364    for(Int_t i = 0; i<nV0s*2;i++){
1365       if(i==nV0*2)     continue;
1366       if(i==(nV0*2)+1) continue;
1367       if(fElectronLabelArray[i] == posLabel){
1368          return kFALSE;}
1369       if(fElectronLabelArray[i] == negLabel){
1370          return kFALSE;}
1371    }
1372
1373    return kTRUE;
1374 }
1375 ///________________________________________________________________________
1376 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
1377    Double_t posX = photon->GetConversionX();
1378    Double_t posY = photon->GetConversionY();
1379    Double_t posZ = photon->GetConversionZ();
1380
1381    for(Int_t i = 0;i<photons->GetEntries();i++){
1382       if(nV0 == i) continue;
1383       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1384       Double_t posCompX = photonComp->GetConversionX();
1385       Double_t posCompY = photonComp->GetConversionY();
1386       Double_t posCompZ = photonComp->GetConversionZ();
1387
1388       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
1389
1390       if(dist < fminV0Dist*fminV0Dist){
1391          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1392          else {
1393             return kFALSE;}
1394       }
1395
1396    }
1397    return kTRUE;
1398 }
1399
1400 ///________________________________________________________________________
1401 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1402 {
1403
1404    if (photon==NULL) return;
1405    Double_t facPBrem = 1.;
1406    Double_t facPSig = 0.;
1407
1408    Double_t phi=0.;
1409    Double_t theta=0.;
1410    Double_t P=0.;
1411
1412
1413    P=photon->P();
1414    phi=photon->Phi();
1415    if( photon->P()!=0){
1416       theta=acos( photon->Pz()/ photon->P());
1417    }
1418
1419    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1420       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1421    }
1422
1423    if( fPBremSmearing != 1.){
1424       if(fBrem!=NULL){
1425          facPBrem = fBrem->GetRandom();
1426       }
1427    }
1428
1429    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1430    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1431    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1432    photon->SetE(photon->P());
1433 }