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 class iostream;
  43
  44 using namespace std;
  45
  46 ClassImp(AliConversionMesonCuts)
  47
  48
  49 const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = {
  50 "MesonKind",
  51 "BackgroundScheme",
  52 "NumberOfBGEvents",
  53 "DegreesForRotationMethod",
  54 "RapidityMesonCut",
  55 "RCut",
  56 "AlphaMesonCut",
  57 "Chi2MesonCut",
  58 "SharedElectronCuts",
  59 "RejectToCloseV0s",
  60 "UseMCPSmearing",
  61 };
  62
  63
  64 //________________________________________________________________________
  65 AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) : AliAnalysisCuts(name,title),
  66     fHistograms(NULL),
  67     fMaxR(200),
  68     fChi2CutMeson(1000),
  69     fAlphaMinCutMeson(0),
  70     fAlphaCutMeson(1),
  71     fRapidityCutMeson(1),
  72     fUseRotationMethodInBG(kFALSE),
  73     fdoBGProbability(kFALSE),
  74     fUseTrackMultiplicityForBG(kFALSE),
  75     fnDegreeRotationPMForBG(0),
  76     fNumberOfBGEvents(0),
  77     fOpeningAngle(0.005),
  78     fDoToCloseV0sCut(kFALSE),
  79     fminV0Dist(200.),
  80     fDoSharedElecCut(kFALSE),
  81     fUseMCPSmearing(kFALSE),
  82     fPBremSmearing(0),
  83     fPSigSmearing(0),
  84     fPSigSmearingCte(0),
  85     fBrem(NULL),
  86     fRandom(0),
  87     fElectronLabelArray(NULL),
  88     fCutString(NULL),
  89     hMesonCuts(NULL),
  90     hMesonBGCuts(NULL)
  91 {
  92         for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
  93         fCutString=new TObjString((GetCutNumber()).Data());
  94         fElectronLabelArray = new Int_t[500];
  95         if (fBrem == NULL){
  96                 fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
  97                 // tests done with 1.0e-14
  98                 fBrem->SetParameter(0,fPBremSmearing);
  99                 fBrem->SetNpx(100000);
 100         }
 101
 102 }
 103
 104 //________________________________________________________________________
 105 AliConversionMesonCuts::~AliConversionMesonCuts() {
 106     // Destructor
 107   //Deleting fHistograms leads to seg fault it it's added to output collection of a task
 108   // if(fHistograms)
 109   //    delete fHistograms;
 110   // fHistograms = NULL;
 111    if(fCutString != NULL){
 112       delete fCutString;
 113       fCutString = NULL;
 114    }
 115    if(fElectronLabelArray){
 116       delete fElectronLabelArray;
 117       fElectronLabelArray = NULL;
 118    }
 119
 120 }
 121
 122 //________________________________________________________________________
 123 void AliConversionMesonCuts::InitCutHistograms(TString name, Bool_t preCut){
 124
 125     // Initialize Cut Histograms for QA (only initialized and filled if function is called)
 126
 127         if(fHistograms != NULL){
 128                 delete fHistograms;
 129                 fHistograms=NULL;
 130    }
 131         if(fHistograms==NULL){
 132                 fHistograms=new TList();
 133                 if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
 134                 else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
 135     }
 136
 137     // Meson Cuts
 138     hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",10,-0.5,9.5);
 139     hMesonCuts->GetXaxis()->SetBinLabel(1,"in");
 140     hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 141     hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 142     hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 143     hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 144     hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 145     hMesonCuts->GetXaxis()->SetBinLabel(7,"out");
 146     fHistograms->Add(hMesonCuts);
 147
 148     hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",10,-0.5,9.5);
 149     hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in");
 150     hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 151     hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 152     hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 153     hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 154     hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 155     hMesonBGCuts->GetXaxis()->SetBinLabel(7,"out");
 156     fHistograms->Add(hMesonBGCuts);
 157 }
 158
 159
 160 //________________________________________________________________________
 161 Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack,Bool_t bMCDaughtersInAcceptance){
 162     // Returns true for all pions within acceptance cuts for decay into 2 photons
 163     // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 164
 165         if(!fMCStack)return kFALSE;
 166
 167         if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
 168                 if(fMCMother->R()>fMaxR)        return kFALSE; // cuts on distance from collision point
 169
 170                 Double_t rapidity = 10.;
 171                 if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 172                         rapidity=8.;
 173                 } else{
 174                         rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
 175                 }
 176
 177                 // Rapidity Cut
 178                 if(TMath::Abs(rapidity)>fRapidityCutMeson)return kFALSE;
 179
 180                 // Select only -> 2y decay channel
 181                 if(fMCMother->GetNDaughters()!=2)return kFALSE;
 182
 183                 for(Int_t i=0;i<2;i++){
 184                         TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
 185                         // Is Daughter a Photon?
 186                         if(MDaughter->GetPdgCode()!=22)return kFALSE;
 187                                         // Is Photon in Acceptance?
 188         //          if(bMCDaughtersInAcceptance){
 189         //              if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
 190         //          }
 191                 }
 192                 return kTRUE;
 193         }
 194         return kFALSE;
 195 }
 196
 197 //________________________________________________________________________
 198 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack){
 199         // Returns true for all pions within acceptance cuts for decay into 2 photons
 200         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 201
 202         if(!fMCStack)return kFALSE;
 203
 204         if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
 205
 206                 if(fMCMother->R()>fMaxR)        return kFALSE; // cuts on distance from collision point
 207
 208                 Double_t rapidity = 10.;
 209                 if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 210                         rapidity=8.;
 211                 }
 212                 else{
 213                         rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())));
 214                 }
 215
 216                 // Rapidity Cut
 217                 if(TMath::Abs(rapidity)>fRapidityCutMeson)return kFALSE;
 218
 219                 // Select only -> Dalitz decay channel
 220                 if(fMCMother->GetNDaughters()!=3)return kFALSE;
 221
 222                 Int_t daughterPDGs[3] = {0,0,0};
 223                 Int_t index = 0;
 224
 225 //                iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter()
 226
 227                 for(Int_t i=fMCMother->GetFirstDaughter(); i<= fMCMother->GetLastDaughter();i++){
 228                         TParticle *MDaughter=fMCStack->Particle(i);
 229                         // Is Daughter a Photon or an electron?
 230                         daughterPDGs[index]=MDaughter->GetPdgCode();
 231                         index++;
 232                 }
 233                 for (Int_t j=0;j<2;j++){
 234
 235                         for (Int_t i=0;i<2;i++){
 236                                 if (daughterPDGs[i] > daughterPDGs[i+1]){
 237                                         Int_t interMed = daughterPDGs[i] ;
 238                                         daughterPDGs[i] = daughterPDGs[i+1];
 239                                         daughterPDGs[i+1] = interMed;
 240                                 }
 241                         }
 242                 }
 243                 if (daughterPDGs[0] == -11 && daughterPDGs[1] == 11 && daughterPDGs[2] == 22) return kTRUE;
 244         }
 245         return kFALSE;
 246 }
 247
 248
 249 ///________________________________________________________________________
 250 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal)
 251 {
 252         // Selection of reconstructed Meson candidates
 253         // Use flag IsSignal in order to fill Fill different
 254         // histograms for Signal and Background
 255         TH1 *hist=0x0;
 256
 257
 258
 259         if(IsSignal){hist=hMesonCuts;}
 260         else{hist=hMesonBGCuts;}
 261
 262         Int_t cutIndex=0;
 263
 264         if(hist)hist->Fill(cutIndex);
 265         cutIndex++;
 266
 267         // Undefined Rapidity -> Floating Point exception
 268         if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
 269                 if(hist)hist->Fill(cutIndex);
 270                 cutIndex++;
 271                 return kFALSE;
 272         }
 273         else{
 274                 // PseudoRapidity Cut --> But we cut on Rapidity !!!
 275                 cutIndex++;
 276                 if(TMath::Abs(pi0->Rapidity())>fRapidityCutMeson){
 277                 //if(TMath::Abs(pi0->PseudoRapidity())>fRapidityCutMeson){
 278                         if(hist)hist->Fill(cutIndex);
 279                         return kFALSE;
 280                 }
 281         }
 282         cutIndex++;
 283
 284         // Opening Angle Cut
 285         //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
 286         if( pi0->GetOpeningAngle() < fOpeningAngle){
 287                 if(hist)hist->Fill(cutIndex);
 288                 return kFALSE;
 289         }
 290         cutIndex++;
 291
 292         // Alpha Max Cut
 293         if(pi0->GetAlpha()>fAlphaCutMeson){
 294                 if(hist)hist->Fill(cutIndex);
 295                 return kFALSE;
 296         }
 297         cutIndex++;
 298
 299         // Alpha Min Cut
 300         if(pi0->GetAlpha()<fAlphaMinCutMeson){
 301                 if(hist)hist->Fill(cutIndex);
 302                 return kFALSE;
 303         }
 304         cutIndex++;
 305
 306         if(hist)hist->Fill(cutIndex);
 307         return kTRUE;
 308 }
 309
 310
 311
 312 ///________________________________________________________________________
 313 ///________________________________________________________________________
 314 Bool_t AliConversionMesonCuts::UpdateCutString(cutIds cutID, Int_t value) {
 315 ///Update the cut string (if it has been created yet)
 316
 317   if(fCutString && fCutString->GetString().Length() == kNCuts) {
 318 //         cout << "Updating cut id in spot number " << cutID << " to " << value << endl;
 319         fCutString->SetString(GetCutNumber());
 320   } else {
 321 //         cout << "fCutString not yet initialized, will not be updated" << endl;
 322         return kFALSE;
 323   }
 324 //   cout << fCutString->GetString().Data() << endl;
 325   return kTRUE;
 326 }
 327
 328 ///________________________________________________________________________
 329 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 330         // Initialize Cuts from a given Cut string
 331         AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
 332         if(analysisCutSelection.Length()!=kNCuts) {
 333                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 334                 return kFALSE;
 335         }
 336         if(!analysisCutSelection.IsDigit()){
 337                 AliError("Cut selection contains characters");
 338                 return kFALSE;
 339         }
 340
 341         const char *cutSelection = analysisCutSelection.Data();
 342         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 343         for(Int_t ii=0;ii<kNCuts;ii++){
 344                         ASSIGNARRAY(ii);
 345         }
 346
 347   // Set Individual Cuts
 348   for(Int_t ii=0;ii<kNCuts;ii++){
 349       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 350   }
 351
 352   //PrintCuts();
 353   return kTRUE;
 354 }
 355 ///________________________________________________________________________
 356 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
 357   ///Set individual cut ID
 358
 359   //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 360         switch (cutID) {
 361         case kMesonKind:
 362                 if( SetMesonKind(value)) {
 363                 fCuts[kMesonKind] = value;
 364                 UpdateCutString(cutID, value);
 365                 return kTRUE;
 366                 } else return kFALSE;
 367
 368
 369         case kchi2MesonCut:
 370                 if( SetChi2MesonCut(value)) {
 371                 fCuts[kchi2MesonCut] = value;
 372                 UpdateCutString(cutID, value);
 373                 return kTRUE;
 374                 } else return kFALSE;
 375
 376         case kalphaMesonCut:
 377                 if( SetAlphaMesonCut(value)) {
 378                 fCuts[kalphaMesonCut] = value;
 379                 UpdateCutString(cutID, value);
 380                 return kTRUE;
 381                 } else return kFALSE;
 382
 383         case kRCut:
 384                 if( SetRCut(value)) {
 385                 fCuts[kRCut] = value;
 386                 UpdateCutString(cutID, value);
 387                 return kTRUE;
 388                 } else return kFALSE;
 389
 390         case kRapidityMesonCut:
 391                 if( SetRapidityMesonCut(value)) {
 392                 fCuts[kRapidityMesonCut] = value;
 393                 UpdateCutString(cutID, value);
 394                 return kTRUE;
 395                 } else return kFALSE;
 396
 397         case kBackgroundScheme:
 398                 if( SetBackgroundScheme(value)) {
 399                 fCuts[kBackgroundScheme] = value;
 400                 UpdateCutString(cutID, value);
 401                 return kTRUE;
 402                 } else return kFALSE;
 403
 404         case kDegreesForRotationMethod:
 405                 if( SetNDegreesForRotationMethod(value)) {
 406                 fCuts[kDegreesForRotationMethod] = value;
 407                 UpdateCutString(cutID, value);
 408                 return kTRUE;
 409                 } else return kFALSE;
 410
 411         case kNumberOfBGEvents:
 412                 if( SetNumberOfBGEvents(value)) {
 413                 fCuts[kNumberOfBGEvents] = value;
 414                 UpdateCutString(cutID, value);
 415                 return kTRUE;
 416                 } else return kFALSE;
 417
 418         case kuseMCPSmearing:
 419                 if( SetMCPSmearing(value)) {
 420                 fCuts[kuseMCPSmearing] = value;
 421                 UpdateCutString(cutID, value);
 422                 return kTRUE;
 423                 } else return kFALSE;
 424
 425
 426
 427         case kElecShare:
 428                 if( SetSharedElectronCut(value)) {
 429                 fCuts[kElecShare] = value;
 430                 UpdateCutString(cutID, value);
 431                 return kTRUE;
 432                 } else return kFALSE;
 433
 434         case kToCloseV0s:
 435                 if( SetToCloseV0sCut(value)) {
 436                 fCuts[kToCloseV0s] = value;
 437                 UpdateCutString(cutID, value);
 438                 return kTRUE;
 439                 } else return kFALSE;
 440
 441         case kNCuts:
 442                 cout << "Error:: Cut id out of range"<< endl;
 443                 return kFALSE;
 444         }
 445
 446   cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 447   return kFALSE;
 448
 449 }
 450
 451
 452 ///________________________________________________________________________
 453 void AliConversionMesonCuts::PrintCuts() {
 454         // Print out current Cut Selection
 455         for(Int_t ic = 0; ic < kNCuts; ic++) {
 456                 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 457         }
 458 }
 459
 460 ///________________________________________________________________________
 461 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
 462         // Set Cut
 463         switch(mesonKind){
 464                 case 0:
 465                         fMesonKind=0;
 466                         break;
 467                 case 1:
 468                         fMesonKind=1;;
 469                         break;
 470                 default:
 471                         cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
 472                         return kFALSE;
 473         }
 474         return kTRUE;
 475 }
 476
 477 ///________________________________________________________________________
 478 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
 479         // Set Cut
 480         switch(rCut){
 481                 case 0:
 482                         fMaxR = 180.;
 483                         break;
 484                 case 1:
 485                         fMaxR = 180.;
 486                         break;
 487                 case 2:
 488                         fMaxR = 180.;
 489                         break;
 490                 case 3:
 491                         fMaxR = 70.;
 492                         break;
 493                 case 4:
 494                         fMaxR = 70.;
 495                         break;
 496                 case 5:
 497                         fMaxR = 180.;
 498                         break;
 499                         // High purity cuts for PbPb
 500                 case 9:
 501                         fMaxR = 180.;
 502                         break;
 503                 default:
 504                         cout<<"Warning: rCut not defined"<<rCut<<endl;
 505                         return kFALSE;
 506         }
 507         return kTRUE;
 508 }
 509
 510 ///________________________________________________________________________
 511 Bool_t AliConversionMesonCuts::SetChi2MesonCut(Int_t chi2MesonCut){
 512         // Set Cut
 513         switch(chi2MesonCut){
 514                 case 0:  // 100.
 515                 fChi2CutMeson = 100.;
 516                 break;
 517                 case 1:  // 50.
 518                 fChi2CutMeson = 50.;
 519                 break;
 520                 case 2:  // 30.
 521                 fChi2CutMeson = 30.;
 522                 break;
 523                 case 3:
 524                 fChi2CutMeson = 200.;
 525                 break;
 526                 case 4:
 527                 fChi2CutMeson = 500.;
 528                 break;
 529                 case 5:
 530                 fChi2CutMeson = 1000.;
 531                 break;
 532                 default:
 533                 cout<<"Warning: Chi2MesonCut not defined "<<chi2MesonCut<<endl;
 534                 return kFALSE;
 535         }
 536         return kTRUE;
 537 }
 538
 539
 540 ///________________________________________________________________________
 541 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 542 {   // Set Cut
 543         switch(alphaMesonCut){
 544                 case 0: // 0- 0.7
 545                 fAlphaMinCutMeson        = 0.0;
 546                 fAlphaCutMeson   = 0.7;
 547                 break;
 548                         case 1: // 0-0.5
 549                 fAlphaMinCutMeson        = 0.0;
 550                 fAlphaCutMeson   = 0.5;
 551                 break;
 552                         case 2: // 0.5-1
 553                 fAlphaMinCutMeson        = 0.5;
 554                 fAlphaCutMeson   = 1.;
 555                 break;
 556                         case 3: // 0.0-1
 557                 fAlphaMinCutMeson        = 0.0;
 558                 fAlphaCutMeson   = 1.;
 559                 break;
 560                         case 4: // 0-0.65
 561                 fAlphaMinCutMeson        = 0.0;
 562                 fAlphaCutMeson   = 0.65;
 563                 break;
 564                         case 5: // 0-0.75
 565                 fAlphaMinCutMeson        = 0.0;
 566                 fAlphaCutMeson   = 0.75;
 567                 break;
 568                         case 6: // 0-0.8
 569                 fAlphaMinCutMeson        = 0.0;
 570                 fAlphaCutMeson   = 0.8;
 571                 break;
 572                         case 7: // 0.0-0.85
 573                 fAlphaMinCutMeson        = 0.0;
 574                 fAlphaCutMeson   = 0.85;
 575                 break;
 576                         case 8: // 0.0-0.6
 577                 fAlphaMinCutMeson        = 0.0;
 578                 fAlphaCutMeson   = 0.6;
 579                 break;
 580                         default:
 581                                 cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
 582                 return kFALSE;
 583         }
 584         return kTRUE;
 585 }
 586
 587 ///________________________________________________________________________
 588 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
 589         // Set Cut
 590         switch(RapidityMesonCut){
 591                 case 0:  //
 592                 fRapidityCutMeson   = 0.9;
 593                 break;
 594                 case 1:  //
 595                 fRapidityCutMeson   = 0.8;
 596                 break;
 597                 case 2:  //
 598                 fRapidityCutMeson   = 0.7;
 599                 break;
 600                 case 3:  //
 601                 fRapidityCutMeson   = 0.6;
 602                 break;
 603                 case 4:  //
 604                 fRapidityCutMeson   = 0.5;
 605                 break;
 606
 607                 default:
 608                         cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
 609                 return kFALSE;
 610         }
 611     return kTRUE;
 612 }
 613
 614
 615 ///________________________________________________________________________
 616 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
 617         // Set Cut
 618         switch(BackgroundScheme){
 619                 case 0: //Rotation
 620                 fUseRotationMethodInBG=kTRUE;
 621                 fdoBGProbability=kFALSE;
 622                 break;
 623                 case 1: // mixed event with V0 multiplicity
 624                 fUseRotationMethodInBG=kFALSE;
 625                 fUseTrackMultiplicityForBG=kFALSE;
 626                 fdoBGProbability=kFALSE;
 627                 break;
 628                 case 2: // mixed event with track multiplicity
 629                 fUseRotationMethodInBG=kFALSE;
 630                 fUseTrackMultiplicityForBG=kTRUE;
 631                 fdoBGProbability=kFALSE;
 632                 break;
 633                 case 3: //Rotation
 634                 fUseRotationMethodInBG=kTRUE;
 635                 fdoBGProbability=kTRUE;
 636                 break;
 637                 default:
 638                 cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
 639                 return kFALSE;
 640         }
 641         return kTRUE;
 642 }
 643
 644
 645 ///________________________________________________________________________
 646 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
 647     // Set Cut
 648         switch(DegreesForRotationMethod){
 649                 case 0:
 650                 fnDegreeRotationPMForBG = 5;
 651                 break;
 652                 case 1:
 653                 fnDegreeRotationPMForBG = 10;
 654                 break;
 655                 case 2:
 656                 fnDegreeRotationPMForBG = 15;
 657                 break;
 658                 case 3:
 659                 fnDegreeRotationPMForBG = 20;
 660                 break;
 661                 default:
 662                         cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
 663         return kFALSE;
 664         }
 665         fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
 666         return kTRUE;
 667 }
 668
 669 ///________________________________________________________________________
 670 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
 671         // Set Cut
 672         switch(NumberOfBGEvents){
 673                 case 0:
 674                 fNumberOfBGEvents = 5;
 675                 break;
 676                 case 1:
 677                 fNumberOfBGEvents = 10;
 678                 break;
 679                 case 2:
 680                 fNumberOfBGEvents = 15;
 681                 break;
 682                 case 3:
 683                 fNumberOfBGEvents = 20;
 684                 break;
 685                 case 4:
 686                 fNumberOfBGEvents = 2;
 687                 break;
 688                 case 5:
 689                 fNumberOfBGEvents = 50;
 690                 break;
 691                 case 6:
 692                 fNumberOfBGEvents = 80;
 693                 break;
 694                 case 7:
 695                 fNumberOfBGEvents = 100;
 696                 break;
 697                 default:
 698                         cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
 699         return kFALSE;
 700         }
 701         return kTRUE;
 702 }
 703
 704 // ///________________________________________________________________________
 705 // Bool_t AliConversionMesonCuts::SetPsiPairCut(Int_t psiCut) {
 706 //
 707 //      switch(psiCut) {
 708 //              case 0:
 709 //                      fPsiPairCut = 10000; //
 710 //                      break;
 711 //              case 1:
 712 //                      fPsiPairCut = 0.1; //
 713 //                      break;
 714 //              case 2:
 715 //                      fPsiPairCut = 0.05; // Standard
 716 //                      break;
 717 //              case 3:
 718 //                      fPsiPairCut = 0.035; //
 719 //                      break;
 720 //              case 4:
 721 //                      fPsiPairCut = 0.15; //
 722 //                      break;
 723 //              case 5:
 724 //                      fPsiPairCut = 0.2; //
 725 //                      break;
 726 //              case 6:
 727 //                      fPsiPairCut = 0.03; //
 728 //                      break;
 729 //              case 7:
 730 //                      fPsiPairCut = 0.025; //
 731 //                      break;
 732 //              case 8:
 733 //                      fPsiPairCut = 0.01; //
 734 //                      break;
 735 //              default:
 736 //                      cout<<"Warning: PsiPairCut not defined "<<fPsiPairCut<<endl;
 737 //                      return kFALSE;
 738 //      }
 739 //
 740 //      return kTRUE;
 741 // }
 742
 743
 744 ///________________________________________________________________________
 745 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
 746
 747         switch(sharedElec){
 748                 case 0:
 749                         fDoSharedElecCut = kFALSE;
 750                 break;
 751                 case 1:
 752                         fDoSharedElecCut = kTRUE;
 753                 break;
 754                 default:
 755                 cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
 756                 return kFALSE;
 757         }
 758
 759         return kTRUE;
 760 }
 761
 762 ///________________________________________________________________________
 763 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
 764
 765         switch(toClose){
 766                 case 0:
 767                         fDoToCloseV0sCut = kFALSE;
 768                         fminV0Dist = 250;
 769                         break;
 770                 case 1:
 771                         fDoToCloseV0sCut = kTRUE;
 772                         fminV0Dist = 1;
 773                         break;
 774                 case 2:
 775                         fDoToCloseV0sCut = kTRUE;
 776                         fminV0Dist = 2;
 777                         break;
 778                 case 3:
 779                         fDoToCloseV0sCut = kTRUE;
 780                         fminV0Dist = 3;
 781                         break;
 782                 default:
 783                 cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
 784                 return kFALSE;
 785         }
 786         return kTRUE;
 787 }
 788
 789 ///________________________________________________________________________
 790 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
 791 {// Set Cut
 792     switch(useMCPSmearing){
 793     case 0:
 794         fUseMCPSmearing=0;
 795         fPBremSmearing=1.;
 796         fPSigSmearing=0.;
 797         fPSigSmearingCte=0.;
 798         break;
 799     case 1:
 800         fUseMCPSmearing=1;
 801         fPBremSmearing=1.0e-14;
 802         fPSigSmearing=0.;
 803         fPSigSmearingCte=0.;
 804         break;
 805     case 2:
 806         fUseMCPSmearing=1;
 807         fPBremSmearing=1.0e-15;
 808         fPSigSmearing=0.0;
 809         fPSigSmearingCte=0.;
 810         break;
 811     case 3:
 812         fUseMCPSmearing=1;
 813         fPBremSmearing=1.;
 814         fPSigSmearing=0.003;
 815         fPSigSmearingCte=0.002;
 816         break;
 817     case 4:
 818         fUseMCPSmearing=1;
 819         fPBremSmearing=1.;
 820         fPSigSmearing=0.003;
 821         fPSigSmearingCte=0.007;
 822         break;
 823     case 5:
 824         fUseMCPSmearing=1;
 825         fPBremSmearing=1.;
 826         fPSigSmearing=0.003;
 827         fPSigSmearingCte=0.016;
 828         break;
 829     case 6:
 830         fUseMCPSmearing=1;
 831         fPBremSmearing=1.;
 832         fPSigSmearing=0.007;
 833         fPSigSmearingCte=0.016;
 834         break;
 835     case 7:
 836         fUseMCPSmearing=1;
 837         fPBremSmearing=1.0e-16;
 838         fPSigSmearing=0.0;
 839         fPSigSmearingCte=0.;
 840         break;
 841     case 8:
 842         fUseMCPSmearing=1;
 843         fPBremSmearing=1.;
 844         fPSigSmearing=0.007;
 845         fPSigSmearingCte=0.014;
 846         break;
 847     case 9:
 848         fUseMCPSmearing=1;
 849         fPBremSmearing=1.;
 850         fPSigSmearing=0.007;
 851         fPSigSmearingCte=0.011;
 852         break;
 853
 854     default:
 855         AliError("Warning: UseMCPSmearing not defined");
 856         return kFALSE;
 857     }
 858     return kTRUE;
 859 }
 860
 861 // Bool_t AliConversionMesonCuts::PsiPairCut(const AliConversionPhotonBase * photon) const {
 862 //
 863 //      if(photon->GetPsiPair() > fPsiPairCut){
 864 //      return kFALSE;}
 865 //      else{return kTRUE;}
 866 //
 867 // }
 868
 869 ///________________________________________________________________________
 870 TString AliConversionMesonCuts::GetCutNumber(){
 871         // returns TString with current cut number
 872         TString a(kNCuts);
 873         for(Int_t ii=0;ii<kNCuts;ii++){
 874                 a.Append(Form("%d",fCuts[ii]));
 875         }
 876         return a;
 877 }
 878
 879 ///________________________________________________________________________
 880 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
 881
 882         Int_t posLabel = photon->GetTrackLabelPositive();
 883         Int_t negLabel = photon->GetTrackLabelNegative();
 884
 885         fElectronLabelArray[nV0*2] = posLabel;
 886         fElectronLabelArray[(nV0*2)+1] = negLabel;
 887 }
 888
 889 ///________________________________________________________________________
 890 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
 891
 892         Int_t posLabel = photon->GetTrackLabelPositive();
 893         Int_t negLabel = photon->GetTrackLabelNegative();
 894
 895         for(Int_t i = 0; i<nV0s*2;i++){
 896                 if(i==nV0*2)     continue;
 897                 if(i==(nV0*2)+1) continue;
 898                 if(fElectronLabelArray[i] == posLabel){
 899                         return kFALSE;}
 900                 if(fElectronLabelArray[i] == negLabel){
 901                         return kFALSE;}
 902         }
 903
 904         return kTRUE;
 905 }
 906 ///________________________________________________________________________
 907 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
 908         Double_t posX = photon->GetConversionX();
 909         Double_t posY = photon->GetConversionY();
 910         Double_t posZ = photon->GetConversionZ();
 911
 912         for(Int_t i = 0;i<photons->GetEntries();i++){
 913                 if(nV0 == i) continue;
 914                 AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
 915                 Double_t posCompX = photonComp->GetConversionX();
 916                 Double_t posCompY = photonComp->GetConversionY();
 917                 Double_t posCompZ = photonComp->GetConversionZ();
 918
 919                 Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
 920
 921                 if(dist < fminV0Dist*fminV0Dist){
 922                         if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
 923                         else {
 924                                 return kFALSE;}
 925                 }
 926
 927         }
 928         return kTRUE;
 929 }
 930
 931 ///________________________________________________________________________
 932 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
 933 {
 934    Double_t facPBrem = 1.;
 935    Double_t facPSig = 0.;
 936
 937    Double_t phi=0.;
 938    Double_t theta=0.;
 939    Double_t P=0.;
 940
 941
 942    P=photon->P();
 943    phi=photon->Phi();
 944    if( photon->P()!=0){
 945       theta=acos( photon->Pz()/ photon->P());
 946    }
 947
 948    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
 949       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
 950    }
 951
 952    if( fPBremSmearing != 1.){
 953       if(fBrem!=NULL){
 954          facPBrem = fBrem->GetRandom();
 955       }
 956    }
 957
 958    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
 959    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
 960    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
 961    photon->SetE(photon->P());
 962 }