]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/GammaConv/AliConversionMesonCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- changed order of filling validated histograms
[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
485 for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
486
487 TParticle* temp = (TParticle*)fMCStack->Particle( index );
488 switch( temp->GetPdgCode() ) {
489 case 211:
490 posPion = temp;
491 labelPosPion = index;
492 break;
493 case -211:
494 negPion = temp;
495 labelNegPion = index;
496 break;
497 case 111:
498 neutPion = temp;
499 labelNeutPion = index;
500 break;
501 }
502 }
503
504 if( posPion && negPion && neutPion ) return kTRUE;
505 return kFALSE;
506}
507
508
0f8c33c1 509//________________________________________________________________________
510Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
511 // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
512 // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
513
514 if(!fMCStack)return kFALSE;
4803eb1f 515 // if(fMCMother->GetPdgCode()==20443 ){
516 // return kFALSE;
517 // }
0f8c33c1 518 if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
4803eb1f 519 if(fMCMother->R()>fMaxR) return kFALSE; // cuts on distance from collision point
520
521 Double_t rapidity = 10.;
522 if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
523 rapidity=8.-fRapidityShift;
524 }
525 else{
526 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
527 }
528
529 // Rapidity Cut
530 if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
531
532 // Select only -> ChiC radiative (JPsi+gamma) decay channel
0f8c33c1 533 if(fMCMother->GetNDaughters()!=2)return kFALSE;
534
4803eb1f 535 TParticle *jpsi = 0x0;
536 TParticle *gamma = 0x0;
537 TParticle *positron = 0x0;
538 TParticle *electron = 0x0;
539
540 Int_t labeljpsiChiC = -1;
541
542 for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
543
544 TParticle* temp = (TParticle*)fMCStack->Particle( index );
545
546 switch( temp->GetPdgCode() ) {
547 case 443:
548 jpsi = temp;
549 labeljpsiChiC = index;
550 break;
551 case 22:
552 gamma = temp;
553 labelgammaChiC = index;
554 break;
555 }
556 }
0f8c33c1 557
558 if ( !jpsi || ! gamma) return kFALSE;
4803eb1f 559 if(jpsi->GetNDaughters()!=2)return kFALSE;
560
561
562 for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
563 TParticle* temp = (TParticle*)fMCStack->Particle( index );
564 switch( temp->GetPdgCode() ) {
565 case -11:
566 electron = temp;
567 labelelectronChiC = index;
568 break;
569 case 11:
570 positron = temp;
571 labelpositronChiC = index;
572 break;
573 }
574 }
575 if( !electron || !positron) return kFALSE;
576 if( positron && electron && gamma) return kTRUE;
0f8c33c1 577 }
578 return kFALSE;
579}
580
2bb2434e 581///________________________________________________________________________
11c1e680 582Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
2bb2434e 583{
4803eb1f 584
e5b6e8a6 585 // Selection of reconstructed Meson candidates
586 // Use flag IsSignal in order to fill Fill different
587 // histograms for Signal and Background
588 TH1 *hist=0x0;
589
590 if(IsSignal){hist=hMesonCuts;}
591 else{hist=hMesonBGCuts;}
592
593 Int_t cutIndex=0;
594
595 if(hist)hist->Fill(cutIndex);
596 cutIndex++;
597
598 // Undefined Rapidity -> Floating Point exception
599 if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
600 if(hist)hist->Fill(cutIndex);
601 cutIndex++;
d9d6352b 602// cout << "undefined rapidity" << endl;
e5b6e8a6 603 return kFALSE;
604 }
605 else{
606 // PseudoRapidity Cut --> But we cut on Rapidity !!!
607 cutIndex++;
11c1e680 608 if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
e5b6e8a6 609 if(hist)hist->Fill(cutIndex);
d9d6352b 610// cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
e5b6e8a6 611 return kFALSE;
612 }
613 }
614 cutIndex++;
615
616 // Opening Angle Cut
617 //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
618 if( pi0->GetOpeningAngle() < fOpeningAngle){
d9d6352b 619// cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl;
e5b6e8a6 620 if(hist)hist->Fill(cutIndex);
621 return kFALSE;
622 }
623 cutIndex++;
624
976b1f89 625 if ( fAlphaPtDepCut == kTRUE ) {
626
627 fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
628 }
629
630
e5b6e8a6 631 // Alpha Max Cut
632 if(pi0->GetAlpha()>fAlphaCutMeson){
d9d6352b 633// cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl;
e5b6e8a6 634 if(hist)hist->Fill(cutIndex);
635 return kFALSE;
636 }
637 cutIndex++;
638
639 // Alpha Min Cut
640 if(pi0->GetAlpha()<fAlphaMinCutMeson){
d9d6352b 641// cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl;
e5b6e8a6 642 if(hist)hist->Fill(cutIndex);
643 return kFALSE;
644 }
645 cutIndex++;
646
4803eb1f 647 if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
648 if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
649
650 if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
d9d6352b 651// cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl;
4803eb1f 652 if(hist)hist->Fill(cutIndex);
653 return kFALSE;
654 }
655 cutIndex++;
656
657 if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
d9d6352b 658// cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl;
4803eb1f 659 if(hist)hist->Fill(cutIndex);
660 return kFALSE;
661 }
662 cutIndex++;
663
664
665 if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
666
667 if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
d9d6352b 668// cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl;
4803eb1f 669 if(hist)hist->Fill(cutIndex);
670 return kFALSE;
671 }
672 cutIndex++;
673
674
675 if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
676 if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
677 if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
678
e5b6e8a6 679 if(hist)hist->Fill(cutIndex);
680 return kTRUE;
2bb2434e 681}
682
683
684
685///________________________________________________________________________
686///________________________________________________________________________
e5b6e8a6 687Bool_t AliConversionMesonCuts::UpdateCutString() {
688 ///Update the cut string (if it has been created yet)
689
690 if(fCutString && fCutString->GetString().Length() == kNCuts) {
e5b6e8a6 691 fCutString->SetString(GetCutNumber());
692 } else {
e5b6e8a6 693 return kFALSE;
694 }
e5b6e8a6 695 return kTRUE;
2bb2434e 696}
697
698///________________________________________________________________________
699Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
e5b6e8a6 700 // Initialize Cuts from a given Cut string
701 AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
702 if(analysisCutSelection.Length()!=kNCuts) {
703 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
704 return kFALSE;
705 }
706 if(!analysisCutSelection.IsDigit()){
707 AliError("Cut selection contains characters");
708 return kFALSE;
709 }
4803eb1f 710
e5b6e8a6 711 const char *cutSelection = analysisCutSelection.Data();
712#define ASSIGNARRAY(i) fCuts[i] = cutSelection[i] - '0'
713 for(Int_t ii=0;ii<kNCuts;ii++){
714 ASSIGNARRAY(ii);
715 }
716
717 // Set Individual Cuts
718 for(Int_t ii=0;ii<kNCuts;ii++){
2bb2434e 719 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
e5b6e8a6 720 }
2bb2434e 721
e656a821 722 PrintCutsWithValues();
e5b6e8a6 723 return kTRUE;
2bb2434e 724}
725///________________________________________________________________________
726Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
e5b6e8a6 727 ///Set individual cut ID
728
729 //cout << "Updating cut " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
730 switch (cutID) {
731 case kMesonKind:
732 if( SetMesonKind(value)) {
733 fCuts[kMesonKind] = value;
734 UpdateCutString();
735 return kTRUE;
736 } else return kFALSE;
42fe34e2 737 case kSelectionCut:
738 if( SetSelectionWindowCut(value)) {
739 fCuts[kSelectionCut] = value;
e5b6e8a6 740 UpdateCutString();
741 return kTRUE;
742 } else return kFALSE;
743 case kalphaMesonCut:
744 if( SetAlphaMesonCut(value)) {
745 fCuts[kalphaMesonCut] = value;
746 UpdateCutString();
747 return kTRUE;
748 } else return kFALSE;
749 case kRCut:
750 if( SetRCut(value)) {
751 fCuts[kRCut] = value;
752 UpdateCutString();
753 return kTRUE;
754 } else return kFALSE;
755
756 case kRapidityMesonCut:
757 if( SetRapidityMesonCut(value)) {
758 fCuts[kRapidityMesonCut] = value;
759 UpdateCutString();
760 return kTRUE;
761 } else return kFALSE;
762
763 case kBackgroundScheme:
764 if( SetBackgroundScheme(value)) {
765 fCuts[kBackgroundScheme] = value;
766 UpdateCutString();
767 return kTRUE;
768 } else return kFALSE;
769
770 case kDegreesForRotationMethod:
771 if( SetNDegreesForRotationMethod(value)) {
772 fCuts[kDegreesForRotationMethod] = value;
773 UpdateCutString();
774 return kTRUE;
775 } else return kFALSE;
776
777 case kNumberOfBGEvents:
778 if( SetNumberOfBGEvents(value)) {
779 fCuts[kNumberOfBGEvents] = value;
780 UpdateCutString();
781 return kTRUE;
782 } else return kFALSE;
783
784 case kuseMCPSmearing:
785 if( SetMCPSmearing(value)) {
786 fCuts[kuseMCPSmearing] = value;
787 UpdateCutString();
788 return kTRUE;
789 } else return kFALSE;
790 case kElecShare:
791 if( SetSharedElectronCut(value)) {
792 fCuts[kElecShare] = value;
793 UpdateCutString();
794 return kTRUE;
795 } else return kFALSE;
796 case kToCloseV0s:
797 if( SetToCloseV0sCut(value)) {
798 fCuts[kToCloseV0s] = value;
799 UpdateCutString();
800 return kTRUE;
801 } else return kFALSE;
4803eb1f 802 case kDcaGammaGamma:
803 if( SetDCAGammaGammaCut(value)) {
804 fCuts[kDcaGammaGamma] = value;
805 UpdateCutString();
806 return kTRUE;
807 } else return kFALSE;
808 case kDcaZPrimVtx:
809 if( SetDCAZMesonPrimVtxCut(value)) {
810 fCuts[kDcaZPrimVtx] = value;
811 UpdateCutString();
812 return kTRUE;
813 } else return kFALSE;
814 case kDcaRPrimVtx:
815 if( SetDCARMesonPrimVtxCut(value)) {
816 fCuts[kDcaRPrimVtx] = value;
817 UpdateCutString();
818 return kTRUE;
819 } else return kFALSE;
e5b6e8a6 820
821 case kNCuts:
822 cout << "Error:: Cut id out of range"<< endl;
823 return kFALSE;
824 }
4803eb1f 825
e5b6e8a6 826 cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
827 return kFALSE;
4803eb1f 828
2bb2434e 829}
830
831
832///________________________________________________________________________
833void AliConversionMesonCuts::PrintCuts() {
e5b6e8a6 834 // Print out current Cut Selection
835 for(Int_t ic = 0; ic < kNCuts; ic++) {
836 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
837 }
2bb2434e 838}
839
e656a821 840///________________________________________________________________________
841void AliConversionMesonCuts::PrintCutsWithValues() {
842 // Print out current Cut Selection with values
843 printf("\nMeson cutnumber \n");
844 for(Int_t ic = 0; ic < kNCuts; ic++) {
845 printf("%d",fCuts[ic]);
846 }
847 printf("\n\n");
848
849 printf("Meson cuts \n");
850 printf("\t |y| < %3.2f \n", fRapidityCutMeson);
851 printf("\t theta_{open} < %3.2f\n", fOpeningAngle);
852 if (!fAlphaPtDepCut) printf("\t %3.2f < alpha < %3.2f\n", fAlphaMinCutMeson, fAlphaCutMeson);
853 printf("\t dca_{gamma,gamma} > %3.2f\n", fDCAGammaGammaCut);
854 printf("\t dca_{R, prim Vtx} > %3.2f\n", fDCARMesonPrimVtxCut);
855 printf("\t dca_{Z, prim Vtx} > %3.2f\n\n", fDCAZMesonPrimVtxCut);
42fe34e2 856 printf("\t Meson selection window for further analysis %3.3f > M_{gamma,gamma} > %3.3f\n\n", fSelectionLow, fSelectionHigh);
e656a821 857
858 printf("Meson BG settings \n");
859 if (!fDoBG){
860 if (!fUseRotationMethodInBG & !fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing V0 mult \n");
861 if (!fUseRotationMethodInBG & fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing track mult \n");
862 if (fUseRotationMethodInBG )printf("\t BG scheme: rotation \n");
863 if (fdoBGProbability) printf("\t -> use BG probability \n");
864 if (fBackgroundHandler) printf("\t -> use new BG handler \n");
865 printf("\t depth of pool: %d\n", fNumberOfBGEvents);
866 if (fUseRotationMethodInBG )printf("\t degree's for BG rotation: %d\n", fnDegreeRotationPMForBG);
867 }
868
869}
870
871
2bb2434e 872///________________________________________________________________________
873Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
e5b6e8a6 874 // Set Cut
875 switch(mesonKind){
876 case 0:
877 fMesonKind=0;
878 break;
879 case 1:
880 fMesonKind=1;;
881 break;
882 default:
883 cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
884 return kFALSE;
885 }
886 return kTRUE;
2bb2434e 887}
888
889///________________________________________________________________________
890Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
e5b6e8a6 891 // Set Cut
892 switch(rCut){
893 case 0:
894 fMaxR = 180.;
895 break;
896 case 1:
897 fMaxR = 180.;
898 break;
899 case 2:
900 fMaxR = 180.;
901 break;
902 case 3:
4803eb1f 903 fMaxR = 70.;
e5b6e8a6 904 break;
905 case 4:
906 fMaxR = 70.;
907 break;
908 case 5:
909 fMaxR = 180.;
910 break;
911 // High purity cuts for PbPb
912 case 9:
913 fMaxR = 180.;
914 break;
915 default:
916 cout<<"Warning: rCut not defined"<<rCut<<endl;
917 return kFALSE;
918 }
919 return kTRUE;
2bb2434e 920}
921
922///________________________________________________________________________
42fe34e2 923Bool_t AliConversionMesonCuts::SetSelectionWindowCut(Int_t selectionCut){
e5b6e8a6 924 // Set Cut
42fe34e2 925 switch(selectionCut){
926 case 0:
927 fSelectionLow = 0.08;
928 fSelectionHigh = 0.145;
e5b6e8a6 929 break;
42fe34e2 930 case 1:
931 fSelectionLow = 0.1;
932 fSelectionHigh = 0.145;
e5b6e8a6 933 break;
42fe34e2 934 case 2:
935 fSelectionLow = 0.11;
936 fSelectionHigh = 0.145;
937 break;
e5b6e8a6 938 case 3:
42fe34e2 939 fSelectionLow = 0.12;
940 fSelectionHigh = 0.145;
e5b6e8a6 941 break;
942 case 4:
42fe34e2 943 fSelectionLow = 0.1;
944 fSelectionHigh = 0.15;
e5b6e8a6 945 break;
946 case 5:
42fe34e2 947 fSelectionLow = 0.11;
948 fSelectionHigh = 0.15;
949 break;
950 case 6:
951 fSelectionLow = 0.12;
952 fSelectionHigh = 0.15;
e5b6e8a6 953 break;
42fe34e2 954
e5b6e8a6 955 default:
42fe34e2 956 cout<<"Warning: SelectionCut not defined "<<selectionCut<<endl;
e5b6e8a6 957 return kFALSE;
958 }
959 return kTRUE;
2bb2434e 960}
961
962
963///________________________________________________________________________
964Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
965{ // Set Cut
e5b6e8a6 966 switch(alphaMesonCut){
967 case 0: // 0- 0.7
968 fAlphaMinCutMeson = 0.0;
969 fAlphaCutMeson = 0.7;
976b1f89 970 fAlphaPtDepCut = kFALSE;
e5b6e8a6 971 break;
976b1f89 972 case 1: // Updated 31 October 2013 before 0.0 - 0.5
973 if( fFAlphaCut ) delete fFAlphaCut;
974 fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
975 fFAlphaCut->SetParameter(0,0.7);
976 fFAlphaCut->SetParameter(1,1.2);
977 fAlphaMinCutMeson = 0.0;
978 fAlphaCutMeson = -1.0;
979 fAlphaPtDepCut = kTRUE;
e5b6e8a6 980 break;
976b1f89 981 case 2: // Updated 31 October 2013 before 0.5-1
982 if( fFAlphaCut ) delete fFAlphaCut;
983 fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
984 fFAlphaCut->SetParameter(0,0.8);
985 fFAlphaCut->SetParameter(1,1.2);
986 fAlphaMinCutMeson = 0.0;
987 fAlphaCutMeson = -1.0;
988 fAlphaPtDepCut = kTRUE;
e5b6e8a6 989 break;
990 case 3: // 0.0-1
991 fAlphaMinCutMeson = 0.0;
992 fAlphaCutMeson = 1.;
976b1f89 993 fAlphaPtDepCut = kFALSE;
e5b6e8a6 994 break;
995 case 4: // 0-0.65
996 fAlphaMinCutMeson = 0.0;
997 fAlphaCutMeson = 0.65;
976b1f89 998 fAlphaPtDepCut = kFALSE;
e5b6e8a6 999 break;
1000 case 5: // 0-0.75
1001 fAlphaMinCutMeson = 0.0;
1002 fAlphaCutMeson = 0.75;
976b1f89 1003 fAlphaPtDepCut = kFALSE;
e5b6e8a6 1004 break;
1005 case 6: // 0-0.8
1006 fAlphaMinCutMeson = 0.0;
1007 fAlphaCutMeson = 0.8;
976b1f89 1008 fAlphaPtDepCut = kFALSE;
e5b6e8a6 1009 break;
1010 case 7: // 0.0-0.85
1011 fAlphaMinCutMeson = 0.0;
1012 fAlphaCutMeson = 0.85;
976b1f89 1013 fAlphaPtDepCut = kFALSE;
e5b6e8a6 1014 break;
1015 case 8: // 0.0-0.6
1016 fAlphaMinCutMeson = 0.0;
1017 fAlphaCutMeson = 0.6;
976b1f89 1018 fAlphaPtDepCut = kFALSE;
e5b6e8a6 1019 break;
49af2ef8 1020 case 9: // Updated 11 November 2013 before 0.0 - 0.3
1021 if( fFAlphaCut ) delete fFAlphaCut;
1022 fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
1023 fFAlphaCut->SetParameter(0,0.65);
1024 fFAlphaCut->SetParameter(1,1.2);
1025 fAlphaMinCutMeson = 0.0;
1026 fAlphaCutMeson = -1.0;
1027 fAlphaPtDepCut = kTRUE;
e5b6e8a6 1028 break;
1029 default:
1030 cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
1031 return kFALSE;
1032 }
1033 return kTRUE;
2bb2434e 1034}
1035
1036///________________________________________________________________________
4803eb1f 1037Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
e5b6e8a6 1038 // Set Cut
1039 switch(RapidityMesonCut){
621d4c94 1040 case 0: // changed from 0.9 to 1.35
1041 fRapidityCutMeson = 1.35;
e5b6e8a6 1042 break;
1043 case 1: //
1044 fRapidityCutMeson = 0.8;
1045 break;
1046 case 2: //
1047 fRapidityCutMeson = 0.7;
1048 break;
1049 case 3: //
1050 fRapidityCutMeson = 0.6;
1051 break;
1052 case 4: //
1053 fRapidityCutMeson = 0.5;
1054 break;
1055 case 5: //
1056 fRapidityCutMeson = 0.85;
1057 break;
1058 case 6: //
1059 fRapidityCutMeson = 0.75;
1060 break;
11c1e680 1061 case 7: //
1062 fRapidityCutMeson = 0.3;
1063 break;
3c56f740 1064 case 8: //changed, before 0.35
1065 fRapidityCutMeson = 0.25;
11c1e680 1066 break;
1067 case 9: //
1068 fRapidityCutMeson = 0.4;
1069 break;
e5b6e8a6 1070 default:
1071 cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
1072 return kFALSE;
1073 }
1074 return kTRUE;
2bb2434e 1075}
1076
1077
1078///________________________________________________________________________
1079Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
e5b6e8a6 1080 // Set Cut
1081 switch(BackgroundScheme){
1082 case 0: //Rotation
1083 fUseRotationMethodInBG=kTRUE;
1084 fdoBGProbability=kFALSE;
1085 break;
1086 case 1: // mixed event with V0 multiplicity
1087 fUseRotationMethodInBG=kFALSE;
1088 fUseTrackMultiplicityForBG=kFALSE;
1089 fdoBGProbability=kFALSE;
1090 break;
1091 case 2: // mixed event with track multiplicity
1092 fUseRotationMethodInBG=kFALSE;
1093 fUseTrackMultiplicityForBG=kTRUE;
1094 fdoBGProbability=kFALSE;
1095 break;
1096 case 3: //Rotation
1097 fUseRotationMethodInBG=kTRUE;
1098 fdoBGProbability=kTRUE;
1099 break;
1100 case 4: //No BG calculation
1101 cout << "no BG calculation should be done" << endl;
1102 fUseRotationMethodInBG=kFALSE;
1103 fdoBGProbability=kFALSE;
1104 fDoBG=kFALSE;
1105 break;
1106 case 5: //Rotation
1107 fUseRotationMethodInBG=kTRUE;
1108 fdoBGProbability=kFALSE;
1109 fBackgroundHandler = 1;
1110 break;
1111 case 6: // mixed event with V0 multiplicity
1112 fUseRotationMethodInBG=kFALSE;
1113 fUseTrackMultiplicityForBG=kFALSE;
1114 fdoBGProbability=kFALSE;
1115 fBackgroundHandler = 1;
1116 break;
1117 case 7: // mixed event with track multiplicity
1118 fUseRotationMethodInBG=kFALSE;
1119 fUseTrackMultiplicityForBG=kTRUE;
1120 fdoBGProbability=kFALSE;
1121 fBackgroundHandler = 1;
1122 break;
1123 case 8: //Rotation
1124 fUseRotationMethodInBG=kTRUE;
1125 fdoBGProbability=kTRUE;
1126 fBackgroundHandler = 1;
1127 break;
1128 default:
1129 cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1130 return kFALSE;
1131 }
1132 return kTRUE;
2bb2434e 1133}
1134
ca91a3e1 1135
2bb2434e 1136///________________________________________________________________________
1137Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
e5b6e8a6 1138 // Set Cut
1139 switch(DegreesForRotationMethod){
1140 case 0:
1141 fnDegreeRotationPMForBG = 5;
1142 break;
1143 case 1:
1144 fnDegreeRotationPMForBG = 10;
1145 break;
1146 case 2:
1147 fnDegreeRotationPMForBG = 15;
1148 break;
1149 case 3:
1150 fnDegreeRotationPMForBG = 20;
1151 break;
1152 default:
1153 cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1154 return kFALSE;
1155 }
1156 fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1157 return kTRUE;
2bb2434e 1158}
1159
1160///________________________________________________________________________
ca91a3e1 1161Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
e5b6e8a6 1162 // Set Cut
1163 switch(NumberOfBGEvents){
1164 case 0:
1165 fNumberOfBGEvents = 5;
1166 break;
1167 case 1:
1168 fNumberOfBGEvents = 10;
1169 break;
1170 case 2:
1171 fNumberOfBGEvents = 15;
1172 break;
1173 case 3:
1174 fNumberOfBGEvents = 20;
1175 break;
1176 case 4:
1177 fNumberOfBGEvents = 2;
1178 break;
1179 case 5:
1180 fNumberOfBGEvents = 50;
1181 break;
1182 case 6:
1183 fNumberOfBGEvents = 80;
1184 break;
1185 case 7:
1186 fNumberOfBGEvents = 100;
1187 break;
1188 default:
1189 cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1190 return kFALSE;
1191 }
1192 return kTRUE;
2bb2434e 1193}
2bb2434e 1194///________________________________________________________________________
1195Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1196
e5b6e8a6 1197 switch(sharedElec){
1198 case 0:
1199 fDoSharedElecCut = kFALSE;
1200 break;
1201 case 1:
1202 fDoSharedElecCut = kTRUE;
1203 break;
1204 default:
1205 cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1206 return kFALSE;
1207 }
4803eb1f 1208
e5b6e8a6 1209 return kTRUE;
2bb2434e 1210}
1211
1212///________________________________________________________________________
1213Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1214
e5b6e8a6 1215 switch(toClose){
1216 case 0:
1217 fDoToCloseV0sCut = kFALSE;
1218 fminV0Dist = 250;
1219 break;
1220 case 1:
1221 fDoToCloseV0sCut = kTRUE;
1222 fminV0Dist = 1;
1223 break;
1224 case 2:
1225 fDoToCloseV0sCut = kTRUE;
1226 fminV0Dist = 2;
1227 break;
1228 case 3:
1229 fDoToCloseV0sCut = kTRUE;
1230 fminV0Dist = 3;
1231 break;
1232 default:
1233 cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1234 return kFALSE;
1235 }
1236 return kTRUE;
2bb2434e 1237}
1238
ca91a3e1 1239///________________________________________________________________________
1240Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1241{// Set Cut
e5b6e8a6 1242 switch(useMCPSmearing){
1243 case 0:
1244 fUseMCPSmearing=0;
1245 fPBremSmearing=1.;
1246 fPSigSmearing=0.;
1247 fPSigSmearingCte=0.;
1248 break;
1249 case 1:
1250 fUseMCPSmearing=1;
1251 fPBremSmearing=1.0e-14;
1252 fPSigSmearing=0.;
1253 fPSigSmearingCte=0.;
1254 break;
1255 case 2:
1256 fUseMCPSmearing=1;
1257 fPBremSmearing=1.0e-15;
1258 fPSigSmearing=0.0;
1259 fPSigSmearingCte=0.;
1260 break;
1261 case 3:
1262 fUseMCPSmearing=1;
1263 fPBremSmearing=1.;
1264 fPSigSmearing=0.003;
1265 fPSigSmearingCte=0.002;
1266 break;
1267 case 4:
1268 fUseMCPSmearing=1;
1269 fPBremSmearing=1.;
1270 fPSigSmearing=0.003;
1271 fPSigSmearingCte=0.007;
1272 break;
1273 case 5:
1274 fUseMCPSmearing=1;
1275 fPBremSmearing=1.;
1276 fPSigSmearing=0.003;
1277 fPSigSmearingCte=0.016;
1278 break;
1279 case 6:
1280 fUseMCPSmearing=1;
1281 fPBremSmearing=1.;
1282 fPSigSmearing=0.007;
1283 fPSigSmearingCte=0.016;
1284 break;
1285 case 7:
1286 fUseMCPSmearing=1;
1287 fPBremSmearing=1.0e-16;
1288 fPSigSmearing=0.0;
1289 fPSigSmearingCte=0.;
1290 break;
1291 case 8:
1292 fUseMCPSmearing=1;
1293 fPBremSmearing=1.;
1294 fPSigSmearing=0.007;
1295 fPSigSmearingCte=0.014;
1296 break;
1297 case 9:
1298 fUseMCPSmearing=1;
1299 fPBremSmearing=1.;
1300 fPSigSmearing=0.007;
1301 fPSigSmearingCte=0.011;
1302 break;
1303
1304 default:
1305 AliError("Warning: UseMCPSmearing not defined");
1306 return kFALSE;
1307 }
1308 return kTRUE;
ca91a3e1 1309}
4803eb1f 1310
1311
1312///________________________________________________________________________
1313Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1314 // Set Cut
1315 switch(DCAGammaGamma){
1316 case 0: //
1317 fDCAGammaGammaCut = 1000;
1318 break;
1319 case 1: //
1320 fDCAGammaGammaCut = 10;
1321 break;
1322 case 2: //
1323 fDCAGammaGammaCut = 5;
1324 break;
1325 case 3: //
1326 fDCAGammaGammaCut = 4;
1327 break;
1328 case 4: //
1329 fDCAGammaGammaCut = 3;
1330 break;
1331 case 5: //
1332 fDCAGammaGammaCut = 2.5;
1333 break;
1334 case 6: //
1335 fDCAGammaGammaCut = 2;
1336 break;
1337 case 7: //
1338 fDCAGammaGammaCut = 1.5;
1339 break;
1340 case 8: //
1341 fDCAGammaGammaCut = 1;
1342 break;
1343 case 9: //
1344 fDCAGammaGammaCut = 0.5;
1345 break;
1346 default:
1347 cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1348 return kFALSE;
1349 }
1350 return kTRUE;
1351}
1352
1353///________________________________________________________________________
1354Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1355 // Set Cut
1356 switch(DCAZMesonPrimVtx){
1357 case 0: //
1358 fDCAZMesonPrimVtxCut = 1000;
1359 break;
1360 case 1: //
1361 fDCAZMesonPrimVtxCut = 10;
1362 break;
1363 case 2: //
1364 fDCAZMesonPrimVtxCut = 5;
1365 break;
1366 case 3: //
1367 fDCAZMesonPrimVtxCut = 4;
1368 break;
1369 case 4: //
1370 fDCAZMesonPrimVtxCut = 3;
1371 break;
1372 case 5: //
1373 fDCAZMesonPrimVtxCut = 2.5;
1374 break;
1375 case 6: //
1376 fDCAZMesonPrimVtxCut = 2;
1377 break;
1378 case 7: //
1379 fDCAZMesonPrimVtxCut = 1.5;
1380 break;
1381 case 8: //
1382 fDCAZMesonPrimVtxCut = 1;
1383 break;
1384 case 9: //
1385 fDCAZMesonPrimVtxCut = 0.5;
1386 break;
1387 default:
1388 cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1389 return kFALSE;
1390 }
1391 return kTRUE;
1392}
1393
1394///________________________________________________________________________
1395Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1396 // Set Cut
1397 switch(DCARMesonPrimVtx){
1398 case 0: //
1399 fDCARMesonPrimVtxCut = 1000;
1400 break;
1401 case 1: //
1402 fDCARMesonPrimVtxCut = 10;
1403 break;
1404 case 2: //
1405 fDCARMesonPrimVtxCut = 5;
1406 break;
1407 case 3: //
1408 fDCARMesonPrimVtxCut = 4;
1409 break;
1410 case 4: //
1411 fDCARMesonPrimVtxCut = 3;
1412 break;
1413 case 5: //
1414 fDCARMesonPrimVtxCut = 2.5;
1415 break;
1416 case 6: //
1417 fDCARMesonPrimVtxCut = 2;
1418 break;
1419 case 7: //
1420 fDCARMesonPrimVtxCut = 1.5;
1421 break;
1422 case 8: //
1423 fDCARMesonPrimVtxCut = 1;
1424 break;
1425 case 9: //
1426 fDCARMesonPrimVtxCut = 0.5;
1427 break;
1428 default:
1429 cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1430 return kFALSE;
1431 }
1432 return kTRUE;
1433}
1434
1435
2bb2434e 1436///________________________________________________________________________
1437TString AliConversionMesonCuts::GetCutNumber(){
e5b6e8a6 1438 // returns TString with current cut number
1439 TString a(kNCuts);
1440 for(Int_t ii=0;ii<kNCuts;ii++){
1441 a.Append(Form("%d",fCuts[ii]));
1442 }
1443 return a;
2bb2434e 1444}
1445
1446///________________________________________________________________________
1447void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
4803eb1f 1448
e5b6e8a6 1449 Int_t posLabel = photon->GetTrackLabelPositive();
1450 Int_t negLabel = photon->GetTrackLabelNegative();
4803eb1f 1451
e5b6e8a6 1452 fElectronLabelArray[nV0*2] = posLabel;
1453 fElectronLabelArray[(nV0*2)+1] = negLabel;
2bb2434e 1454}
1455
1456///________________________________________________________________________
1457Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1458
e5b6e8a6 1459 Int_t posLabel = photon->GetTrackLabelPositive();
1460 Int_t negLabel = photon->GetTrackLabelNegative();
2bb2434e 1461
e5b6e8a6 1462 for(Int_t i = 0; i<nV0s*2;i++){
1463 if(i==nV0*2) continue;
1464 if(i==(nV0*2)+1) continue;
1465 if(fElectronLabelArray[i] == posLabel){
1466 return kFALSE;}
1467 if(fElectronLabelArray[i] == negLabel){
1468 return kFALSE;}
1469 }
2bb2434e 1470
e5b6e8a6 1471 return kTRUE;
2bb2434e 1472}
1473///________________________________________________________________________
1474Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
e5b6e8a6 1475 Double_t posX = photon->GetConversionX();
1476 Double_t posY = photon->GetConversionY();
1477 Double_t posZ = photon->GetConversionZ();
1478
1479 for(Int_t i = 0;i<photons->GetEntries();i++){
1480 if(nV0 == i) continue;
1481 AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1482 Double_t posCompX = photonComp->GetConversionX();
1483 Double_t posCompY = photonComp->GetConversionY();
1484 Double_t posCompZ = photonComp->GetConversionZ();
4803eb1f 1485
e5b6e8a6 1486 Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
2bb2434e 1487
e5b6e8a6 1488 if(dist < fminV0Dist*fminV0Dist){
1489 if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1490 else {
1491 return kFALSE;}
1492 }
4803eb1f 1493
e5b6e8a6 1494 }
1495 return kTRUE;
2bb2434e 1496}
ca91a3e1 1497
1498///________________________________________________________________________
1499void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1500{
4803eb1f 1501
0a2b2b4b 1502 if (photon==NULL) return;
ca91a3e1 1503 Double_t facPBrem = 1.;
1504 Double_t facPSig = 0.;
1505
1506 Double_t phi=0.;
1507 Double_t theta=0.;
1508 Double_t P=0.;
1509
4803eb1f 1510
ca91a3e1 1511 P=photon->P();
1512 phi=photon->Phi();
1513 if( photon->P()!=0){
1514 theta=acos( photon->Pz()/ photon->P());
1515 }
1516
4803eb1f 1517 if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
ca91a3e1 1518 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1519 }
4803eb1f 1520
ca91a3e1 1521 if( fPBremSmearing != 1.){
1522 if(fBrem!=NULL){
1523 facPBrem = fBrem->GetRandom();
1524 }
1525 }
1526
1527 photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1528 photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1529 photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1530 photon->SetE(photon->P());
1531}
a1c2f90c 1532///________________________________________________________________________
1533void AliConversionMesonCuts::SmearVirtualPhoton(AliAODConversionPhoton* photon)
1534{
1535
1536 if (photon==NULL) return;
1537 Double_t facPBrem = 1.;
1538 Double_t facPSig = 0.;
1539
1540 Double_t phi=0.;
1541 Double_t theta=0.;
1542 Double_t P=0.;
1543
1544
1545 P=photon->P();
1546 phi=photon->Phi();
1547 if( photon->P()!=0){
1548 theta=acos( photon->Pz()/ photon->P());
1549 }
1550
1551 if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1552 facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1553 }
1554
1555 if( fPBremSmearing != 1.){
1556 if(fBrem!=NULL){
1557 facPBrem = fBrem->GetRandom();
1558 }
1559 }
1560
1561 photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1562 photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1563 photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1564
7ab8e97e 1565}
1566///________________________________________________________________________
1567TLorentzVector AliConversionMesonCuts::SmearElectron(TLorentzVector particle)
1568{
1569
1570 //if (particle==0) return;
1571 Double_t facPBrem = 1.;
1572 Double_t facPSig = 0.;
1573
1574 Double_t phi=0.;
1575 Double_t theta=0.;
1576 Double_t P=0.;
1577
1578
1579 P=particle.P();
1580
1581
1582 phi=particle.Phi();
1583 if (phi < 0.) phi += 2. * TMath::Pi();
1584
1585 if( particle.P()!=0){
1586 theta=acos( particle.Pz()/ particle.P());
1587 }
1588
1589
1590 Double_t fPSigSmearingHalf = fPSigSmearing / 2.0; //The parameter was set for gammas with 2 particles and here we have just one electron
1591 Double_t sqrtfPSigSmearingCteHalf = fPSigSmearingCte / 2.0 ; //The parameter was set for gammas with 2 particles and here we have just one electron
1592
1593
1594
1595 if( fPSigSmearingHalf != 0. || sqrtfPSigSmearingCteHalf!=0. ){
1596 facPSig = TMath::Sqrt(sqrtfPSigSmearingCteHalf*sqrtfPSigSmearingCteHalf+fPSigSmearingHalf*fPSigSmearingHalf*P*P)*fRandom.Gaus(0.,1.);
1597 }
1598
1599 if( fPBremSmearing != 1.){
1600 if(fBrem!=NULL){
1601 facPBrem = fBrem->GetRandom();
1602 }
1603 }
1604
1605 TLorentzVector SmearedParticle;
1606
1607 SmearedParticle.SetXYZM( facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) , facPBrem* (1+facPSig)* P*sin(theta)*sin(phi) ,
1608 facPBrem* (1+facPSig)* P*cos(theta) , TDatabasePDG::Instance()->GetParticle( ::kElectron )->Mass()) ;
1609
1610 //particle.SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1611 //particle.SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1612 //particle.SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1613
1614 return SmearedParticle;
a1c2f90c 1615
1616}
1617