1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Pedro González, Pedro Ladrón de Guevara, Ernesto López Torres, *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 // Analysis task for pi0->e+e-gamma (Dalitz decay)
21 #include "TParticle.h"
23 #include "TMCProcess.h"
24 #include "TDatabasePDG.h"
27 #include "TDirectory.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDInputHandler.h"
32 #include "AliESDtrack.h"
33 #include "AliMCEvent.h"
35 #include "AliMCEventHandler.h"
38 #include "AliESDtrackCuts.h"
39 #include "AliESDpidCuts.h"
40 #include "AliMCEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDpid.h"
44 #include "AliKFParticle.h"
45 #include "AliMCEventHandler.h"
46 #include "AliGammaConversionHistograms.h"
47 #include "AliV0Reader.h"
48 #include "AliKFVertex.h"
50 #include "AliAnalysisTaskGammaConvDalitz.h"
53 ClassImp( AliAnalysisTaskGammaConvDalitz )
55 //-----------------------------------------------------------------------------------------------
56 AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz():
61 fEposCandidateIndex(),
62 fEnegCandidateIndex(),
63 fGammaCandidatePosIndex(),
64 fGammaCandidateNegIndex(),
82 fUseTrackIndexCut(kTRUE),
83 fUsePsiPairCut(kTRUE),
86 fReadMagFieldSign(kTRUE),
89 fkElectronMass(0.00051099891),
92 fDeltaPhiCutMax(0.12),
95 fNSigmaBelowElecTPCbethe(-2.),
96 fNSigmaAboveElecTPCbethe(3.),
97 fNSigmaAbovePionTPCbethe(3.),
98 fNSigmaAboveKaonTPCbethe(3.),
99 fNSigmaAboveProtonTPCbethe(3.),
100 fTrkSelectionCriteria(kGlobalTrack)
103 // Default constructor
105 AdoptITSsaTrackCuts();
109 fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
110 fGammaPool->SetOwner(kTRUE);
113 //-----------------------------------------------------------------------------------------------
114 AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz( const char* name ):
115 AliAnalysisTaskSE( name ),
119 fEposCandidateIndex(),
120 fEnegCandidateIndex(),
121 fGammaCandidatePosIndex(),
122 fGammaCandidateNegIndex(),
139 fUseBayesPID(kFALSE),
140 fUseTrackIndexCut(kTRUE),
141 fUsePsiPairCut(kTRUE),
143 fUseGammaCut(kFALSE),
144 fReadMagFieldSign(kTRUE),
147 fkElectronMass(0.00051099891),
150 fDeltaPhiCutMax(0.12),
153 fNSigmaBelowElecTPCbethe(-2.),
154 fNSigmaAboveElecTPCbethe(3.),
155 fNSigmaAbovePionTPCbethe(3.),
156 fNSigmaAboveKaonTPCbethe(3.),
157 fNSigmaAboveProtonTPCbethe(3.),
158 fTrkSelectionCriteria(kGlobalTrack)
160 // Common I/O in slot 0
161 DefineInput (0, TChain::Class());
163 // Your private output
164 DefineOutput(1, TList::Class());
165 // DefineOutput(2, AliCFContainer::Class()); // for CF
167 AdoptITSsaTrackCuts();
170 // fkElectronMass = TDatabasePDG::Instance()->GetParticle( ::kElectron )->Mass(); //
172 fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
173 fGammaPool->SetOwner(kTRUE);
176 //-----------------------------------------------------------------------------------------------
177 AliAnalysisTaskGammaConvDalitz::~AliAnalysisTaskGammaConvDalitz()
180 // virtual destructor
183 if( fOutputContainer ) delete fOutputContainer;
184 if( fHistograms ) delete fHistograms;
185 if( fStandalone && fV0Reader ) delete fV0Reader;
186 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
187 if( fESDtrackCuts ) delete fESDtrackCuts;
188 if( fESDpidCuts ) delete fESDpidCuts;
189 if( fGammaCandidates) delete fGammaCandidates;
190 if( fGammaPool ) delete fGammaPool;
193 //-----------------------------------------------------------------------------------------------
194 void AliAnalysisTaskGammaConvDalitz::ConnectInputData(Option_t *option)
197 // Connect Input Data
199 if( fDebug ) AliInfo("=> ConnectInputData");
201 AliAnalysisTaskSE::ConnectInputData(option);
205 AliFatal("There is not pointer to AliV0Reader object!!!");
209 //-----------------------------------------------------------------------------------------------
210 void AliAnalysisTaskGammaConvDalitz::UserCreateOutputObjects()
213 // Create ouput objects
215 if( fDebug ) AliInfo("=> UserCreateOutputObjects");
217 // Create the output container
218 if( fOutputContainer != 0 )
220 delete fOutputContainer;
223 fOutputContainer = new TList();
225 // Add the histograms to the output container
226 fHistograms->GetOutputContainer( fOutputContainer );
227 fOutputContainer->SetOwner(kTRUE);
229 PostData( 1, fOutputContainer );
232 //-----------------------------------------------------------------------------------------------
233 void AliAnalysisTaskGammaConvDalitz::UserExec(Option_t */*option*/)
236 // Execute analysis for current event
238 if( fDebug ) AliInfo("=> UserExec");
242 AliFatal("no pointer to AliV0Reader");
246 // Create list of gamma candidates in standalone mode
247 // otherwise use the created ones by AliAnalysisTaskGammaConversion
251 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
252 AliESDInputHandler *esdHandler=0;
253 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
254 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
256 //load esd pid bethe bloch parameters depending on the existance of the MC handler
257 // yes: MC parameters
258 // no: data parameters
259 if (!AliV0Reader::GetESDpid()){
261 AliV0Reader::InitESDpid();
263 AliV0Reader::InitESDpid(1);
271 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
272 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
275 AliError("Could not retrive MC event handler!");
279 if (!mcHandler->InitOk() ){
283 if (!mcHandler->TreeK() ){
286 if (!mcHandler->TreeTR() ) {
294 fV0Reader->SetInputAndMCEvent( InputEvent(), MCEvent() );
295 fV0Reader->Initialize();
298 if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
300 if( fDebug ) AliInfo("no contributors to primary vertex");
304 if( fV0Reader->CheckForPrimaryVertexZ() == kFALSE )
307 if( fDebug ) AliInfo("z vertex out of range");
312 fBGEventHandler = fV0Reader->GetBGHandler();
313 fESDpid = fV0Reader->GetESDpid();
314 fESDEvent = fV0Reader->GetESDEvent();
315 if(fDoMC && MCEvent())
317 fStack= MCEvent()->Stack();
318 fGCMCEvent=MCEvent();
321 // Read the magnetic field sign from ESD
322 if ( fReadMagFieldSign == kTRUE )
324 fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
327 // Process MC information
334 while(fV0Reader->NextV0()){}; //SelectGammas
335 fV0Reader->ResetV0IndexNumber();
338 CreateListOfDalitzPairCandidates();
339 ProcessGammaElectronsForDalitzAnalysis();
343 fV0Reader->UpdateEventByEventData();
347 PostData( 1, fOutputContainer );
351 void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
354 if( fDebug ) AliInfo("Not to do anything in Terminate");
357 //-----------------------------------------------------------------------------------------------
358 void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
361 // set user ITSsa track cuts
363 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
367 fITSsaTrackCuts = esdCuts;
372 fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
374 fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
375 fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
377 fITSsaTrackCuts->SetMinNClustersITS(2);
378 fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
379 fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
381 fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
382 fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
384 fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
388 //-----------------------------------------------------------------------------------------------
389 void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
392 // set user global track cuts
394 if( fESDtrackCuts ) delete fESDtrackCuts;
398 fESDtrackCuts = esdCuts;
403 fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
405 fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
406 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
408 fESDtrackCuts->SetMinNClustersITS(2);
409 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
410 fESDtrackCuts->SetRequireITSRefit(kTRUE);
412 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
413 fESDtrackCuts->SetMaxNsigmaToVertex(3);
415 fESDtrackCuts->SetMinNClustersTPC(80);
416 fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
417 fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
418 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
422 //-----------------------------------------------------------------------------------------------
423 void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
428 if( fESDpidCuts ) delete fESDpidCuts;
431 fESDpidCuts = esdPIDCuts;
435 fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
436 // fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
437 fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
441 //-----------------------------------------------------------------------------------------------
442 void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
445 // Process generation
447 if( fDebug ) AliInfo("=> ProcessMCData");
449 fHistograms->FillTable("Table_Generation", 0); //number of events
451 for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
453 TParticle* iParticle = fStack->Particle( i );
454 if( !iParticle ) continue;
456 if ( i >= fStack->GetNprimary() ) continue; // is primary?
457 if ( iParticle->GetPdgCode() != ::kPi0 ) continue; // is Pi0?
459 if( iParticle->GetNDaughters() == 2 &&
460 fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
461 fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
463 fHistograms->FillTable("Table_Generation", 1); // pi0 -> gg
466 if ( iParticle->GetNDaughters() != 3 ) continue; // Num == 3 (e+,e-,gamma)
468 // Check for Pi0 Dalitz decay
469 TParticle* eposPi0 = 0;
470 TParticle* enegPi0 = 0;
471 TParticle* gammaPi0 = 0;
473 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
475 switch(fStack->Particle(idxPi0)->GetPdgCode())
478 eposPi0 = fStack->Particle(idxPi0);
481 enegPi0 = fStack->Particle(idxPi0);
484 gammaPi0 = fStack->Particle(idxPi0);
489 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
491 // found a Pi0 Dalitz decay
493 fHistograms->FillTable("Table_Generation", 2);
494 fHistograms->FillHistogram("MC_Pi0Dalitz_P", iParticle->P());
495 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt", iParticle->Pt());
496 fHistograms->FillHistogram("MC_Pi0Dalitz_Eta", iParticle->Eta());
497 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
498 fHistograms->FillHistogram("MC_EposDalitz_Pt", eposPi0->Pt());
499 fHistograms->FillHistogram("MC_EposDalitz_Eta", eposPi0->Eta());
500 fHistograms->FillHistogram("MC_EnegDalitz_Pt", enegPi0->Pt());
501 fHistograms->FillHistogram("MC_EnegDalitz_Eta", enegPi0->Eta());
502 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Pt", gammaPi0->Pt());
503 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Eta", gammaPi0->Eta());
505 // Angle between the gamma and the plane e+e-
506 TVector3 ePosMom( eposPi0->Px(), eposPi0->Py(), eposPi0->Pz() );
507 TVector3 eNegMom( enegPi0->Px(), enegPi0->Py(), enegPi0->Pz() );
508 TVector3 gamMom( gammaPi0->Px(), gammaPi0->Py() , gammaPi0->Pz() );
509 TVector3 planeEposEneg = eNegMom.Cross( ePosMom );
510 Double_t anglePlaneGamma = planeEposEneg.Angle(gamMom);
512 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
514 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle", anglePlaneGamma);
515 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_P", anglePlaneGamma, gammaPi0->P());
516 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_Pt", anglePlaneGamma, gammaPi0->Pt());
519 // check for gamma conversion
520 Bool_t daugGammaElectron = kFALSE;
521 Bool_t daugGammaPositron = kFALSE; // acceptance
522 Bool_t daugGammaElectronAll = kFALSE;
523 Bool_t daugGammaPositronAll = kFALSE;
525 // is the gamma converted? -> has 2 daughter e+e-
526 // are e+ e- from gamma in the acceptance for the V0s
527 if( gammaPi0->GetNDaughters() >= 2 )
529 for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
531 TParticle* tmpDaughter = fStack->Particle(tIndex);
533 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
534 if( tmpDaughter->GetPdgCode() == ::kElectron )
536 daugGammaElectronAll = kTRUE;
538 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
539 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
540 (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
542 daugGammaElectron = kTRUE;
545 else if( tmpDaughter->GetPdgCode() == ::kPositron )
547 daugGammaPositronAll = kTRUE;
548 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
549 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
550 (tmpDaughter->R() < fV0Reader->GetMaxRCut() ) )
552 daugGammaPositron = kTRUE;
559 if( daugGammaElectronAll && daugGammaPositronAll )
561 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
562 fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
565 Float_t etaMin, etaMax;
566 fESDtrackCuts->GetEtaRange( etaMin, etaMax );
568 // e+e- pair within acceptance
569 if ( TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
571 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt", eposPi0->Pt());
572 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta", eposPi0->Eta());
573 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt", enegPi0->Pt());
574 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta", enegPi0->Eta());
575 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_EposPt_vs_EnegPt", eposPi0->Pt(), enegPi0->Pt());
578 // Pi0 (e+e-gamma) within acceptance
579 //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
581 if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
582 TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
584 fHistograms->FillTable("Table_Generation",3); //
585 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt",iParticle->Pt());
586 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Eta",iParticle->Eta());
587 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Radius",iParticle->R());
588 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Pt",gammaPi0->Pt());
589 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Eta",gammaPi0->Eta());
590 fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Pt",eposPi0->Pt());
591 fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Eta",eposPi0->Eta());
592 fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Pt",enegPi0->Pt());
593 fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Eta",enegPi0->Eta());
594 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
595 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
597 // Pi0 within acceptance with gamma converted
599 if ( daugGammaElectron && daugGammaPositron )
601 fHistograms->FillTable("Table_Generation",4); //
603 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
604 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
605 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Pt",eposPi0->Pt());
606 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Eta",eposPi0->Eta());
607 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Pt",enegPi0->Pt());
608 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Eta",enegPi0->Eta());
609 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Pt",gammaPi0->Pt());
610 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Eta",gammaPi0->Eta());
611 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle",anglePlaneGamma);
612 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle_vs_Pt",anglePlaneGamma,gammaPi0->Pt());
613 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
614 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
615 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
621 //-----------------------------------------------------------------------------------------------
622 void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
625 // Dalitz pair candidates
627 if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
629 fEposCandidateIndex.clear();
630 fEnegCandidateIndex.clear();
632 fHistograms->FillTable("Table_Cuts", 0);
634 for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
636 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
637 if ( !iTrack ) continue;
642 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
644 TVector3 iMom(p[0],p[1],p[2]);
647 // Check track cuts and find track type
650 Bool_t isTrackAccepted = 0;
651 Int_t trackType = -1;
652 switch(fTrkSelectionCriteria)
655 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
656 trackType = kITSsaTrack;
660 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
661 trackType = kGlobalTrack;
664 case kITSsaGlobalTrack:
665 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
667 isTrackAccepted = kTRUE;
668 if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
669 else trackType = kGlobalTrack;
674 if(!isTrackAccepted) continue;
685 pid = GetBayesPid(iTrack,trackType);
689 pid = GetNSigmaPid(iTrack,trackType);
694 pidMC = GetMonteCarloPid(iTrack);
696 Int_t iLabel = TMath::Abs(iTrack->GetLabel());
697 TParticle* iParticle = fStack->Particle(iLabel);
698 FillPidTable(iParticle, pid);
701 // ITS standalone tracks
702 if( trackType == kITSsaTrack)
704 Double_t mom = iTrack->GetP();
705 Double_t signal = iTrack->GetITSsignal();
707 fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
709 if( pid == AliPID::kElectron )
712 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
713 if(fDoMC && pid == pidMC)
715 fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
719 if( fDoMC && pidMC == AliPID::kElectron)
721 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
725 else // global tracks
727 const AliExternalTrackParam *in = iTrack->GetInnerParam();
728 Double_t mom = in->GetP();
729 Double_t signal = iTrack->GetTPCsignal();
731 fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
733 if( fDoMC && pidMC == AliPID::kElectron )
735 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
738 if( pid == AliPID::kElectron )
740 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
741 if(fDoMC && pid == pidMC)
743 fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
748 if( AliPID::kElectron != pid) continue;
750 // electron track candidates from here
752 if( iTrack->GetSign() > 0 )
754 fEposCandidateIndex.push_back(i);
758 fEnegCandidateIndex.push_back(i);
763 GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
767 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
768 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
772 if(fUseTrackIndexCut) // remove repeated tracks
774 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
778 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
779 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
784 if(fUsePsiPairCut) // remove electrons from gamma conversions
786 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
790 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
791 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
798 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
802 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
803 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
810 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
814 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
815 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
821 //-----------------------------------------------------------------------------------------------
822 void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
825 // Process gamma and electrons for pi0 Dalitz decay
827 if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
829 fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
831 TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
833 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
835 AliKFParticle* epos = (AliKFParticle*) ePosCandidates->At(i);
836 fHistograms->FillHistogram("ESD_EposCandidates_Pt", epos->GetPt());
837 fHistograms->FillHistogram("ESD_EposCandidates_Eta", epos->GetEta());
838 fHistograms->FillTable( "Table_Reconstruction", 1);
841 TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
843 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
845 AliKFParticle* eneg = (AliKFParticle*) eNegCandidates->At(i);
846 fHistograms->FillHistogram("ESD_EnegCandidates_Pt", eneg->GetPt());
847 fHistograms->FillHistogram("ESD_EnegCandidates_Eta", eneg->GetEta());
848 fHistograms->FillTable( "Table_Reconstruction", 2);
851 TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
852 for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
854 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
856 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
857 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Mass", dalitz->M());
861 for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
863 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
864 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
865 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
868 // psi pair for all candidates
870 FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
872 // Angle epos,eneg gamma
873 FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
874 FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
876 TClonesArray* pi0Candidates = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaCandidates);
877 for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
879 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
881 fHistograms->FillHistogram("ESD_Pi0_P", pi0->P());
882 fHistograms->FillHistogram("ESD_Pi0_Pt", pi0->Pt());
883 fHistograms->FillHistogram("ESD_Pi0_Eta", pi0->Eta());
884 fHistograms->FillHistogram("ESD_Pi0_Y", pi0->Rapidity());
885 fHistograms->FillHistogram("ESD_Pi0_Phi", pi0->Phi());
886 fHistograms->FillHistogram("ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
887 fHistograms->FillHistogram("ESD_Pi0_Mass", pi0->M());
888 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
889 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Y", pi0->Rapidity(),pi0->M());
890 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
893 delete dalitzPairCandidates;
894 delete pi0Candidates;
898 // 1) e+e- with with gammas used in the signal
899 TClonesArray* pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaPool);
901 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
903 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
904 fHistograms->FillHistogram("ESD_BKG_PrevGamma", pi0->M());
905 fHistograms->FillHistogram("ESD_BKG_PrevGamma_vs_Pt", pi0->Pt(), pi0->M());
908 if(ePosCandidates->GetEntriesFast() > 0 &&
909 eNegCandidates->GetEntriesFast() > 0 &&
910 fGammaCandidates->GetEntriesFast() > 0)
912 UpdateGammaPool(fGammaCandidates);
917 // 2) e+e- with gammas from a pool of events
918 TClonesArray* gammaBGHandler = GammasFromBGHandler();
919 pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, gammaBGHandler);
921 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
923 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
924 fHistograms->FillHistogram("ESD_BKG_BGHandler", pi0->M());
925 fHistograms->FillHistogram("ESD_BKG_BGHandler_vs_Pt", pi0->Pt(), pi0->M());
930 // 3) e+ with e-, gamma from a pool of events
931 TClonesArray* elecBGHandler = ElectronFromBGHandler();
932 pi0Bkg = FindPi0Dalitz(ePosCandidates, elecBGHandler, gammaBGHandler);
934 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
936 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
937 fHistograms->FillHistogram("ESD_BKG_Electron", pi0->M());
938 fHistograms->FillHistogram("ESD_BKG_Electron_vs_Pt", pi0->Pt(), pi0->M());
941 if(eNegCandidates->GetEntriesFast() > 0)
943 UpdateElectronPool(eNegCandidates);
946 delete gammaBGHandler;
947 delete elecBGHandler;
952 delete ePosCandidates;
953 delete eNegCandidates;
957 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
958 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
960 AliKFParticle* epos = (AliKFParticle*) ePosPi0Dalitz->At(i);
961 fHistograms->FillHistogram("MC_ESD_EposDalitz_Pt", epos->GetPt());
962 fHistograms->FillHistogram("MC_ESD_EposDalitz_Eta", epos->GetEta());
963 fHistograms->FillTable( "Table_Reconstruction", 3);
966 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
967 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
969 AliKFParticle* eneg = (AliKFParticle*) eNegPi0Dalitz->At(i);
970 fHistograms->FillHistogram("MC_ESD_EnegDalitz_Pt", eneg->GetPt());
971 fHistograms->FillHistogram("MC_ESD_EnegDalitz_Eta", eneg->GetEta());
972 fHistograms->FillTable( "Table_Reconstruction", 4);
975 TClonesArray* dalitzPair = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex);
976 for(Int_t i=0; i < dalitzPair->GetEntriesFast(); ++i)
978 TLorentzVector* dalitz = (TLorentzVector*) dalitzPair->At(i);
979 fHistograms->FillHistogram("MC_ESD_DalitzPair_Pt", dalitz->Pt());
980 fHistograms->FillHistogram("MC_ESD_DalitzPair_Mass", dalitz->M());
981 fHistograms->FillHistogram( "Table_Reconstruction", 5 );
984 // psi pair for dalitz pairs
986 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_DalitzPair_PsiPair_vs_DPhi");
988 delete ePosPi0Dalitz;
989 delete eNegPi0Dalitz;
993 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
994 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
996 AliKFParticle* iGamma = (AliKFParticle*) gamma->At(i);
997 fHistograms->FillHistogram("MC_ESD_Gamma_Pt", iGamma->GetPt());
998 fHistograms->FillHistogram("MC_ESD_Gamma_Eta", iGamma->GetEta());
1003 // gamma from pi0 dalitz
1004 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1005 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
1007 AliKFParticle* iGamma = (AliKFParticle*) gammaPi0Dalitz->At(i);
1008 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Pt", iGamma->GetPt());
1009 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Eta", iGamma->GetEta());
1010 fHistograms->FillTable( "Table_Reconstruction", 6);
1013 delete gammaPi0Dalitz;
1015 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1017 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
1019 TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
1021 fHistograms->FillHistogram("MC_ESD_Pi0_P", pi0->P());
1022 fHistograms->FillHistogram("MC_ESD_Pi0_Pt", pi0->Pt());
1023 fHistograms->FillHistogram("MC_ESD_Pi0_Eta", pi0->Eta());
1024 fHistograms->FillHistogram("MC_ESD_Pi0_Y", pi0->Rapidity());
1025 fHistograms->FillHistogram("MC_ESD_Pi0_Phi", pi0->Phi());
1026 fHistograms->FillHistogram("MC_ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
1027 fHistograms->FillHistogram("MC_ESD_Pi0_Mass", pi0->M());
1028 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
1029 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Y", pi0->Rapidity(),pi0->M());
1030 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
1031 fHistograms->FillHistogram( "Table_Reconstruction", 7);
1035 // psi pair for electrons from gamma conversions assuming they came from main vertex
1036 // if(fUsePsiPairCut)
1037 for(UInt_t i=0; i < fEposCandidateIndex.size(); ++i)
1039 AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
1040 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1042 for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
1044 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
1045 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1047 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
1049 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1050 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1052 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
1055 // FIXME: eta -> e+e-gamma
1059 //--------------------------------------------------------------------------
1060 Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
1065 const double kEPSILON=1.e-16;
1067 if(p->Energy() - TMath::Abs(p->Pz()) < kEPSILON )
1071 return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
1074 //--------------------------------------------------------------------------
1075 void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1078 // Fill histogram with psipair(pos,neg)
1080 for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1082 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1083 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1085 AliKFParticle* negKF = (AliKFParticle*) neg->At(j);
1086 Double_t psiPair = GetPsiPair(posKF, negKF);
1087 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negKF->GetPhi() - posKF->GetPhi());
1088 fHistograms->FillHistogram(hName, deltaPhi, psiPair);
1093 //--------------------------------------------------------------------------
1094 void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1097 // Fill histogram with angle(x,y)
1099 for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1101 AliKFParticle* xKF = (AliKFParticle*) x->At(i);
1102 TVector3 xMom(xKF->Px(),xKF->Py(),xKF->Pz());
1103 for( Int_t j=0; j < y->GetEntriesFast(); ++j )
1105 AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1106 TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1107 fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1112 //--------------------------------------------------------------------------
1113 void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1116 // Fill table with pid info
1119 switch(TMath::Abs(p->GetPdgCode()))
1121 case ::kElectron: iGen=0; break;
1122 case ::kMuonMinus: iGen=1; break;
1123 case ::kPiPlus: iGen=2; break;
1124 case ::kKPlus: iGen=3; break;
1125 case ::kProton: iGen=4; break;
1130 if(pid > -1 && pid < 5) jRec = pid;
1132 if ((iGen > -1) && (jRec > -1))
1134 fHistograms->FillTable("Table_PID", iGen, jRec);
1136 fHistograms->FillTable("Table_PID", iGen, 5);
1137 fHistograms->FillTable("Table_PID", 5, jRec);
1141 //--------------------------------------------------------------------------
1142 void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1145 // Make a copy of gamma candidates from V0reader
1150 if(gamma) delete gamma;
1152 TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1154 gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1155 gamma->SetOwner(kTRUE);
1158 for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1160 AliKFParticle* gamKF = (AliKFParticle*)gammaV0->At(i);
1161 new ((*gamma)[i]) AliKFParticle(*gamKF);
1162 posIndex.push_back(fV0Reader->GetPindex(i));
1163 negIndex.push_back(fV0Reader->GetNindex(i));
1167 //--------------------------------------------------------------------------
1168 TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1171 // Convert track index vector to AliKFParticle array
1173 TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1174 indexKF->SetOwner(kTRUE);
1176 for(UInt_t i = 0; i < index.size(); ++i)
1178 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1179 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1185 //--------------------------------------------------------------------------
1186 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1189 // Find true electrons from pi0 Dalitz decay candidates with MC
1191 TClonesArray* elec = new TClonesArray("AliKFParticle");
1192 elec->SetOwner(kTRUE);
1194 for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1196 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1197 Int_t trackLabel = TMath::Abs(track->GetLabel());
1199 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1200 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1202 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1208 //--------------------------------------------------------------------------
1209 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1212 // Find true gammas from pi0 Dalitz decay candidates with MC
1214 TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1215 gammaPi0->SetOwner(kTRUE);
1217 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1219 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1221 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1222 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1224 if( !HaveSameMother(labelv1,labelv2) ) continue;
1226 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1228 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1230 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1232 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1238 //--------------------------------------------------------------------------
1239 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1242 // Find true gammas from gamma candidates with MC
1244 TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1245 gammaConv->SetOwner(kTRUE);
1247 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1249 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1251 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1252 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1254 if( !HaveSameMother(labelv1,labelv2) ) continue;
1256 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1258 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1260 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1266 //--------------------------------------------------------------
1267 void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1270 // Remove repeated electron candidate tracks
1271 // according to the gamma candidate array
1273 vector<Bool_t> posTag(pos.size(),kTRUE);
1274 vector<Bool_t> negTag(neg.size(),kTRUE);
1276 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1278 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1279 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1281 for( UInt_t j=0; j < pos.size(); ++j )
1283 if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1286 for( UInt_t j=0; j < neg.size(); ++j )
1288 if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1292 CleanArray(pos, posTag);
1293 CleanArray(neg, negTag);
1296 //--------------------------------------------------------------------------
1297 void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1300 // Remove electron candidates from gamma conversions
1301 // according to the Psi pair angle
1303 vector<Bool_t> posTag(pos.size(), kTRUE);
1304 vector<Bool_t> negTag(neg.size(), kTRUE);
1306 for( UInt_t i=0; i < pos.size(); ++i )
1308 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1310 for( UInt_t j=0; j < neg.size(); ++j )
1312 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1314 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1315 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1317 if(IsFromGammaConversion( psiPair, deltaPhi ))
1325 CleanArray(pos, posTag);
1326 CleanArray(neg, negTag);
1329 //-----------------------------------------------------------------------------------
1330 void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1333 // Remove electron candidates pairs
1334 // which have mass not in the range (fMassCutMin,fMassCutMax)
1336 vector<Bool_t> posTag(pos.size(), kTRUE);
1337 vector<Bool_t> negTag(neg.size(), kTRUE);
1339 for( UInt_t i=0; i < pos.size(); ++i )
1341 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1343 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1344 TLorentzVector posLV;
1345 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1347 for( UInt_t j=0; j < neg.size(); ++j )
1349 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1351 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1352 TLorentzVector negLV;
1353 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1355 TLorentzVector posnegLV = posLV + negLV;
1357 if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1365 CleanArray(pos, posTag);
1366 CleanArray(neg, negTag);
1369 //-----------------------------------------------------------------------------------------------
1370 void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1373 // Clean the x array according to the tag parameter
1377 for(UInt_t i=0; i< x.size(); ++i)
1379 if(tag[i]) tmp.push_back(x[i]);
1385 //--------------------------------------------------------------------------
1386 void AliAnalysisTaskGammaConvDalitz::AngleEposEnegGammaCut( const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* candidates, TClonesArray*& gamma, vector<Int_t>& posGamIdx, vector<Int_t>& negGamIdx)
1389 // Remove gamma candidates according to
1390 // the angle between the plane e+,e- and the gamma
1392 vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1394 for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1396 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1397 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1398 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1400 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1402 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1403 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1404 TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1406 // normal vector to x+y- plane
1407 TVector3 planePosNeg = xMom.Cross(yMom);
1408 for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1410 AliKFParticle* gamKF = (AliKFParticle*)candidates->At(i);
1411 TVector3 gamMom(gamKF->Px(),gamKF->Py(),gamKF->Pz());
1412 if (planePosNeg.Angle(gamMom) < 1. || planePosNeg.Angle(gamMom) > 2.)
1414 gammaTag[i] = kFALSE;
1420 // Rebuild gamma candidates array
1422 if(gamma) delete gamma;
1423 gamma = new TClonesArray("AliKFParticle");
1424 gamma->SetOwner(kTRUE);
1429 for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1431 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1434 new ((*gamma)[j++]) AliKFParticle(*iGamma);
1435 posGamIdx.push_back(fV0Reader->GetPindex(i));
1436 negGamIdx.push_back(fV0Reader->GetNindex(i));
1441 //--------------------------------------------------------------------------
1442 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1445 // Find Dalitz pair candidates
1447 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1448 dalitz->SetOwner(kTRUE);
1450 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1452 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1454 TLorentzVector posLV;
1455 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1457 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1459 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1461 TLorentzVector negLV;
1462 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1466 AliKFParticle posNeg( *posKF, *negKF);
1468 TLorentzVector posNegLV;
1469 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1470 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1474 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1482 //--------------------------------------------------------------------------
1483 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma)
1486 // Find pi0 Dalitz decay candidates
1488 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1489 pi0->SetOwner(kTRUE);
1491 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1493 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1495 TLorentzVector posLV;
1496 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1498 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1500 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1502 TLorentzVector negLV;
1503 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1505 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1507 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1511 AliKFParticle posNegGam( *posKF, *negKF, *gamKF );
1512 TLorentzVector posNegGamLV;
1513 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1514 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1518 TLorentzVector gamLV;
1519 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1521 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1530 //--------------------------------------------------------------------------
1531 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1534 // Find true Dalitz pairs from Dalitz pair candidats with MC
1536 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1537 dalitz->SetOwner(kTRUE);
1539 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1541 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1542 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1543 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1545 TLorentzVector posLV;
1546 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1548 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1550 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1552 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1553 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1555 if(!IsDalitzPair(posLabel,negLabel)) continue;
1559 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1560 AliKFParticle posNeg( posKF, negKF);
1562 TLorentzVector posNegLV;
1563 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1565 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1567 else // TLorentzVector
1569 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1571 TLorentzVector negLV;
1572 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1574 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1582 //--------------------------------------------------------------------------
1583 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam)
1586 // Find true pi0 Dalitz decay from pi0 candidates with MC
1588 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1589 pi0->SetOwner(kTRUE);
1591 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1593 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1594 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1595 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1597 TLorentzVector posLV;
1598 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1600 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1602 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1604 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1605 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1607 if(!IsDalitzPair(posLabel,negLabel)) continue;
1609 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1611 TLorentzVector negLV;
1612 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1614 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1616 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1618 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1620 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1621 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1623 if( !HaveSameMother(labelv1,labelv2) ) continue;
1625 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1627 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1628 if( !HaveSameMother(labelGamma, posLabel) ) continue;
1632 AliKFParticle posNegGam( posKF, negKF, *gamKF );
1633 TLorentzVector posNegGamLV;
1634 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1635 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1637 else // TLorentzVector
1639 TLorentzVector gamLV;
1640 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1642 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1651 //-----------------------------------------------------------------------------------------------
1652 void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
1655 // Update gamma event pool for background computation
1657 if( fDebug ) AliInfo("=> UpdateGammaPool");
1660 for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
1662 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
1663 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
1665 if(fGamPoolPos == fPoolMaxSize)
1672 void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
1675 // Update electron event pool for background computation
1677 Int_t multiplicity = fV0Reader->CountESDTracks();
1680 fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
1683 //-----------------------------------------------------------------------------------------------
1684 TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
1687 // Gamma copy from events with same multiplicity and Z
1689 if( fDebug ) AliInfo("=> GammasFromBGHandler");
1691 Int_t zbin = fBGEventHandler->GetZBinIndex(fV0Reader->GetVertexZ());
1692 Int_t mbin = fBGEventHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1696 TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
1697 gammaPool->SetOwner(kTRUE);
1699 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1701 AliGammaConversionKFVector* gammaV0s = fBGEventHandler->GetBGGoodV0s(zbin,mbin,iEventBG);
1702 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
1704 new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
1711 //-----------------------------------------------------------------------------------------------
1712 TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
1715 // Electron copy from events with same multiplicity and Z
1717 if( fDebug ) AliInfo("=> ElectronFromBGHandler");
1719 TClonesArray* electronPool = new TClonesArray("AliKFParticle");
1720 electronPool->SetOwner(kTRUE);
1722 Int_t multiplicity = fV0Reader->CountESDTracks();
1727 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1729 AliGammaConversionKFVector* electronNeg = fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
1730 for (UInt_t i = 0; i < electronNeg->size(); ++i )
1732 new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
1736 return electronPool;
1739 //-----------------------------------------------------------------------------------------------
1740 Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
1743 // Get track pid according to MC
1745 Int_t label = TMath::Abs(t->GetLabel());
1746 Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
1750 case ::kElectron: return AliPID::kElectron;
1751 case ::kMuonMinus: return AliPID::kMuon;
1752 case ::kPiPlus: return AliPID::kPion;
1753 case ::kKPlus: return AliPID::kKaon;
1754 case ::kProton: return AliPID::kProton;
1760 //-----------------------------------------------------------------------------------------------
1762 // NOTE prior should be estimated from data
1763 // NOTE: move to config
1765 Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
1768 // Get track pid according to Bayes' formula
1770 double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
1771 Double_t detectoProb[AliPID::kSPECIES];
1773 if( trackType == kITSsaTrack ) // ITS standalone pid
1775 t->GetITSpid( detectoProb );
1777 else // global track
1779 t->GetESDpid( detectoProb );
1782 AliPID bayesPID( detectoProb );
1783 return bayesPID.GetMostProbable( priors );
1786 //-----------------------------------------------------------------------------------------------
1787 Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
1790 // Get track pid according to a n-sigma cut around ITS and/or TPC signals
1792 if( trackType == kITSsaTrack) // ITS standalone tracks
1794 Double_t mom = t->GetP();
1797 // ITS Number of sigmas (FIXME: add new fESDpidCuts)
1798 // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
1799 Double_t nElecSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kElectron );
1800 Double_t nPionSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kPion );
1802 if( nElecSigma < 4. && nElecSigma > -3. && mom < .2 && nPionSigma < -1.5 )
1804 return AliPID::kElectron;
1807 else // global track
1809 Double_t nElecSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kElectron );
1810 Double_t nPionSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kPion );
1811 Double_t nKaonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kKaon ));
1812 Double_t nProtonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kProton));
1813 if( nElecSigma > fNSigmaBelowElecTPCbethe && nElecSigma < fNSigmaAboveElecTPCbethe &&
1814 nPionSigma > fNSigmaAbovePionTPCbethe && //NOTE mom > 0.5
1815 nKaonSigma > fNSigmaAboveKaonTPCbethe &&
1816 nProtonSigma > fNSigmaAboveProtonTPCbethe )
1818 return AliPID::kElectron;
1820 // NOTE: add other particle types
1826 //-----------------------------------------------------------------------------------------------
1827 Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel ) const
1830 // Returns true if the two particles is a Dalitz pair
1832 if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
1834 TParticle* pos = fStack->Particle( posLabel );
1835 TParticle* neg = fStack->Particle( negLabel );
1837 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1838 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1840 //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
1841 //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
1843 Int_t motherLabel = pos->GetMother(0);
1844 if( motherLabel < 0 ) return kFALSE;
1846 TParticle* mother = fStack->Particle( motherLabel );
1848 if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1849 if( mother->GetNDaughters() != 3) return kFALSE;
1850 // NOTE: one of them must be a photon
1855 //-----------------------------------------------------------------------------------------------
1856 Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
1859 // Returns true if the particle comes from Pi0 -> e+ e- gamma
1861 Bool_t ePlusFlag = kFALSE;
1862 Bool_t eMinusFlag = kFALSE;
1863 Bool_t gammaFlag = kFALSE;
1865 Int_t motherLabel = fStack->Particle( label )->GetMother(0);
1867 if( motherLabel < 0 ) return kFALSE;
1869 TParticle* mother = fStack->Particle( motherLabel );
1871 if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1873 if ( mother->GetNDaughters() != 3 ) return kFALSE;
1875 for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
1877 switch( fStack->Particle(idx)->GetPdgCode())
1891 return ( ePlusFlag && eMinusFlag && gammaFlag );
1894 //--------------------------------------------------------------------------
1895 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
1898 // Returns true if it is a gamma conversion according to psi pair value
1900 return ( (deltaPhi > fDeltaPhiCutMin && deltaPhi < fDeltaPhiCutMax) &&
1901 TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
1904 //--------------------------------------------------------------------------
1905 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
1908 // Returns true if it is a gamma conversion according to MC
1910 if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
1912 TParticle* pos = fStack->Particle( posLabel );
1913 TParticle* neg = fStack->Particle( negLabel );
1915 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1916 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1918 if( pos->GetUniqueID() != kPPair ) return kFALSE;
1920 Int_t motherLabel = pos->GetMother(0);
1921 if( motherLabel < 0 ) return kFALSE;
1923 TParticle* mother = fStack->Particle( motherLabel );
1925 if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
1930 //-----------------------------------------------------------------------------------------------
1931 Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
1934 // Returns true if the two particle have the same mother
1936 if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
1937 return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
1940 //-----------------------------------------------------------------------------------------------
1941 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
1944 // This angle is a measure for the contribution of the opening in polar
1945 // direction Δ0 to the opening angle ξ Pair
1947 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
1948 // Master Thesis. Thorsten Dahms. 2005
1949 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
1953 if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
1954 if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
1956 TVector3 posDaughter;
1957 TVector3 negDaughter;
1959 posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1960 negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
1962 Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1963 Double_t openingAngle = posDaughter.Angle( negDaughter ); //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
1964 if( openingAngle < 1e-20 ) return 0.;
1965 Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
1970 //-----------------------------------------------------------------------------------------------
1971 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const AliKFParticle* xPos, const AliKFParticle* yNeg ) const
1974 // Get psi pair value
1976 TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
1977 TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
1979 Double_t deltaTheta = neg.Theta() - pos.Theta();
1980 Double_t openingAngle = pos.Angle( neg );
1982 if( openingAngle < 1e-20 ) return 0.;
1984 return TMath::ASin( deltaTheta/openingAngle );
1987 //-----------------------------------------------------------------------------------------------
1988 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
1991 // Get psi pair value
1993 Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
1994 Double_t openingAngle = xPos->Angle( yNeg->Vect() );
1996 if( openingAngle < 1e-20 ) return 0.;
1998 return TMath::ASin( deltaTheta/openingAngle );;
2001 // tmp NOTE: Should go to AliV0Reader
2002 //-----------------------------------------------------------------------------------------------
2004 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiAngleV0(
2005 Double_t radiusVO, // radius at XY of VO vertex
2006 const AliExternalTrackParam* trackPos, // pos. track parm. at V0 vertex
2007 const AliExternalTrackParam* trackNeg // neg. track parm. at V0 vertex
2010 // This angle is a measure for the contribution of the opening in polar
2011 // direction Δ0 to the opening angle ξ Pair
2013 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2014 // Master Thesis. Thorsten Dahms. 2005
2015 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2017 const Double_t kForceDeltaPhi = 1.2;
2018 static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
2020 Double_t psiAngle = 0.; // Default value
2022 static Double_t MagnField = fESDEvent->GetMagneticField();
2023 static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
2025 if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
2027 // compute propagation radius for a fixed angle
2028 Double_t Rpos = radiusVO + trackPos->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2029 Double_t Rneg = radiusVO + trackNeg->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2033 if( trackPos->GetPxPyPzAt( Rpos, MagnField, MomPos ) == 0 ) trackPos->GetPxPyPz( MomPos );
2034 if( trackNeg->GetPxPyPzAt( Rneg, MagnField, MomNeg ) == 0 ) trackNeg->GetPxPyPz( MomNeg );
2036 TVector3 PosDaughter;
2037 TVector3 NegDaughter;
2038 PosDaughter.SetXYZ( MomPos[0], MomPos[1], MomPos[2] );
2039 NegDaughter.SetXYZ( MomNeg[0], MomNeg[1], MomNeg[2] );
2041 Double_t deltaTheta = NegDaughter.Theta() - PosDaughter.Theta();
2042 Double_t chiPar = PosDaughter.Angle( NegDaughter ); //TMath::ACos( PosDaughter.Dot(NegDaughter) / (NegDaughter.Mag()*PosDaughter.Mag()) );
2043 if( chiPar < 1e-20 ) return psiAngle;
2045 psiAngle = TMath::ASin( deltaTheta / chiPar );