1 #include "AliV0ReaderV1.h"
2 #include "AliKFParticle.h"
5 #include "AliAODEvent.h"
6 #include "AliESDEvent.h"
7 #include "AliKFParticle.h"
8 #include "AliKFConversionPhoton.h"
9 #include "AliAODConversionPhoton.h"
10 #include "AliConversionPhotonBase.h"
12 #include "AliKFVertex.h"
13 #include "AliAODTrack.h"
14 #include "AliESDtrack.h"
15 #include "AliAnalysisManager.h"
16 #include "AliInputEventHandler.h"
17 #include "AliAODHandler.h"
18 #include "AliPIDResponse.h"
29 ClassImp(AliV0ReaderV1)
31 //________________________________________________________________________
32 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
33 fConversionGammas(NULL),
38 fCurrentMotherKFCandidate(NULL),
39 fCurrentPositiveKFParticle(NULL),
40 fCurrentNegativeKFParticle(NULL),
41 fCurrentTrackLabels(NULL),
50 fMaxR(180),// 100 meter(outside of ALICE)
51 fMinR(5),// 100 meter(outside of ALICE)
61 fLineCutZRSlopeMin(0.),
62 fLineCutZValueMin(-2),
63 fChi2CutConversion(30),
64 fPIDProbabilityCutNegativeParticle(0),
65 fPIDProbabilityCutPositiveParticle(0),
66 fDodEdxSigmaCut(kTRUE),
67 fDoTOFsigmaCut(kFALSE), // RRnewTOF
68 fPIDTRDEfficiency(0.95),
70 fPIDnSigmaAboveElectronLine(5),
71 fPIDnSigmaBelowElectronLine(-3),
72 fTofPIDnSigmaAboveElectronLine(100), // RRnewTOF
73 fTofPIDnSigmaBelowElectronLine(-100), // RRnewTOF
74 fPIDnSigmaAbovePionLine(0),
75 fPIDnSigmaAbovePionLineHighPt(-100),
76 fPIDMinPnSigmaAbovePionLine(1),
77 fPIDMaxPnSigmaAbovePionLine(3),
78 fDoKaonRejectionLowP(kTRUE),
79 fDoProtonRejectionLowP(kTRUE),
80 fDoPionRejectionLowP(kTRUE),
81 fPIDnSigmaAtLowPAroundKaonLine(0),
82 fPIDnSigmaAtLowPAroundProtonLine(0),
83 fPIDnSigmaAtLowPAroundPionLine(0),
84 fPIDMinPKaonRejectionLowP(1.5),
85 fPIDMinPProtonRejectionLowP(2),
86 fPIDMinPPionRejectionLowP(0.5),
87 fDoQtGammaSelection(kTRUE),
88 fDoHighPtQtGammaSelection(kTRUE), // RRnew
90 fHighPtQtMax(0.06), // RRnew
91 fPtBorderForQt(2.5), // RRnew
93 fUseImprovedVertex(kTRUE),
94 fUseOwnXYZCalculation(kTRUE),
95 fUseConstructGamma(kFALSE),
96 fUseEtaMinCut(kFALSE),
97 fUseOnFlyV0Finder(kTRUE),
98 fDoPhotonAsymmetryCut(kTRUE),
99 fMinPPhotonAsymmetryCut(100.),
100 fMinPhotonAsymmetry(0.),
101 kUseAODConversionPhoton(kFALSE),
104 fDeltaAODFilename("AliAODGammaConversion.root")
107 fLineCutZRSlope=tan(2*atan(exp(-fEtaCut)));
108 // Input slot #0 works with a TChain
109 DefineInput(0, TChain::Class());
110 // Output slot #0 id reserved by the base class for AOD
111 // Output slot #1 writes into a TH1 container
112 DefineOutput(1, TList::Class());
114 fCurrentTrackLabels=new Int_t[2];
118 //________________________________________________________________________
119 AliV0ReaderV1::~AliV0ReaderV1()
121 if(fConversionGammas){
122 fConversionGammas->Delete();// Clear Objects
123 delete fConversionGammas;
124 fConversionGammas=0x0;
127 if(fCurrentTrackLabels){
128 delete[] fCurrentTrackLabels;
129 fCurrentTrackLabels=NULL;}
133 //________________________________________________________________________
134 void AliV0ReaderV1::UserCreateOutputObjects()
137 if(fOutputList != NULL){
141 if(fOutputList == NULL){
142 fOutputList = new TList();
143 fOutputList->SetOwner(kTRUE);
146 TList *fCutList=new TList();
147 fCutList->SetName("GammaReconstruction");
148 fCutList->SetOwner(kTRUE);
149 fOutputList->Add(fCutList);
152 Int_t kGCnXBinsGammaMass = 4000;
153 Double_t kGCfirstXBinGammaMass= 0.;
154 Double_t kGClastXBinGammaMass = 1.;
156 Int_t kGCnYBinsSpectra = 250;
157 Double_t kGCfirstYBinSpectra = 0.;
158 Double_t kGClastYBinSpectra = 25.;
160 // Process Gammas Histograms
162 hV0CurrentFinder=new TH1F("ESD_V0sCurrentFinder_InvMass","V0sCurrentFinder",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
163 fCutList->Add(hV0CurrentFinder);
164 hV0AllArmenteros=new TH2F("ESD_V0sCurrentFinder_Armenteros","Armenteros Alpha Qt",200,-1,1,250,0,0.25);
165 fCutList->Add(hV0AllArmenteros);
166 hV0Good=new TH1F("ESD_GoodV0s_InvMass","GoodV0s",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
167 fCutList->Add(hV0Good);
168 hV0GoodArmenteros=new TH2F("ESD_GoodV0s_Armenteros","Armenteros Alpha Qt",200,-1,1,250,0,0.25);
169 fCutList->Add(hV0GoodArmenteros);
173 hV0CutLikeSign=new TH1F("ESD_CutLikeSign_InvMass","LikeSign",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
174 fCutList->Add(hV0CutLikeSign);
175 hV0CutRefit=new TH1F("ESD_CutRefit_InvMass","No TPC refit",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
176 fCutList->Add(hV0CutRefit);
177 hV0CutKinks=new TH1F("ESD_CutKink_InvMass","Kinks",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
178 fCutList->Add(hV0CutKinks);
179 hV0CutMinNclsTPCToF=new TH1F("ESD_CutMinNClsTPCToF_InvMass","Min Ncls TPC ToF",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
180 fCutList->Add(hV0CutMinNclsTPCToF);
183 /*hV0CutNContributors=new TH1F("ESD_CutNContributors_InvMass","NContributors<=0",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
184 fCutList->Add(hV0CutNContributors);
185 hV0CutVertexZ=new TH1F("ESD_CutVertexZ_InvMass","VertexZ",kGCnXBinsGammaMass,kGCfirstXBinGammaMass,kGClastXBinGammaMass);
186 fCutList->Add(hV0CutVertexZ);
192 hV0CutdEdxElectron=new TH1F("ESD_CutdEdxSigmaElectronLine_InvMass" ,"dedx ElectronLine" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
193 fCutList->Add(hV0CutdEdxElectron);
194 hV0CutdEdxPion=new TH1F("ESD_CutdEdxSigmaPionLine_InvMass" ,"dedx PionLine" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
195 fCutList->Add(hV0CutdEdxPion);
196 hV0CutdEdxKaonLowP=new TH1F("ESD_CutKaonRejectionLowP_InvMass" ,"dedx KaonRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
197 fCutList->Add(hV0CutdEdxKaonLowP);
198 hV0CutdEdxProtonLowP=new TH1F("ESD_CutProtonRejectionLowP_InvMass" ,"dedx ProtonRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
199 fCutList->Add(hV0CutdEdxProtonLowP);
200 hV0CutdEdxPionLowP=new TH1F("ESD_CutPionRejectionLowP_InvMass" ,"dedx PionRejection LowP" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
201 fCutList->Add(hV0CutdEdxPionLowP);
202 hV0CutdEdxTOFElectron=new TH1F("ESD_CutTOFsigmaElec_InvMass", "ESD_CutTOFsigmaElec_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
203 fCutList->Add(hV0CutdEdxTOFElectron);
204 hV0CutdEdxTRD=new TH1F("ESD_CutTRD_InvMass", "ESD_CutTRD_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
205 fCutList->Add(hV0CutdEdxTRD);
206 hGammadEdxbefore=new TH2F("Gamma_dEdx_before","dEdx Gamma before" ,kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,400, 0,200);
207 fCutList->Add(hGammadEdxbefore);
208 hGammadEdxafter=new TH2F("Gamma_dEdx_after","dEdx Gamma after" ,kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,400, 0,200);
209 fCutList->Add(hGammadEdxafter);
214 hV0CutQt=new TH1F("ESD_CutQt_InvMass","ESD_CutQt_InvMass",kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
215 fCutList->Add(hV0CutQt);
219 hV0CutR=new TH1F("ESD_CutR_InvMass" ,"Above RMax" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
220 fCutList->Add(hV0CutR);
221 hV0CutMinR=new TH1F("ESD_CutMinR_InvMass" ,"Above RMax" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
222 fCutList->Add(hV0CutMinR);
223 hV0CutLine=new TH1F("ESD_CutLine_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
224 fCutList->Add(hV0CutLine);
225 hV0CutZ=new TH1F("ESD_CutZ_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
226 fCutList->Add(hV0CutZ);
227 hV0CutEta=new TH1F("ESD_CutEta_InvMass" ,"Above #eta max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
228 fCutList->Add(hV0CutEta);
231 hV0CutSinglePt=new TH1F("ESD_CutSinglePt_InvMass" ,"Below p_{t} min" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
232 fCutList->Add(hV0CutSinglePt);
233 hV0CutNDF=new TH1F("ESD_CutNDF_InvMass" ,"#chi^{2} > Max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
234 fCutList->Add(hV0CutNDF);
235 hV0CutChi2=new TH1F("ESD_CutChi2_InvMass" ,"#chi^{2} > Max" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
236 fCutList->Add(hV0CutChi2);
237 hV0CutPt= new TH1F("ESD_CutPt_InvMass" ,"Below p_{t} min" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
238 fCutList->Add(hV0CutPt);
242 hV0CutAsymmetry=new TH1F("ESD_CutPhotonAsymmetry_InvMass" ,"Out of reconstruction area" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
243 fCutList->Add(hV0CutAsymmetry);
247 hV0CutPIDProb=new TH1F("ESD_CutPIDProb_InvMass" ,"wrong TPC PID" , kGCnXBinsGammaMass, kGCfirstXBinGammaMass, kGClastXBinGammaMass);
248 fCutList->Add(hV0CutPIDProb);
253 TList *fOtherList=new TList();
254 fOtherList->SetName("EventInfo");
255 fOtherList->SetOwner(kTRUE);
256 fOutputList->Add(fOtherList);
259 hV0EventCuts=new TH1F("ESD_EventCuts","Event Cuts",10,-0.5,9.5);
260 fOtherList->Add(hV0EventCuts);
261 hNEvents=new TH1F("NEvents_vs_Centrality","NEvents vs Centrality",fNCentralityBins,-0.5,fNCentralityBins-0.5);
262 fOtherList->Add(hNEvents);
263 hCentrality=new TH1F("Centrality","Centrality",100,0,100);
264 fOtherList->Add(hCentrality);
265 hVertexZ=new TH1F("VertexZ","VertexZ",1000,-50,50);
266 fOtherList->Add(hVertexZ);
270 TList *fGammaList=new TList();
271 fGammaList->SetName("GammaInfo");
272 fGammaList->SetOwner(kTRUE);
273 fOutputList->Add(fGammaList);
275 hGammaPt=new TH1F*[fNCentralityBins];
277 for(int i=0;i<fNCentralityBins;i++){
278 hGammaPt[i]=new TH1F(Form("GammaSpectrum_RECO_Pt_%d",i),"Reco Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
279 fGammaList->Add(hGammaPt[i]);
282 hGammaPhi=new TH1F("Phi_Gamma","Phi Gamma" ,36, 0, 2*TMath::Pi());
283 fGammaList->Add(hGammaPhi);
284 hGammaConversionMapXY=new TH2F("Gamma_ConversionMap_XY","Conversion Point xy",400,-200,200,400,-200,200);
285 fGammaList->Add(hGammaConversionMapXY);
286 hGammaConversionMapZR=new TH2F("Gamma_ConversionMap_ZR","Conversion Point zr",500,-250,250,180,0,180);
287 fGammaList->Add(hGammaConversionMapZR);
290 // FILL MC PART only if MC is available
291 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
293 hMCPtTRUE=new TH1F*[fNCentralityBins];
294 hMCPtRECOTRUE=new TH1F*[fNCentralityBins];
296 for(int i=0;i<fNCentralityBins;i++){
297 hMCPtTRUE[i]=new TH1F(Form("GammaSpectrum_MC_Pt_%d",i),"TRUE Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
298 fGammaList->Add(hMCPtTRUE[i]);
299 hMCPtRECOTRUE[i]=new TH1F(Form("GammaSpectrum_RECOTRUE_Pt_%d",i),"True Reco Gamma Pt",kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
300 fGammaList->Add(hMCPtRECOTRUE[i]);
304 for (Int_t i=0; i<fGammaList->GetEntries(); i++) {
305 TH1 *h1 = dynamic_cast<TH1*>(fGammaList->At(i));
306 if (h1){h1->Sumw2();}
309 hMCPtResolution=new TH2F("Resolution_Gamma_dPt_Pt","dPt vs Pt", kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra,200,-10,10);
310 fGammaList->Add(hMCPtResolution);
311 hMCPtResolutionPhi=new TH2F("Resolution_Gamma_dPt_Phi","dPt vs Phi",180,0,2*TMath::Pi(),200,-10,10);
312 fGammaList->Add(hMCPtResolutionPhi);
313 hMCRResolutionvsR=new TH2F("Resolution_dRAbs_VS_R","dR vs R", 720,0,360,100,-5,5);
314 fGammaList->Add(hMCRResolutionvsR);
315 hMCZResolutionvsZ=new TH2F("Resolution_dZAbs_VS_Z","dZ vs Z", 200,-50,50,100,-5,5);
316 fGammaList->Add(hMCZResolutionvsZ);
321 if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
323 if(fConversionGammas == NULL){
324 if(kUseAODConversionPhoton){
325 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
327 fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
329 fConversionGammas->Delete();//Reset the TClonesArray
334 fConversionGammas->SetName(Form("GammaConv_gamma"));
336 AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
337 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
340 PostData(1, fOutputList);
343 //________________________________________________________________________
344 void AliV0ReaderV1::UserExec(Option_t *){
346 fAODEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
347 fESDEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
348 if(!fAODEvent&&!fESDEvent) {
349 AliError("No Input event");
355 fMCStack = fMCEvent->Stack();}
358 if(fESDEvent)AliKFParticle::SetField(fESDEvent->GetMagneticField());
359 if(fAODEvent)AliKFParticle::SetField(fAODEvent->GetMagneticField());
361 fConversionGammas->Delete();//Reset the TClonesArray
367 if(EventIsSelected()){
370 for(fCurrentV0Index=0;fCurrentV0Index<fInputEvent->GetNumberOfV0s();fCurrentV0Index++){
376 ProcessMCGammasForEfficiency();
380 ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
382 AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
383 if (aodhandler && aodhandler->GetFillAOD()) {
384 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
388 PostData(1, fOutputList);
390 ///________________________________________________________________________
391 Bool_t AliV0ReaderV1::IsGammaCandidate(AliConversionPhotonBase *fPhotonCandidate)
394 // Fill Histos before Cuts
395 hV0CurrentFinder->Fill(fPhotonCandidate->GetPhotonMass());
396 hV0AllArmenteros->Fill(fPhotonCandidate->GetArmenterosAlpha(),fPhotonCandidate->GetArmenterosQt());
398 Bool_t passcuts=kTRUE;
400 // Gamma selection based on QT from Armenteros
401 if(fDoQtGammaSelection == kTRUE){
402 if(!ArmenterosQtCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
407 if(fPhotonCandidate->GetChi2perNDF() > fChi2CutConversion || fPhotonCandidate->GetChi2perNDF() <=0){
408 hV0CutChi2->Fill(fPhotonCandidate->GetPhotonMass());
412 // Reconstruction Acceptance Cuts
413 if(!AcceptanceCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
417 if(!TrackCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
421 if(!dEdxCuts(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
424 if(fDoPhotonAsymmetryCut == kTRUE){
425 if(!AsymmetryCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
428 //Check the pid probability
430 if(!PIDProbabilityCut(fPhotonCandidate))return kFALSE;//passcuts=kFALSE;
436 ///________________________________________________________________________
437 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(Int_t charge){
440 if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
445 // Check for sign flip
448 AliESDv0 *fCurrentV0=dynamic_cast<AliESDv0*>(fESDEvent->GetV0(fCurrentV0Index));
450 if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
451 if(!fESDEvent->GetTrack(fCurrentV0->GetNindex())||!fESDEvent->GetTrack(fCurrentV0->GetPindex()))return 0x0;
452 if((fESDEvent->GetTrack(fCurrentV0->GetPindex()))->Charge()==charge){
453 fCurrentTrackLabels[label]=fCurrentV0->GetPindex();
454 return fCurrentV0->GetParamP();}
455 if((fESDEvent->GetTrack(fCurrentV0->GetNindex()))->Charge()==charge){
456 fCurrentTrackLabels[label]=fCurrentV0->GetNindex();
457 return fCurrentV0->GetParamN();}
463 /* AliAODv0 *fCurrentV0=dynamic_cast<AliAODv0*>(fAODEvent->GetV0(fCurrentV0Index));
465 if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
466 if(!fAODEvent->GetTrack(fCurrentV0->GetNegID())||!fAODEvent->GetTrack(fCurrentV0->GetPosID()))return 0x0;
467 if((fAODEvent->GetTrack(fCurrentV0->GetPosID()))->Charge()==1){
468 fCurrentTrackLabels[label]=fCurrentV0->GetPosID();
469 return fCurrentV0->GetParamP();}
470 if((fAODEvent->GetTrack(fCurrentV0->GetNegID()))->Charge()==1){
471 fCurrentTrackLabels[label]=fCurrentV0->GetNegID();
472 return fCurrentV0->GetParamN();}
478 ///________________________________________________________________________
479 Bool_t AliV0ReaderV1::CheckV0Status()
482 AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(fCurrentV0Index));
484 printf("Requested V0 does not exist");
486 //checks if on the fly mode is set
487 if(fCurrentV0->GetOnFlyStatus()==fUseOnFlyV0Finder)return kTRUE;
491 AliAODv0 *fCurrentV0=dynamic_cast<AliAODv0*>(fAODEvent->GetV0(fCurrentV0Index));
493 AliWarning("Requested V0 does not exist");
496 //checks if on the fly mode is set
497 if(fCurrentV0->GetOnFlyStatus()==fUseOnFlyV0Finder)return kTRUE;
505 ///________________________________________________________________________
506 void AliV0ReaderV1::ProcessV0(){
509 fCurrentTrackLabels[0]=-1;
510 fCurrentTrackLabels[1]=-1;
513 // cout<<"V0ReaderV1 ProcessV0 "<<fCurrentV0Index<<endl;
515 // Get Daughter KF Particles
517 const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP();
518 const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN();
520 if(fCurrentExternalTrackParamPositive&&fCurrentExternalTrackParamNegative){
522 fCurrentNegativeKFParticle=new AliKFParticle(*(fCurrentExternalTrackParamNegative),11);
523 fCurrentPositiveKFParticle=new AliKFParticle(*(fCurrentExternalTrackParamPositive),-11);
526 //else{hV0CutLikeSign->Fill(fCurrentMotherKFCandidate->M());}// Like Sign error is already catched here
531 if(fCurrentNegativeKFParticle&&fCurrentPositiveKFParticle){
533 if(fUseConstructGamma==kTRUE){
535 fCurrentMotherKFCandidate = new AliKFConversionPhoton();
536 fCurrentMotherKFCandidate->ConstructGamma(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
538 fCurrentMotherKFCandidate = new AliKFConversionPhoton(*fCurrentNegativeKFParticle,*fCurrentPositiveKFParticle);
539 fCurrentMotherKFCandidate->SetMassConstraint(0,fNSigmaMass);
542 if(fCurrentNegativeKFParticle){delete fCurrentNegativeKFParticle;
543 fCurrentNegativeKFParticle=0x0;}
544 if(fCurrentPositiveKFParticle){ delete fCurrentPositiveKFParticle;
545 fCurrentPositiveKFParticle=0x0;}
549 if(fUseImprovedVertex == kTRUE){
550 AliKFVertex primaryVertexImproved(*GetPrimaryVertex());
551 primaryVertexImproved+=*fCurrentMotherKFCandidate;
552 fCurrentMotherKFCandidate->SetProductionVertex(primaryVertexImproved);
557 fCurrentMotherKFCandidate->SetV0Index(fCurrentV0Index);
558 fCurrentMotherKFCandidate->SetTrackLabels(fCurrentTrackLabels[0],fCurrentTrackLabels[1]);
563 Int_t labeln=TMath::Abs(GetTrack(fCurrentMotherKFCandidate->GetTrackLabelPositive())->GetLabel());
564 Int_t labelp=TMath::Abs(GetTrack(fCurrentMotherKFCandidate->GetTrackLabelNegative())->GetLabel());
566 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
567 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
569 if(fPositiveMCParticle&&fNegativeMCParticle){
570 fCurrentMotherKFCandidate->SetMCLabelPositive(labelp);
571 fCurrentMotherKFCandidate->SetMCLabelNegative(labeln);
577 //Add PID information with ESD tender (AOD implementation is not complete)
580 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
581 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
582 AliPIDResponse *fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
585 Int_t labelp=((AliESDv0*)fESDEvent->GetV0(fCurrentMotherKFCandidate->GetV0Index()))->GetNindex();
586 Int_t labeln=((AliESDv0*)fESDEvent->GetV0(fCurrentMotherKFCandidate->GetV0Index()))->GetNindex();
588 AliESDtrack *trackpos=fESDEvent->GetTrack(labelp);
589 AliESDtrack *trackneg=fESDEvent->GetTrack(labeln);
591 if(trackpos&&trackneg){
593 Float_t fNSigmadEdxPositive[5];
594 Float_t fNSigmadEdxNegative[5];
596 fNSigmadEdxPositive[0]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kElectron);
597 fNSigmadEdxPositive[1]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kMuon);
598 fNSigmadEdxPositive[2]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kPion);
599 fNSigmadEdxPositive[3]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kKaon);
600 fNSigmadEdxPositive[4]=fPIDResponse->NumberOfSigmasTPC(trackpos,AliPID::kProton);
602 fNSigmadEdxNegative[0]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kElectron);
603 fNSigmadEdxNegative[1]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kMuon);
604 fNSigmadEdxNegative[2]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kPion);
605 fNSigmadEdxNegative[3]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kKaon);
606 fNSigmadEdxNegative[4]=fPIDResponse->NumberOfSigmasTPC(trackneg,AliPID::kProton);
608 fCurrentMotherKFCandidate->SetNSigmadEdx(fNSigmadEdxPositive,fNSigmadEdxNegative);
612 // Calculate ConversionPoint
614 if(fUseOwnXYZCalculation){
616 Double_t convpos[3]={0,0,0};
617 GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos);
618 fCurrentMotherKFCandidate->SetConversionPoint(convpos);
623 if(IsGammaCandidate(fCurrentMotherKFCandidate)){
624 // Fill Histos after Cuts
628 ProcessMC(fCurrentMotherKFCandidate);
631 hV0Good->Fill(fCurrentMotherKFCandidate->M());
632 hV0GoodArmenteros->Fill(fCurrentMotherKFCandidate->GetArmenterosAlpha(),fCurrentMotherKFCandidate->GetArmenterosQt());
635 // Set Mass Zero for Gammas
638 // Add Gamma to the TClonesArray
640 if(kUseAODConversionPhoton){
641 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
644 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
648 hGammaPhi->Fill(fCurrentMotherKFCandidate->Phi());
649 hGammaPt[fCentralityBin]->Fill(fCurrentMotherKFCandidate->Pt());
650 hGammaConversionMapXY->Fill(fCurrentMotherKFCandidate->GetConversionX(),fCurrentMotherKFCandidate->GetConversionY());
651 hGammaConversionMapZR->Fill(fCurrentMotherKFCandidate->GetConversionZ(),fCurrentMotherKFCandidate->GetConversionRadius());
654 delete fCurrentMotherKFCandidate;
655 fCurrentMotherKFCandidate=NULL;
660 ///________________________________________________________________________
661 Bool_t AliV0ReaderV1::ArmenterosQtCut(AliConversionPhotonBase *fPhotonCandidate)
663 if(fDoHighPtQtGammaSelection){
664 if(fPhotonCandidate->GetPhotonPt() < fPtBorderForQt){
665 if(fPhotonCandidate->GetArmenterosQt()>fQtMax){
666 hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
670 if(fPhotonCandidate->GetArmenterosQt()>fHighPtQtMax){
671 hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
677 if(fPhotonCandidate->GetArmenterosQt()>fQtMax){
678 hV0CutQt->Fill(fPhotonCandidate->GetPhotonMass());
685 ///________________________________________________________________________
686 Bool_t AliV0ReaderV1::AcceptanceCuts(AliConversionPhotonBase *fPhotonCandidate)
688 AliVTrack *fCurrentNegativeTrack=GetTrack(fPhotonCandidate->GetTrackLabelNegative());
689 AliVTrack *fCurrentPositiveTrack=GetTrack(fPhotonCandidate->GetTrackLabelPositive());
691 if(fPhotonCandidate->GetConversionRadius()>fMaxR){ // cuts on distance from collision point
692 hV0CutR->Fill(fPhotonCandidate->GetPhotonMass());
696 if(fPhotonCandidate->GetConversionRadius()<fMinR){ // cuts on distance from collision point
697 hV0CutMinR->Fill(fPhotonCandidate->GetPhotonMass());
701 if(fPhotonCandidate->GetConversionRadius() <= ((TMath::Abs(fPhotonCandidate->GetConversionZ())*fLineCutZRSlope)-fLineCutZValue)){
702 hV0CutLine->Fill(fPhotonCandidate->GetPhotonMass());
705 /*else if (fUseEtaMinCut && fPhotonCandidate->GetConversionRadius() >= ((TMath::Abs(fPhotonCandidate->GetConversionZ())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
706 hV0CutLine->Fill(fPhotonCandidate->GetPhotonMass());
711 if(TMath::Abs(fPhotonCandidate->GetConversionZ()) > fMaxZ ){ // cuts out regions where we do not reconstruct
712 hV0CutZ->Fill(fPhotonCandidate->GetPhotonMass());
716 if(TMath::Abs(fPhotonCandidate->GetPhotonEta())> fEtaCut || TMath::Abs(fPhotonCandidate->GetPhotonEta())< fEtaCutMin){
717 hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
721 if(TMath::Abs(fCurrentNegativeTrack->Eta())> fEtaCut || TMath::Abs(fCurrentNegativeTrack->Eta())< fEtaCutMin){
722 hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
726 if(TMath::Abs(fCurrentPositiveTrack->Eta())> fEtaCut || TMath::Abs(fCurrentPositiveTrack->Eta())< fEtaCutMin){
727 hV0CutEta->Fill(fPhotonCandidate->GetPhotonMass());
731 if( fCurrentNegativeTrack->Pt()< fSinglePtCut || fCurrentNegativeTrack->Pt()< fSinglePtCut){
732 hV0CutSinglePt->Fill(fPhotonCandidate->GetPhotonMass());
737 if(fPhotonCandidate->GetPhotonPt()<fPtCut){
738 hV0CutPt->Fill(fPhotonCandidate->GetPhotonMass());
746 ///________________________________________________________________________
747 Bool_t AliV0ReaderV1::TrackCuts(AliConversionPhotonBase *fPhotonCandidate){
749 Bool_t passtrackcuts=kTRUE;
755 AliESDtrack *fCurrentNegativeESDTrack=(AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
756 AliESDtrack *fCurrentPositiveESDTrack=(AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
758 if(!fCurrentNegativeESDTrack||!fCurrentPositiveESDTrack)return kFALSE;
761 if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){
763 hV0CutLikeSign->Fill(fPhotonCandidate->GetPhotonMass());
764 passtrackcuts=kFALSE;
767 if( (!(fCurrentNegativeESDTrack->IsOn(AliESDtrack::kTPCrefit))||(!(fCurrentPositiveESDTrack->IsOn(AliESDtrack::kTPCrefit))))){
768 hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
769 passtrackcuts=kFALSE;
773 if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
774 fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
775 passtrackcuts=kFALSE;
778 if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
779 passtrackcuts=kFALSE;
780 hV0CutKinks->Fill(fPhotonCandidate->GetPhotonMass());
783 /* Double_t negclsToF = 0.;
784 if (!fUseCorrectedTPCClsInfo ){
785 if(fCurrentNegativeESDTrack->GetTPCNclsF()!=0 ){
786 negclsToF = (Double_t)fCurrentNegativeESDTrack->GetNcls(1)/(Double_t)fCurrentNegativeESDTrack->GetTPCNclsF();
789 negclsToF = fCurrentNegativeESDTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
792 Double_t posclsToF = 0.;
793 if (!fUseCorrectedTPCClsInfo ){
794 if(fCurrentTrack->GetTPCNclsF()!=0 ){
795 posclsToF = (Double_t)fCurrentTrack->GetNcls(1)/(Double_t)fCurrentTrack->GetTPCNclsF();
798 posclsToF = fCurrentTrack->GetTPCClusterInfo(2,0,GetFirstTPCRow(GetXYRadius()));
801 if( negclsToF < fMinClsTPCToF || posclsToF < fMinClsTPCToF ){
802 hV0CutMinNclsTPCToF->Fill(fPhotonCandidate->GetPhotonMass());
803 passtrackcuts=kFALSE; }
813 AliAODTrack *fCurrentNegativeESDTrack=(AliAODTrack*)fAODEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
814 AliAODTrack *fCurrentPositiveESDTrack=(AliAODTrack*)fAODEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
816 if(!fCurrentNegativeESDTrack||!fCurrentPositiveESDTrack)return kFALSE;
819 if(fCurrentNegativeESDTrack->Charge() == fCurrentPositiveESDTrack->Charge()){
821 hV0CutLikeSign->Fill(fPhotonCandidate->GetPhotonMass());
822 passtrackcuts=kFALSE;
825 if( !(fCurrentNegativeESDTrack->IsOn(AliESDtrack::kTPCrefit))){
826 hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
827 passtrackcuts=kFALSE;
830 if( !(fCurrentPositiveESDTrack->IsOn(AliESDtrack::kTPCrefit))){
831 hV0CutRefit->Fill(fPhotonCandidate->GetPhotonMass());
832 passtrackcuts=kFALSE;
837 if( fCurrentNegativeESDTrack->GetKinkIndex(0) > 0 ||
838 fCurrentPositiveESDTrack->GetKinkIndex(0) > 0) {
841 if(fCurrentNegativeESDTrack->GetNcls(1) < fMinClsTPC || fCurrentPositiveESDTrack->GetNcls(1) < fMinClsTPC ){
842 passtrackcuts=kFALSE;}
847 return passtrackcuts;
850 ///________________________________________________________________________
851 Bool_t AliV0ReaderV1::dEdxCuts(AliConversionPhotonBase *fPhotonCandidate){
853 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
854 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
855 AliPIDResponse *fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
857 AliVTrack *fCurrentTrack=0x0;
859 for(Int_t ilabel=0;ilabel<2;ilabel++){
861 fCurrentTrack=GetTrack(fPhotonCandidate->GetTrackLabel(ilabel));
863 if(!fCurrentTrack)return kFALSE;
865 // Fill dEdx before cuts
867 hGammadEdxbefore->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
869 if(fDodEdxSigmaCut == kTRUE){
870 if( fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaBelowElectronLine ||
871 fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaAboveElectronLine){
873 hV0CutdEdxElectron->Fill(fPhotonCandidate->GetPhotonMass());
877 if( fCurrentTrack->P()>fPIDMinPnSigmaAbovePionLine && fCurrentTrack->P()<fPIDMaxPnSigmaAbovePionLine ){
878 if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
879 fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
880 fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLine){
882 hV0CutdEdxPion->Fill(fPhotonCandidate->GetPhotonMass());
888 if( fCurrentTrack->P()>fPIDMaxPnSigmaAbovePionLine ){
889 if(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)>fPIDnSigmaBelowElectronLine &&
890 fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kElectron)<fPIDnSigmaAboveElectronLine&&
891 fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion)<fPIDnSigmaAbovePionLineHighPt){
893 hV0CutdEdxPion->Fill(fPhotonCandidate->GetPhotonMass());
899 if(fDoKaonRejectionLowP == kTRUE){
900 if(fCurrentTrack->P()<fPIDMinPKaonRejectionLowP ){
901 if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kKaon))<fPIDnSigmaAtLowPAroundKaonLine){
902 hV0CutdEdxKaonLowP->Fill(fPhotonCandidate->GetPhotonMass());
908 if(fDoProtonRejectionLowP == kTRUE){
909 if( fCurrentTrack->P()<fPIDMinPProtonRejectionLowP ){
910 if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kProton))<fPIDnSigmaAtLowPAroundProtonLine){
911 hV0CutdEdxProtonLowP->Fill(fPhotonCandidate->GetPhotonMass());
917 if(fDoPionRejectionLowP == kTRUE){
918 if( fCurrentTrack->P()<fPIDMinPPionRejectionLowP ){
919 if( TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fCurrentTrack,AliPID::kPion))<fPIDnSigmaAtLowPAroundPionLine){
920 hV0CutdEdxPionLowP->Fill(fPhotonCandidate->GetPhotonMass());
927 if( fDoTOFsigmaCut == kTRUE ){ // RRnewTOF start /////////////////////////////////////////////////////////////////////////////
929 if((fPIDResponse->NumberOfSigmasTOF(fCurrentTrack,AliPID::kElectron)>fTofPIDnSigmaAboveElectronLine) || (fPIDResponse->NumberOfSigmasTOF(fCurrentTrack,AliPID::kElectron)<fTofPIDnSigmaBelowElectronLine)){
930 hV0CutdEdxTOFElectron->Fill(fPhotonCandidate->GetPhotonMass());
933 } /////////////////////////////// RRnewTOF end ///////////////////////////////////////////////////////////////////////////////
937 if(!fPIDResponse->IdentifiedAsElectronTRD(fCurrentTrack,fPIDTRDEfficiency)){
938 hV0CutdEdxTRD->Fill(fPhotonCandidate->GetPhotonMass());
943 // Fill dEdx Histogram after Cuts
945 hGammadEdxafter->Fill(fCurrentTrack->P(),fCurrentTrack->GetTPCsignal());
952 ///________________________________________________________________________
953 Bool_t AliV0ReaderV1::AsymmetryCut(AliConversionPhotonBase *fPhotonCandidate)
956 for(Int_t ii=0;ii<2;ii++){
957 AliVTrack *fCurrentTrack=GetTrack(fPhotonCandidate->GetTrackLabel(ii));
959 if( fCurrentTrack->P()>fMinPPhotonAsymmetryCut ){
960 Double_t trackNegAsy=0;
961 if (fPhotonCandidate->GetPhotonP()!=0.){
962 trackNegAsy= fCurrentTrack->P()/fPhotonCandidate->GetPhotonP();
964 if( trackNegAsy<fMinPhotonAsymmetry ||trackNegAsy>(1.- fMinPhotonAsymmetry)){
965 hV0CutAsymmetry->Fill(fPhotonCandidate->GetPhotonMass());
975 ///________________________________________________________________________
977 Int_t AliV0ReaderV1::GetNumberOfContributorsVtx(){
978 // returns number of contributors to the vertex
980 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()>0) {
981 return fESDEvent->GetPrimaryVertexTracks()->GetNContributors();
984 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors()<1) {
986 //-AM test pi0s without SPD only vertex
987 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
988 return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
991 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
992 // cout<<"number of contributors from bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
998 if(fAODEvent->GetPrimaryVertex()->GetNContributors()>0) {
999 return fESDEvent->GetPrimaryVertex()->GetNContributors();
1001 if(fAODEvent->GetPrimaryVertex()->GetNContributors()<1) {
1002 if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
1003 return fAODEvent->GetPrimaryVertexSPD()->GetNContributors();
1006 if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()<1) {
1007 AliWarning(Form("Number of contributors from bad vertex type:: %s",fAODEvent->GetPrimaryVertex()->GetName()));
1017 ///________________________________________________________________________
1018 Double_t AliV0ReaderV1::GetCentrality(){
1020 if(fAODEvent->GetHeader()){return fAODEvent->GetHeader()->GetCentrality();}
1024 return ((AliCentrality *)fESDEvent->GetCentrality())->GetCentralityPercentile("V0M");
1030 ///________________________________________________________________________
1032 AliEventplane *AliV0ReaderV1::GetEventPlane(){
1034 return fESDEvent->GetEventplane();
1037 return fAODEvent->GetEventplane();
1042 //________________________________________________________________________
1043 Bool_t AliV0ReaderV1::SetEventPlane()
1045 if(fMCStack||!fIsHeavyIon){
1046 // NO EP in MC mode and pp mode
1052 AliEventplane *fEP=GetEventPlane();
1054 if(!fEP)return kFALSE;
1056 fEPAngle=fEP->GetEventplane("Q");
1063 ///________________________________________________________________________
1064 Bool_t AliV0ReaderV1::EventCuts(){
1066 fEventIsSelected=kTRUE;
1068 // Z Vertex Position Cut
1071 hV0EventCuts->Fill(0);
1072 fEventIsSelected=kFALSE;;
1074 // Number of Contributors Cut
1076 if(GetNumberOfContributorsVtx()<=0) {
1077 hV0EventCuts->Fill(1);
1078 fEventIsSelected=kFALSE;;
1080 // Centrality Selection
1082 if(!CentralitySelection()){
1083 hV0EventCuts->Fill(2);
1084 fEventIsSelected=kFALSE;;
1088 if(!SetEventPlane()){hV0EventCuts->Fill(4);}
1090 // Fill Event Histograms
1092 if(fEventIsSelected){
1093 // Fill Event Histograms
1094 hV0EventCuts->Fill(9);
1095 hVertexZ->Fill(fVertexZ);
1096 hCentrality->Fill(fCentrality);
1097 hNEvents->Fill(fCentralityBin);
1099 else{hV0EventCuts->Fill(8);}
1101 return fEventIsSelected;
1104 ///________________________________________________________________________
1105 Bool_t AliV0ReaderV1::VertexZCut(){
1107 fVertexZ=GetPrimaryVertex()->GetZ();
1110 if(fBGHandler->GetZBinIndex(fVertexZ)<0)return kFALSE;
1113 if(fVertexZ>fMaxVertexZ)return kFALSE;
1118 ///________________________________________________________________________
1119 AliVTrack *AliV0ReaderV1::GetTrack(Int_t label){
1121 return (AliESDtrack*)fESDEvent->GetTrack(label);
1123 if(fAODEvent)return (AliAODTrack*)fAODEvent->GetTrack(label);
1127 ///________________________________________________________________________
1128 Bool_t AliV0ReaderV1::PIDProbabilityCut(AliConversionPhotonBase *fPhotonCandidate){
1132 Bool_t iResult=kFALSE;
1134 Double_t *posProbArray = new Double_t[AliPID::kSPECIES];
1135 Double_t *negProbArray = new Double_t[AliPID::kSPECIES];
1137 AliESDtrack* negTrack = (AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelNegative());
1138 AliESDtrack* posTrack = (AliESDtrack*)fESDEvent->GetTrack(fPhotonCandidate->GetTrackLabelPositive());
1140 if(negProbArray && posProbArray){
1142 negTrack->GetTPCpid(negProbArray);
1143 posTrack->GetTPCpid(posProbArray);
1145 if(negProbArray[AliPID::kElectron]>=fPIDProbabilityCutNegativeParticle && posProbArray[AliPID::kElectron]>=fPIDProbabilityCutPositiveParticle){
1148 else{hV0CutPIDProb->Fill(fPhotonCandidate->GetPhotonMass());}
1151 delete [] posProbArray;
1152 delete [] negProbArray;
1157 // not possible to implement
1162 ///________________________________________________________________________
1163 const AliVertex *AliV0ReaderV1::GetPrimaryVertex() {
1164 if(fESDEvent)return fESDEvent->GetPrimaryVertex();
1165 if(fAODEvent)return const_cast<const AliVertex*>(dynamic_cast<AliVertex*>(fAODEvent->GetPrimaryVertex()));
1169 ///________________________________________________________________________
1170 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track, Double_t b,Int_t charge, Double_t center[2]){
1171 // see header file for documentation
1174 track->GetHelixParameters(helix,b);
1176 Double_t xpos = helix[5];
1177 Double_t ypos = helix[0];
1178 Double_t radius = TMath::Abs(1./helix[4]);
1179 Double_t phi = helix[2];
1182 phi = phi + 2*TMath::Pi();
1185 phi -= TMath::Pi()/2.;
1186 Double_t xpoint = radius * TMath::Cos(phi);
1187 Double_t ypoint = radius * TMath::Sin(phi);
1211 center[0] = xpos + xpoint;
1212 center[1] = ypos + ypoint;
1216 ///________________________________________________________________________
1217 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3]){
1219 if(!pparam||!nparam)return kFALSE;
1221 Double_t helixcenterpos[2];
1222 GetHelixCenter(pparam,GetMagneticField(),pparam->Charge(),helixcenterpos);
1224 Double_t helixcenterneg[2];
1225 GetHelixCenter(nparam,GetMagneticField(),nparam->Charge(),helixcenterneg);
1227 Double_t helixpos[6];
1228 pparam->GetHelixParameters(helixpos,GetMagneticField());
1229 Double_t posradius = TMath::Abs(1./helixpos[4]);
1231 Double_t helixneg[6];
1232 nparam->GetHelixParameters(helixneg,GetMagneticField());
1233 Double_t negradius = TMath::Abs(1./helixneg[4]);
1235 // Calculate xy-position
1237 Double_t xpos = helixcenterpos[0];
1238 Double_t ypos = helixcenterpos[1];
1239 Double_t xneg = helixcenterneg[0];
1240 Double_t yneg = helixcenterneg[1];
1242 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
1243 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
1246 // Calculate z-position
1248 Double_t deltaXPos = convpos[0] - xpos;
1249 Double_t deltaYPos = convpos[1] - ypos;
1251 Double_t deltaXNeg = convpos[0] - xneg;
1252 Double_t deltaYNeg = convpos[1] - yneg;
1254 Double_t alphaPos = TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
1255 Double_t alphaNeg = TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
1257 Double_t vertexXNeg = xneg + TMath::Abs(negradius)*
1258 TMath::Cos(alphaNeg);
1259 Double_t vertexYNeg = yneg + TMath::Abs(negradius)*
1260 TMath::Sin(alphaNeg);
1262 Double_t vertexXPos = xpos + TMath::Abs(posradius)*
1263 TMath::Cos(alphaPos);
1264 Double_t vertexYPos = ypos + TMath::Abs(posradius)*
1265 TMath::Sin(alphaPos);
1267 Double_t x0neg = helixneg[5];
1268 Double_t y0neg = helixneg[0];
1270 Double_t x0pos = helixpos[5];
1271 Double_t y0pos = helixpos[0];
1273 Double_t dNeg = TMath::Sqrt((vertexXNeg - x0neg)*(vertexXNeg - x0neg)
1274 +(vertexYNeg - y0neg)*(vertexYNeg - y0neg));
1276 Double_t dPos = TMath::Sqrt((vertexXPos - x0pos)*(vertexXPos - x0pos)
1277 +(vertexYPos - y0pos)*(vertexYPos - y0pos));
1279 Double_t rNeg = TMath::Sqrt(negradius*negradius -
1282 Double_t rPos = TMath::Sqrt(posradius*posradius -
1285 Double_t deltabetaNeg = 2*(TMath::Pi() + TMath::ATan2(-dNeg/2.,-rNeg));
1286 Double_t deltabetaPos = 2*(TMath::Pi() + TMath::ATan2(-dPos/2.,-rPos));
1288 Double_t deltaUNeg = negradius*deltabetaNeg;
1289 Double_t deltaUPos = posradius*deltabetaPos;
1291 Double_t zphaseNeg = nparam->GetZ() + deltaUNeg * nparam->GetTgl();
1292 Double_t zphasePos = pparam->GetZ() + deltaUPos * pparam->GetTgl();
1294 convpos[2] = (zphasePos*negradius+zphaseNeg*posradius)/(negradius+posradius);
1298 ///________________________________________________________________________
1299 void AliV0ReaderV1::ProcessMC(AliKFConversionPhoton *fCurrentReconstructedGamma){
1301 if(!fMCStack)return;
1303 TParticle *fMotherMCParticle=NULL;
1304 TParticle *fNegativeMCParticle=NULL;
1305 TParticle *fPositiveMCParticle=NULL;
1308 fMotherMCParticle = fCurrentReconstructedGamma->GetMCParticle(fMCStack);
1309 fNegativeMCParticle = fCurrentReconstructedGamma->GetNegativeMCDaughter(fMCStack);
1310 fPositiveMCParticle = fCurrentReconstructedGamma->GetPositiveMCDaughter(fMCStack);
1312 if(fPositiveMCParticle&&fNegativeMCParticle&&fMotherMCParticle){
1314 // Check if it is a true photon
1316 if(fMotherMCParticle->GetPdgCode()==22){
1318 hMCPtRECOTRUE[fCentralityBin]->Fill(fCurrentReconstructedGamma->GetPt());
1322 Double_t mcpt = fMotherMCParticle->Pt();
1323 Double_t esdpt = fCurrentReconstructedGamma->GetPt();
1324 Double_t resdPt = 0.;
1326 resdPt = ((esdpt - mcpt)/mcpt)*100.;
1327 } else if(mcpt < 0){
1328 AliWarning("Pt of MC particle is negative, this will cause wrong calculation of resPt");
1332 hMCPtResolution->Fill(mcpt,resdPt);
1333 hMCPtResolutionPhi->Fill(fMotherMCParticle->Phi(),resdPt);
1335 // Conversion Point Resolution
1337 Double_t resdR = 0.;
1338 if(fNegativeMCParticle->R() != 0){
1339 resdR = ((fCurrentReconstructedGamma->GetConversionRadius() - fNegativeMCParticle->R())/fNegativeMCParticle->R())*100.;
1341 Double_t resdRAbs = 0.;
1342 resdRAbs = (fCurrentReconstructedGamma->GetConversionRadius() - fNegativeMCParticle->R());
1344 hMCRResolutionvsR->Fill(fNegativeMCParticle->R(),resdRAbs);
1346 // fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1347 // fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1348 // fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1349 // fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1351 Double_t resdZ = 0.;
1352 if(fNegativeMCParticle->Vz() != 0){
1353 resdZ = ((fCurrentReconstructedGamma->GetZ() -fNegativeMCParticle->Vz())/fNegativeMCParticle->Vz())*100.;
1355 Double_t resdZAbs = 0.;
1356 resdZAbs = fCurrentReconstructedGamma->GetZ() -fNegativeMCParticle->Vz();
1358 hMCZResolutionvsZ->Fill( fNegativeMCParticle->Vz(), resdZAbs);
1364 ///________________________________________________________________________
1365 void AliV0ReaderV1::ProcessMCGammasForEfficiency(){
1367 if(!fMCStack)return;
1369 for (Int_t iTracks = 0; iTracks < fMCStack->GetNprimary(); iTracks++) {
1370 TParticle* particle = (TParticle *)fMCStack->Particle(iTracks);
1372 //process the gammas
1374 if(IsMCConversionGammaInAcceptance(particle)){
1376 hMCPtTRUE[fCentralityBin]->Fill(particle->Pt());
1383 ///________________________________________________________________________
1384 Bool_t AliV0ReaderV1::IsMCConversionGammaInAcceptance(TParticle *particle){
1385 if(!fMCStack)return kFALSE;
1387 if (particle->GetPdgCode() == 22){
1388 if(TMath::Abs(particle->Eta())> fEtaCut || TMath::Abs(particle->Eta())< fEtaCutMin) return kFALSE;
1390 if(particle->GetMother(0) >-1 && fMCStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
1391 return kFALSE; // no photon as mothers!
1394 if(particle->GetMother(0) >= fMCStack->GetNprimary()){
1395 return kFALSE; // the gamma has a mother, and it is not a primary particle
1398 // looking for conversion (electron + positron from pairbuilding (= 5) )
1399 TParticle* ePos = NULL;
1400 TParticle* eNeg = NULL;
1402 if(particle->GetNDaughters() >= 2){
1403 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1404 TParticle *tmpDaughter = fMCStack->Particle(daughterIndex);
1405 if(tmpDaughter->GetUniqueID() == 5){
1406 if(tmpDaughter->GetPdgCode() == 11){
1408 } else if(tmpDaughter->GetPdgCode() == -11){
1415 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
1419 if(AcceptanceCut(particle,ePos,eNeg))return kTRUE;
1424 ///________________________________________________________________________
1425 Bool_t AliV0ReaderV1::AcceptanceCut(TParticle *particle, TParticle * ePos,TParticle* eNeg){
1427 // cuts on distance from collision point
1429 if(particle->R()>fMaxR){
1432 if(ePos->R()>fMaxR){
1436 if(ePos->R()<fMinR){
1440 if( ePos->R() <= ((TMath::Abs(ePos->Vz())*fLineCutZRSlope)-fLineCutZValue)){
1443 /*else if (fUseEtaMinCut && ePos->R() >= ((TMath::Abs(ePos->Vz())*fLineCutZRSlopeMin)-fLineCutZValueMin )){
1447 if(TMath::Abs(eNeg->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
1451 if(eNeg->Vz()!=ePos->Vz()||eNeg->R()!=ePos->R()){
1455 if(TMath::Abs(ePos->Vz()) > fMaxZ ){ // cuts out regions where we do not reconstruct
1459 if(TMath::Abs(particle->Eta())> fEtaCut || TMath::Abs(particle->Eta())< fEtaCutMin){
1463 if(TMath::Abs(ePos->Eta())> fEtaCut || TMath::Abs(ePos->Eta())< fEtaCutMin){
1467 if(TMath::Abs(eNeg->Eta())> fEtaCut || TMath::Abs(eNeg->Eta())< fEtaCutMin){
1471 if( ePos->Pt()< fSinglePtCut || eNeg->Pt()< fSinglePtCut){
1475 if(particle->Pt()<fPtCut){
1481 ///________________________________________________________________________
1482 void AliV0ReaderV1::PrintCuts(){
1484 cout<<"V0 Reader initialized with following settings"<<endl;
1487 cout<<"Acceptance Eta:"<<endl;
1488 cout<<fEtaCutMin<<" < eta < "<<fEtaCut<<endl;
1489 cout<<"Conversion Point"<<endl;
1490 cout<<"Z <"<<fMaxZ<<endl;
1491 cout<<fMinR<<" < R < "<<fMaxR<<endl;
1492 cout<<"Line Cut Slope"<<fLineCutZRSlope<<" ZValue "<<fLineCutZValue<<endl;
1494 cout<<"Pt Gamma > "<<fPtCut<<endl;
1495 cout<<"Pt Daughters > "<<fSinglePtCut<<endl;
1497 cout<<"Armenteros Qt Cut"<<endl;
1499 if(fDoHighPtQtGammaSelection){
1500 cout<<" qt < "<<fQtMax<<" for pt < "<<fPtBorderForQt<<endl;
1501 cout<<" qt < "<<fHighPtQtMax<<" for pt > "<<fPtBorderForQt<<endl;
1504 cout<<" qt < "<<fQtMax<<endl;
1506 cout<<"Chi2perNDF > "<<fChi2CutConversion<<endl;
1512 //_______________________________________________________________________
1514 Bool_t AliV0ReaderV1::CentralitySelection(){
1518 if(!fIsHeavyIon){fCentralityBin=0;fCentrality=0;return kTRUE;}
1520 fCentrality=GetCentrality();
1525 fCentralityBin=fBGHandler->GetCentralityBinIndex(Int_t(fCentrality));
1528 Double_t fCentralityBins[fNCentralityBins+1];
1529 for(int i=0;i<fNCentralityBins;i++){
1530 fCentralityBins[i]=i*100/Double_t(fNCentralityBins);
1533 for(int i=0;i<fNCentralityBins;i++){
1534 if(fCentrality>fCentralityBins[i]&&fCentrality<fCentralityBins[i+1]){
1541 if(fCentralityBin>=0&&fCentrality>=0){
1546 AliWarning("Centrality not defined");
1551 //________________________________________________________________________
1552 void AliV0ReaderV1::Terminate(Option_t *)
1555 printf(Form("V0ReaderV1: V0s processed: %4.0f reconstructed photons: %4.0f \n",hV0CurrentFinder->GetEntries(),hV0Good->GetEntries()));