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