]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/GammaConv/AliAnalysisTaskGammaConvDalitz.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix Coverity
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConvDalitz.cxx
CommitLineData
7f3c7cc6 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Pedro González, Pedro Ladrón de Guevara, Ernesto López Torres, *
5 * Eulogio Serradilla *
6 * Version 2 *
7 * *
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 **************************************************************************/
16
17// Analysis task for pi0->e+e-gamma (Dalitz decay)
18
19#include <vector>
20
21#include "TParticle.h"
22#include "TPDGCode.h"
23#include "TMCProcess.h"
24#include "TDatabasePDG.h"
25#include "TList.h"
26#include "TChain.h"
27#include "TDirectory.h"
28
29#include "AliStack.h"
30#include "AliAnalysisManager.h"
31#include "AliESDInputHandler.h"
32#include "AliESDtrack.h"
33#include "AliMCEvent.h"
34#include "AliStack.h"
35#include "AliMCEventHandler.h"
36#include "AliPID.h"
37#include "AliLog.h"
38#include "AliESDtrackCuts.h"
39#include "AliESDpidCuts.h"
40#include "AliMCEvent.h"
41#include "AliESDv0.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"
49
50#include "AliAnalysisTaskGammaConvDalitz.h"
51#include "TH1.h"
52
53ClassImp( AliAnalysisTaskGammaConvDalitz )
54
55//-----------------------------------------------------------------------------------------------
56AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz():
57 AliAnalysisTaskSE(),
58 fStack(0),
7811d09e 59 fGCMCEvent(0),
7f3c7cc6 60 fESDEvent(0),
61 fEposCandidateIndex(),
62 fEnegCandidateIndex(),
63 fGammaCandidatePosIndex(),
64 fGammaCandidateNegIndex(),
65 fGammaCandidates(0),
66 fGammaPool(0),
67 fPoolMaxSize(10),
68 fGamPoolPos(0),
69 fBGEventHandler(0),
70 fOutputContainer(0),
71 fMCTruth(0),
72 fV0Reader(0),
73 fESDpid(0),
74 fESDtrackCuts(0),
75 fITSsaTrackCuts(0),
76 fESDpidCuts(0),
77 fHistograms(0),
78 fStandalone(kFALSE),
79 fDoMC(kFALSE),
80 fComputeBkg(kTRUE),
81 fUseBayesPID(kFALSE),
82 fUseTrackIndexCut(kTRUE),
83 fUsePsiPairCut(kTRUE),
84 fUseMassCut(kFALSE),
85 fUseGammaCut(kFALSE),
86 fReadMagFieldSign(kTRUE),
87 fUseAliKF(kFALSE),
88 fMagFieldSign(1),
b9295e41 89 fkElectronMass(0.00051099891),
7f3c7cc6 90 fPsiPairCut(0.45),
91 fDeltaPhiCutMin(0.),
92 fDeltaPhiCutMax(0.12),
93 fMassCutMin(0.),
94 fMassCutMax(0.1),
95 fNSigmaBelowElecTPCbethe(-2.),
96 fNSigmaAboveElecTPCbethe(3.),
97 fNSigmaAbovePionTPCbethe(3.),
98 fNSigmaAboveKaonTPCbethe(3.),
99 fNSigmaAboveProtonTPCbethe(3.),
100 fTrkSelectionCriteria(kGlobalTrack)
101{
7811d09e 102//
103// Default constructor
104//
7f3c7cc6 105 AdoptITSsaTrackCuts();
106 AdoptESDtrackCuts();
107 AdoptESDpidCuts();
108
109 fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
110 fGammaPool->SetOwner(kTRUE);
111}
112
113//-----------------------------------------------------------------------------------------------
114AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz( const char* name ):
115 AliAnalysisTaskSE( name ),
116 fStack(0),
7811d09e 117 fGCMCEvent(0),
7f3c7cc6 118 fESDEvent(0),
119 fEposCandidateIndex(),
120 fEnegCandidateIndex(),
121 fGammaCandidatePosIndex(),
122 fGammaCandidateNegIndex(),
123 fGammaCandidates(0),
124 fGammaPool(0),
125 fPoolMaxSize(10),
126 fGamPoolPos(0),
127 fBGEventHandler(0),
128 fOutputContainer(0),
129 fMCTruth(0),
130 fV0Reader(0),
131 fESDpid(0),
132 fESDtrackCuts(0),
133 fITSsaTrackCuts(0),
134 fESDpidCuts(0),
135 fHistograms(0),
136 fStandalone(kFALSE),
137 fDoMC(kFALSE),
138 fComputeBkg(kTRUE),
139 fUseBayesPID(kFALSE),
140 fUseTrackIndexCut(kTRUE),
141 fUsePsiPairCut(kTRUE),
142 fUseMassCut(kFALSE),
143 fUseGammaCut(kFALSE),
144 fReadMagFieldSign(kTRUE),
145 fUseAliKF(kFALSE),
146 fMagFieldSign(1),
b9295e41 147 fkElectronMass(0.00051099891),
7f3c7cc6 148 fPsiPairCut(0.45),
149 fDeltaPhiCutMin(0.),
150 fDeltaPhiCutMax(0.12),
151 fMassCutMin(0.),
152 fMassCutMax(0.1),
153 fNSigmaBelowElecTPCbethe(-2.),
154 fNSigmaAboveElecTPCbethe(3.),
155 fNSigmaAbovePionTPCbethe(3.),
156 fNSigmaAboveKaonTPCbethe(3.),
157 fNSigmaAboveProtonTPCbethe(3.),
158 fTrkSelectionCriteria(kGlobalTrack)
159{
160 // Common I/O in slot 0
161 DefineInput (0, TChain::Class());
162
163 // Your private output
164 DefineOutput(1, TList::Class());
165// DefineOutput(2, AliCFContainer::Class()); // for CF
166
167 AdoptITSsaTrackCuts();
168 AdoptESDtrackCuts();
169 AdoptESDpidCuts();
b9295e41 170 // fkElectronMass = TDatabasePDG::Instance()->GetParticle( ::kElectron )->Mass(); //
7f3c7cc6 171
172 fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
173 fGammaPool->SetOwner(kTRUE);
174}
175
176//-----------------------------------------------------------------------------------------------
177AliAnalysisTaskGammaConvDalitz::~AliAnalysisTaskGammaConvDalitz()
178{
7811d09e 179//
180// virtual destructor
181//
182
7f3c7cc6 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;
191}
192
193//-----------------------------------------------------------------------------------------------
194void AliAnalysisTaskGammaConvDalitz::ConnectInputData(Option_t *option)
195{
7811d09e 196//
197// Connect Input Data
198//
7f3c7cc6 199 if( fDebug ) AliInfo("=> ConnectInputData");
200
201 AliAnalysisTaskSE::ConnectInputData(option);
202
203 if( fV0Reader == 0 )
204 {
205 AliFatal("There is not pointer to AliV0Reader object!!!");
206 }
207}
208
209//-----------------------------------------------------------------------------------------------
210void AliAnalysisTaskGammaConvDalitz::UserCreateOutputObjects()
211{
7811d09e 212//
213// Create ouput objects
214//
7f3c7cc6 215 if( fDebug ) AliInfo("=> UserCreateOutputObjects");
216
217 // Create the output container
218 if( fOutputContainer != 0 )
219 {
220 delete fOutputContainer;
221 }
222
223 fOutputContainer = new TList();
224
225 // Add the histograms to the output container
226 fHistograms->GetOutputContainer( fOutputContainer );
227 fOutputContainer->SetOwner(kTRUE);
228
229 PostData( 1, fOutputContainer );
230}
231
232//-----------------------------------------------------------------------------------------------
7f3c7cc6 233void AliAnalysisTaskGammaConvDalitz::UserExec(Option_t */*option*/)
234{
7811d09e 235//
236// Execute analysis for current event
237//
2eedd4ed 238
239
240
7f3c7cc6 241 if( fDebug ) AliInfo("=> UserExec");
242
243 if( fV0Reader == 0 )
244 {
245 AliFatal("no pointer to AliV0Reader");
246 return;
247 }
248
7f3c7cc6 249 // Create list of gamma candidates in standalone mode
250 // otherwise use the created ones by AliAnalysisTaskGammaConversion
251 if( fStandalone )
252 {
7811d09e 253
254 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
255 AliESDInputHandler *esdHandler=0;
256 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
257 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
258 } else {
259 //load esd pid bethe bloch parameters depending on the existance of the MC handler
260 // yes: MC parameters
261 // no: data parameters
262 if (!AliV0Reader::GetESDpid()){
9cb96e64 263 if (MCEvent() ) {
7811d09e 264 AliV0Reader::InitESDpid();
265 } else {
266 AliV0Reader::InitESDpid(1);
267 }
268 }
269 }
270
271
d1758e4f 272 if (MCEvent() ) {
7811d09e 273
274 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
275 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
276
277 if (!mcHandler){
278 AliError("Could not retrive MC event handler!");
279 return;
280 }
281
282 if (!mcHandler->InitOk() ){
283
284 return;
285 }
286 if (!mcHandler->TreeK() ){
287 return;
288 }
289 if (!mcHandler->TreeTR() ) {
290 return;
291 }
292 }
293
294
295
296
297 fV0Reader->SetInputAndMCEvent( InputEvent(), MCEvent() );
7f3c7cc6 298 fV0Reader->Initialize();
299 }
300
301 if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
302 {
303 if( fDebug ) AliInfo("no contributors to primary vertex");
304 return;
305 }
306
7811d09e 307 if( fV0Reader->CheckForPrimaryVertexZ() == kFALSE )
308 {
309
310 if( fDebug ) AliInfo("z vertex out of range");
311 return;
312 }
313
7f3c7cc6 314 // Get Pointers
315 fBGEventHandler = fV0Reader->GetBGHandler();
316 fESDpid = fV0Reader->GetESDpid();
317 fESDEvent = fV0Reader->GetESDEvent();
18922343 318 if(fDoMC && MCEvent())
7f3c7cc6 319 {
7811d09e 320 fStack= MCEvent()->Stack();
321 fGCMCEvent=MCEvent();
7f3c7cc6 322 }
323
324 // Read the magnetic field sign from ESD
325 if ( fReadMagFieldSign == kTRUE )
326 {
327 fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
328 }
329
330 // Process MC information
331 if(fDoMC)
332 {
333 ProcessMCData();
334 }
335
336 if( fStandalone ){
337 while(fV0Reader->NextV0()){}; //SelectGammas
338 fV0Reader->ResetV0IndexNumber();
339 }
7811d09e 340
7f3c7cc6 341 CreateListOfDalitzPairCandidates();
342 ProcessGammaElectronsForDalitzAnalysis();
7811d09e 343
344 if ( fStandalone ){
345
346 fV0Reader->UpdateEventByEventData();
347
348 }
7f3c7cc6 349
350 PostData( 1, fOutputContainer );
351}
352
353
354void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
355{
7811d09e 356//
7f3c7cc6 357 if( fDebug ) AliInfo("Not to do anything in Terminate");
358}
359
360//-----------------------------------------------------------------------------------------------
7f3c7cc6 361void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
362{
7811d09e 363//
364// set user ITSsa track cuts
365//
7f3c7cc6 366 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
367
368 if( esdCuts )
369 {
370 fITSsaTrackCuts = esdCuts;
371 }
372 else
373 {
7811d09e 374 // default cuts
7f3c7cc6 375 fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
376
377 fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
378 fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
379
380 fITSsaTrackCuts->SetMinNClustersITS(2);
381 fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
382 fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
383
384 fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
385 fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
386
387 fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
388 }
389}
390
391//-----------------------------------------------------------------------------------------------
7f3c7cc6 392void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
393{
7811d09e 394//
395// set user global track cuts
396//
7f3c7cc6 397 if( fESDtrackCuts ) delete fESDtrackCuts;
398
399 if( esdCuts )
400 {
401 fESDtrackCuts = esdCuts;
402 }
403 else
404 {
7811d09e 405 //default cuts
7f3c7cc6 406 fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
407
408 fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
409 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
410
411 fESDtrackCuts->SetMinNClustersITS(2);
412 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
413 fESDtrackCuts->SetRequireITSRefit(kTRUE);
414
415 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
416 fESDtrackCuts->SetMaxNsigmaToVertex(3);
417
418 fESDtrackCuts->SetMinNClustersTPC(80);
419 fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
420 fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
421 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
422 }
423}
424
425//-----------------------------------------------------------------------------------------------
426void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
427{
7811d09e 428//
429// set user pid cuts
430//
7f3c7cc6 431 if( fESDpidCuts ) delete fESDpidCuts;
432 if( esdPIDCuts )
433 {
434 fESDpidCuts = esdPIDCuts;
435 }
7811d09e 436 else // default cuts
7f3c7cc6 437 {
438 fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
439 // fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
440 fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
441 }
442}
443
444//-----------------------------------------------------------------------------------------------
7f3c7cc6 445void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
446{
7811d09e 447//
448// Process generation
449//
7f3c7cc6 450 if( fDebug ) AliInfo("=> ProcessMCData");
451
452 fHistograms->FillTable("Table_Generation", 0); //number of events
453
454 for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
455 {
456 TParticle* iParticle = fStack->Particle( i );
457 if( !iParticle ) continue;
458
459 if ( i >= fStack->GetNprimary() ) continue; // is primary?
460 if ( iParticle->GetPdgCode() != ::kPi0 ) continue; // is Pi0?
461
462 if( iParticle->GetNDaughters() == 2 &&
463 fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
464 fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
465 {
466 fHistograms->FillTable("Table_Generation", 1); // pi0 -> gg
467 }
468
469 if ( iParticle->GetNDaughters() != 3 ) continue; // Num == 3 (e+,e-,gamma)
470
471 // Check for Pi0 Dalitz decay
472 TParticle* eposPi0 = 0;
473 TParticle* enegPi0 = 0;
474 TParticle* gammaPi0 = 0;
475
476 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
477 {
478 switch(fStack->Particle(idxPi0)->GetPdgCode())
479 {
480 case ::kPositron:
481 eposPi0 = fStack->Particle(idxPi0);
482 break;
483 case ::kElectron:
484 enegPi0 = fStack->Particle(idxPi0);
485 break;
486 case ::kGamma:
487 gammaPi0 = fStack->Particle(idxPi0);
488 break;
489 }
490 }
491
492 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
493
494 // found a Pi0 Dalitz decay
495
496 fHistograms->FillTable("Table_Generation", 2);
497 fHistograms->FillHistogram("MC_Pi0Dalitz_P", iParticle->P());
498 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt", iParticle->Pt());
499 fHistograms->FillHistogram("MC_Pi0Dalitz_Eta", iParticle->Eta());
500 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
501 fHistograms->FillHistogram("MC_EposDalitz_Pt", eposPi0->Pt());
502 fHistograms->FillHistogram("MC_EposDalitz_Eta", eposPi0->Eta());
503 fHistograms->FillHistogram("MC_EnegDalitz_Pt", enegPi0->Pt());
504 fHistograms->FillHistogram("MC_EnegDalitz_Eta", enegPi0->Eta());
505 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Pt", gammaPi0->Pt());
506 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Eta", gammaPi0->Eta());
507
508 // Angle between the gamma and the plane e+e-
509 TVector3 ePosMom( eposPi0->Px(), eposPi0->Py(), eposPi0->Pz() );
510 TVector3 eNegMom( enegPi0->Px(), enegPi0->Py(), enegPi0->Pz() );
511 TVector3 gamMom( gammaPi0->Px(), gammaPi0->Py() , gammaPi0->Pz() );
512 TVector3 planeEposEneg = eNegMom.Cross( ePosMom );
513 Double_t anglePlaneGamma = planeEposEneg.Angle(gamMom);
514
515 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
516
517 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle", anglePlaneGamma);
518 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_P", anglePlaneGamma, gammaPi0->P());
519 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_Pt", anglePlaneGamma, gammaPi0->Pt());
520
521
522 // check for gamma conversion
523 Bool_t daugGammaElectron = kFALSE;
524 Bool_t daugGammaPositron = kFALSE; // acceptance
525 Bool_t daugGammaElectronAll = kFALSE;
526 Bool_t daugGammaPositronAll = kFALSE;
527
528 // is the gamma converted? -> has 2 daughter e+e-
529 // are e+ e- from gamma in the acceptance for the V0s
530 if( gammaPi0->GetNDaughters() >= 2 )
531 {
532 for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
533 {
534 TParticle* tmpDaughter = fStack->Particle(tIndex);
7811d09e 535
7f3c7cc6 536 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
537 if( tmpDaughter->GetPdgCode() == ::kElectron )
538 { // e+
539 daugGammaElectronAll = kTRUE;
7811d09e 540
7f3c7cc6 541 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
542 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
543 (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
544 {
545 daugGammaElectron = kTRUE;
546 }
547 }
548 else if( tmpDaughter->GetPdgCode() == ::kPositron )
549 {
550 daugGammaPositronAll = kTRUE;
551 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
552 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
553 (tmpDaughter->R() < fV0Reader->GetMaxRCut() ) )
554 {
555 daugGammaPositron = kTRUE;
556 }
557 }
558 }
559 }
560
561
562 if( daugGammaElectronAll && daugGammaPositronAll )
563 {
564 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
565 fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
566 }
567
568 Float_t etaMin, etaMax;
569 fESDtrackCuts->GetEtaRange( etaMin, etaMax );
570
571 // e+e- pair within acceptance
572 if ( TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
573 {
574 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt", eposPi0->Pt());
575 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta", eposPi0->Eta());
576 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt", enegPi0->Pt());
577 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta", enegPi0->Eta());
578 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_EposPt_vs_EnegPt", eposPi0->Pt(), enegPi0->Pt());
579 }
580
581 // Pi0 (e+e-gamma) within acceptance
582 //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
583
584 if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
585 TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
586 {
587 fHistograms->FillTable("Table_Generation",3); //
588 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt",iParticle->Pt());
589 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Eta",iParticle->Eta());
590 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Radius",iParticle->R());
591 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Pt",gammaPi0->Pt());
592 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Eta",gammaPi0->Eta());
7811d09e 593 fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Pt",eposPi0->Pt());
594 fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Eta",eposPi0->Eta());
595 fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Pt",enegPi0->Pt());
596 fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Eta",enegPi0->Eta());
7f3c7cc6 597 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
598 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
599
600 // Pi0 within acceptance with gamma converted
601
602 if ( daugGammaElectron && daugGammaPositron )
603 {
604 fHistograms->FillTable("Table_Generation",4); //
605
606 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
607 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
7811d09e 608 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Pt",eposPi0->Pt());
609 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Eta",eposPi0->Eta());
610 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Pt",enegPi0->Pt());
611 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Eta",enegPi0->Eta());
7f3c7cc6 612 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Pt",gammaPi0->Pt());
613 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Eta",gammaPi0->Eta());
614 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle",anglePlaneGamma);
615 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle_vs_Pt",anglePlaneGamma,gammaPi0->Pt());
616 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
617 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
618 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
619 }
620 }
621 }
622}
623
624//-----------------------------------------------------------------------------------------------
7f3c7cc6 625void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
626{
7811d09e 627//
628// Dalitz pair candidates
629//
7f3c7cc6 630 if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
631
632 fEposCandidateIndex.clear();
633 fEnegCandidateIndex.clear();
634
635 fHistograms->FillTable("Table_Cuts", 0);
636
637 for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
638 {
639 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
640 if ( !iTrack ) continue;
641
642
643 Double_t p[3];
644
645 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
646
647 TVector3 iMom(p[0],p[1],p[2]);
648
649 //
650 // Check track cuts and find track type
651 //
652
653 Bool_t isTrackAccepted = 0;
654 Int_t trackType = -1;
655 switch(fTrkSelectionCriteria)
656 {
657 case kITSsaTrack:
658 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
659 trackType = kITSsaTrack;
660 break;
661
662 case kGlobalTrack:
663 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
664 trackType = kGlobalTrack;
665 break;
666
667 case kITSsaGlobalTrack:
668 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
669 {
670 isTrackAccepted = kTRUE;
671 if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
672 else trackType = kGlobalTrack;
673 }
674 break;
675 }
676
677 if(!isTrackAccepted) continue;
678
679 //
680 // PID
681 //
682
683 Int_t pid=-1;
684 Int_t pidMC=-1;
685
686 if(fUseBayesPID)
687 {
688 pid = GetBayesPid(iTrack,trackType);
689 }
690 else
691 {
692 pid = GetNSigmaPid(iTrack,trackType);
693 }
694
695 if( fDoMC )
696 {
697 pidMC = GetMonteCarloPid(iTrack);
698 // pid table
699 Int_t iLabel = TMath::Abs(iTrack->GetLabel());
700 TParticle* iParticle = fStack->Particle(iLabel);
701 FillPidTable(iParticle, pid);
702 }
7811d09e 703
7f3c7cc6 704 // ITS standalone tracks
705 if( trackType == kITSsaTrack)
706 {
707 Double_t mom = iTrack->GetP();
708 Double_t signal = iTrack->GetITSsignal();
709
710 fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
711
712 if( pid == AliPID::kElectron )
713 {
714
715 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
716 if(fDoMC && pid == pidMC)
717 {
718 fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
719 }
720 }
721
722 if( fDoMC && pidMC == AliPID::kElectron)
723 {
724 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
725 }
726 }
727
728 else // global tracks
729 {
730 const AliExternalTrackParam *in = iTrack->GetInnerParam();
731 Double_t mom = in->GetP();
732 Double_t signal = iTrack->GetTPCsignal();
733
734 fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
735
736 if( fDoMC && pidMC == AliPID::kElectron )
737 {
738 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
739 }
740
741 if( pid == AliPID::kElectron )
742 {
743 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
744 if(fDoMC && pid == pidMC)
745 {
746 fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
747 }
748 }
749 }
750
751 if( AliPID::kElectron != pid) continue;
752
753 // electron track candidates from here
754
755 if( iTrack->GetSign() > 0 )
756 {
757 fEposCandidateIndex.push_back(i);
758 }
759 else
760 {
761 fEnegCandidateIndex.push_back(i);
762 }
763 }
764
765 // gamma candidates
766 GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
767
768 if(fDoMC)
769 {
2eedd4ed 770 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 771 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
772 delete pi0Dalitz;
773 }
774
775 if(fUseTrackIndexCut) // remove repeated tracks
776 {
777 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
778
779 if(fDoMC)
780 {
2eedd4ed 781 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 782 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
783 delete pi0Dalitz;
784 }
785 }
786
787 if(fUsePsiPairCut) // remove electrons from gamma conversions
788 {
789 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
790
791 if(fDoMC)
792 {
2eedd4ed 793 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 794 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
795 delete pi0Dalitz;
796 }
797 }
798
799 if( fUseMassCut )
800 {
801 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
802
803 if(fDoMC)
804 {
2eedd4ed 805 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 806 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
807 delete pi0Dalitz;
808 }
809 }
810
811 if(fUseGammaCut)
812 {
813 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
814
815 if(fDoMC)
816 {
2eedd4ed 817 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 818 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
819 delete pi0Dalitz;
820 }
821 }
822}
823
824//-----------------------------------------------------------------------------------------------
825void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
826{
7811d09e 827//
828// Process gamma and electrons for pi0 Dalitz decay
829//
7f3c7cc6 830 if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
2eedd4ed 831
7f3c7cc6 832
833 fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
2eedd4ed 834
835
7f3c7cc6 836
837 TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
838
839 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
840 {
841 AliKFParticle* epos = (AliKFParticle*) ePosCandidates->At(i);
842 fHistograms->FillHistogram("ESD_EposCandidates_Pt", epos->GetPt());
843 fHistograms->FillHistogram("ESD_EposCandidates_Eta", epos->GetEta());
844 fHistograms->FillTable( "Table_Reconstruction", 1);
845 }
846
847 TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
848
849 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
850 {
851 AliKFParticle* eneg = (AliKFParticle*) eNegCandidates->At(i);
852 fHistograms->FillHistogram("ESD_EnegCandidates_Pt", eneg->GetPt());
853 fHistograms->FillHistogram("ESD_EnegCandidates_Eta", eneg->GetEta());
854 fHistograms->FillTable( "Table_Reconstruction", 2);
855 }
856
857 TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
858 for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
859 {
860 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
861
862 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
2eedd4ed 863 fHistograms->FillHistogram("ESD_DalitzPairCandidates_InvMass", dalitz->M());
864 fHistograms->FillHistogram("ESD_DalitzPairCandidates_InvMass_vs_Pt",dalitz->M(),dalitz->Pt());
7f3c7cc6 865 }
866
867 // gamma candidates
868 for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
869 {
870 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
871 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
872 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
873 }
874
875 // psi pair for all candidates
876 //if(fUsePsiPairCut)
877 FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
878
879 // Angle epos,eneg gamma
880 FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
881 FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
882
2eedd4ed 883 TClonesArray* pi0Candidates = FindParticleDalitz(ePosCandidates, eNegCandidates, fGammaCandidates,0);
7f3c7cc6 884 for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
885 {
886 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
887
888 fHistograms->FillHistogram("ESD_Pi0_P", pi0->P());
889 fHistograms->FillHistogram("ESD_Pi0_Pt", pi0->Pt());
890 fHistograms->FillHistogram("ESD_Pi0_Eta", pi0->Eta());
891 fHistograms->FillHistogram("ESD_Pi0_Y", pi0->Rapidity());
892 fHistograms->FillHistogram("ESD_Pi0_Phi", pi0->Phi());
2eedd4ed 893 fHistograms->FillHistogram("ESD_Pi0_Pt_vs_Y",pi0->Pt(),pi0->Rapidity());
894 fHistograms->FillHistogram("ESD_Pi0_InvMass", pi0->M());
895 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Pt", pi0->M(),pi0->Pt());
896 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Y",pi0->M(),pi0->Rapidity());
897 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Eta",pi0->M(),pi0->Eta());
7f3c7cc6 898 }
2eedd4ed 899
900 for(Int_t iPos=0; iPos < ePosCandidates->GetEntriesFast(); ++iPos)
901 {
902 AliKFParticle* lPosKF = (AliKFParticle*)ePosCandidates->At(iPos);
903
904 for(Int_t iNeg=0; iNeg < eNegCandidates->GetEntriesFast(); ++iNeg)
905 {
906 AliKFParticle* lNegKF = (AliKFParticle*)eNegCandidates->At(iNeg);
907 AliKFParticle lPosNeg(*lPosKF,*lNegKF );
908
909 for(Int_t iGam=0; iGam < fGammaCandidates->GetEntriesFast(); ++iGam)
910 {
911 AliKFParticle* lGamKF = (AliKFParticle*)fGammaCandidates->At(iGam);
912
913 AliKFParticle lPosNegGam( *lPosKF, *lNegKF, *lGamKF );
914
915 Double_t lDiffMass = lPosNegGam.GetMass() - lPosNeg.GetMass();
916
917 fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_Diff",lDiffMass );
918 fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_vs_Pt_Diff",lDiffMass,lPosNegGam.GetPt());
919 fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass",lPosNegGam.GetMass());
920 fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_vs_Pt",lPosNegGam.GetMass(),lPosNegGam.GetPt());
921
922 }
923 }
924 }
925
926
927
7f3c7cc6 928
929 delete dalitzPairCandidates;
930 delete pi0Candidates;
931
932 if(fComputeBkg)
933 {
2eedd4ed 934
935 // 1) e+e- dalitz
936 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
937 {
938 AliKFParticle* epos1 = (AliKFParticle*) ePosCandidates->At(i);
939
940 for(Int_t j=i+1; j < ePosCandidates->GetEntriesFast(); ++j)
941 {
942 AliKFParticle* epos2 = (AliKFParticle*) ePosCandidates->At(j);
943 AliKFParticle ePosePos( *epos1,*epos2 );
944 fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass",ePosePos.GetMass());
945 fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass_vs_Pt",ePosePos.GetMass(),ePosePos.GetPt());
946
947
948 }
949
950
951 }
952 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
953 {
954 AliKFParticle* eneg1 = (AliKFParticle*) eNegCandidates->At(i);
955
956 for(Int_t j=i+1; j < eNegCandidates->GetEntriesFast(); ++j)
957 {
958 AliKFParticle* eneg2 = (AliKFParticle*) eNegCandidates->At(j);
959 AliKFParticle eNegeNeg( *eneg1,*eneg2 );
960 fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass",eNegeNeg.GetMass());
961 fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass_vs_Pt",eNegeNeg.GetMass(),eNegeNeg.GetPt());
962 }
963
964
965 }
966
7f3c7cc6 967 // 1) e+e- with with gammas used in the signal
2eedd4ed 968 TClonesArray* pi0Bkg = FindParticleDalitz(ePosCandidates, eNegCandidates, fGammaPool,1);
7f3c7cc6 969
970 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
971 {
972 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
2eedd4ed 973 fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass", pi0->M());
974 fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_vs_Pt",pi0->M(),pi0->Pt());
7f3c7cc6 975 }
2eedd4ed 976 ///////////////////////////////Temporal for Dalitz
977
978
979
980
981
982
7f3c7cc6 983
984 if(ePosCandidates->GetEntriesFast() > 0 &&
985 eNegCandidates->GetEntriesFast() > 0 &&
986 fGammaCandidates->GetEntriesFast() > 0)
987 {
988 UpdateGammaPool(fGammaCandidates);
989 }
990
991 delete pi0Bkg;
992
993 // 2) e+e- with gammas from a pool of events
7811d09e 994 TClonesArray* gammaBGHandler = GammasFromBGHandler();
2eedd4ed 995 pi0Bkg = FindParticleDalitz(ePosCandidates, eNegCandidates, gammaBGHandler,2);
7f3c7cc6 996
997 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
998 {
999 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
2eedd4ed 1000 fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass", pi0->M());
1001 fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_vs_Pt",pi0->M(), pi0->Pt());
7f3c7cc6 1002 }
1003
1004 delete pi0Bkg;
1005
1006 // 3) e+ with e-, gamma from a pool of events
1007 TClonesArray* elecBGHandler = ElectronFromBGHandler();
2eedd4ed 1008 pi0Bkg = FindParticleDalitz(ePosCandidates, elecBGHandler, gammaBGHandler,3);
7f3c7cc6 1009
1010 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
1011 {
1012 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
2eedd4ed 1013 fHistograms->FillHistogram("ESD_BKG_Electron_InvMass", pi0->M());
1014 fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_vs_Pt",pi0->M(), pi0->Pt());
7f3c7cc6 1015 }
1016
1017 if(eNegCandidates->GetEntriesFast() > 0)
1018 {
1019 UpdateElectronPool(eNegCandidates);
1020 }
1021
1022 delete gammaBGHandler;
1023 delete elecBGHandler;
1024 delete pi0Bkg;
7811d09e 1025
7f3c7cc6 1026 }
1027
1028 delete ePosCandidates;
1029 delete eNegCandidates;
1030
1031 if(fDoMC)
1032 {
1033 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
1034 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
1035 {
1036 AliKFParticle* epos = (AliKFParticle*) ePosPi0Dalitz->At(i);
2eedd4ed 1037 fHistograms->FillHistogram("MC_ESD_Pi0_EposDalitz_Pt", epos->GetPt());
1038 fHistograms->FillHistogram("MC_ESD_Pi0_EposDalitz_Eta", epos->GetEta());
7f3c7cc6 1039 fHistograms->FillTable( "Table_Reconstruction", 3);
1040 }
1041
1042 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
1043 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
1044 {
1045 AliKFParticle* eneg = (AliKFParticle*) eNegPi0Dalitz->At(i);
2eedd4ed 1046 fHistograms->FillHistogram("MC_ESD_Pi0_EnegDalitz_Pt", eneg->GetPt());
1047 fHistograms->FillHistogram("MC_ESD_Pi0_EnegDalitz_Eta", eneg->GetEta());
7f3c7cc6 1048 fHistograms->FillTable( "Table_Reconstruction", 4);
1049 }
1050
2eedd4ed 1051 TClonesArray* dalitzPairPi0 = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex,1);
1052 for(Int_t i=0; i < dalitzPairPi0->GetEntriesFast(); ++i)
7f3c7cc6 1053 {
2eedd4ed 1054 TLorentzVector* dalitz = (TLorentzVector*) dalitzPairPi0->At(i);
1055 fHistograms->FillHistogram("MC_ESD_Pi0_DalitzPair_Pt", dalitz->Pt());
1056 fHistograms->FillHistogram("MC_ESD_Pi0_DalitzPair_Mass", dalitz->M());
7f3c7cc6 1057 fHistograms->FillHistogram( "Table_Reconstruction", 5 );
1058 }
2eedd4ed 1059
1060 TClonesArray* dalitzPairEta = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex,2);
1061 for(Int_t i=0; i < dalitzPairEta->GetEntriesFast(); ++i)
1062 {
1063 TLorentzVector* dalitz = (TLorentzVector*) dalitzPairEta->At(i);
1064 fHistograms->FillHistogram("MC_ESD_Eta0_DalitzPair_Pt", dalitz->Pt());
1065 fHistograms->FillHistogram("MC_ESD_Eta0_DalitzPair_InvMass", dalitz->M());
1066
1067 }
1068
1069 TClonesArray* lJpsiAll = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,-1);
1070
1071 for(Int_t i=0; i < lJpsiAll->GetEntriesFast(); ++i)
1072 {
1073 TLorentzVector* jpsi = (TLorentzVector*) lJpsiAll->At(i);
1074 fHistograms->FillHistogram("MC_ESD_Jpsi_Pt",jpsi->Pt());
1075 fHistograms->FillHistogram("MC_ESD_Jpsi_InvMass",jpsi->M());
1076 fHistograms->FillHistogram("MC_ESD_Jpsi_InvMass_vs_Pt",jpsi->M(),jpsi->Pt());
1077
1078 }
1079
1080
1081 TClonesArray* lJpsiChic0 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,0);
1082
1083 for(Int_t i=0; i < lJpsiChic0->GetEntriesFast(); ++i)
1084 {
1085 TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic0->At(i);
1086 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic0_Pt",jpsi->Pt());
1087 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic0_InvMass",jpsi->M());
1088
1089 }
1090 TClonesArray* lJpsiChic1 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,1);
1091
1092 for(Int_t i=0; i < lJpsiChic1->GetEntriesFast(); ++i)
1093 {
1094 TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic1->At(i);
1095 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic1_Pt",jpsi->Pt());
1096 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic1_InvMass",jpsi->M());
1097
1098 }
1099 TClonesArray* lJpsiChic2 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,2);
1100
1101 for(Int_t i=0; i < lJpsiChic2->GetEntriesFast(); ++i)
1102 {
1103 TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic2->At(i);
1104 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic2_Pt",jpsi->Pt());
1105 fHistograms->FillHistogram("MC_ESD_Jpsi_Chic2_InvMass",jpsi->M());
1106
1107 }
7f3c7cc6 1108
1109 // psi pair for dalitz pairs
1110 //if(fUsePsiPairCut)
2eedd4ed 1111 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_Pi0_DalitzPair_PsiPair_vs_DPhi");
7f3c7cc6 1112
1113 delete ePosPi0Dalitz;
1114 delete eNegPi0Dalitz;
2eedd4ed 1115 delete dalitzPairPi0;
1116 delete dalitzPairEta;
1117 delete lJpsiAll;
1118 delete lJpsiChic0;
1119 delete lJpsiChic1;
1120 delete lJpsiChic2;
7f3c7cc6 1121 // all gammas
1122 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1123 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1124 {
1125 AliKFParticle* iGamma = (AliKFParticle*) gamma->At(i);
1126 fHistograms->FillHistogram("MC_ESD_Gamma_Pt", iGamma->GetPt());
1127 fHistograms->FillHistogram("MC_ESD_Gamma_Eta", iGamma->GetEta());
1128 }
1129
1130 delete gamma;
1131
1132 // gamma from pi0 dalitz
1133 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1134 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
1135 {
1136 AliKFParticle* iGamma = (AliKFParticle*) gammaPi0Dalitz->At(i);
1137 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Pt", iGamma->GetPt());
1138 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Eta", iGamma->GetEta());
1139 fHistograms->FillTable( "Table_Reconstruction", 6);
1140 }
1141
1142 delete gammaPi0Dalitz;
1143
2eedd4ed 1144 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
7f3c7cc6 1145
1146 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
1147 {
1148 TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
1149
1150 fHistograms->FillHistogram("MC_ESD_Pi0_P", pi0->P());
1151 fHistograms->FillHistogram("MC_ESD_Pi0_Pt", pi0->Pt());
1152 fHistograms->FillHistogram("MC_ESD_Pi0_Eta", pi0->Eta());
1153 fHistograms->FillHistogram("MC_ESD_Pi0_Y", pi0->Rapidity());
1154 fHistograms->FillHistogram("MC_ESD_Pi0_Phi", pi0->Phi());
2eedd4ed 1155 fHistograms->FillHistogram("MC_ESD_Pi0_Y_vs_Pt",pi0->Pt(), pi0->Rapidity());
1156 fHistograms->FillHistogram("MC_ESD_Pi0_InvMass", pi0->M());
1157 fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Pt",pi0->M(),pi0->Pt());
1158 fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Y", pi0->M(), pi0->Rapidity());
1159 fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Eta", pi0->M(),pi0->Eta());
7f3c7cc6 1160 fHistograms->FillHistogram( "Table_Reconstruction", 7);
1161 }
1162 delete pi0Dalitz;
2eedd4ed 1163
1164
1165 TClonesArray* eta0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,2);
1166
1167
1168 for(Int_t i=0; i < eta0Dalitz->GetEntriesFast(); ++i)
1169 {
1170 TLorentzVector* eta0 = (TLorentzVector*) eta0Dalitz->At(i);
1171
1172 fHistograms->FillHistogram("MC_ESD_Eta0_P", eta0->P());
1173 fHistograms->FillHistogram("MC_ESD_Eta0_Pt", eta0->Pt());
1174 fHistograms->FillHistogram("MC_ESD_Eta0_Eta", eta0->Eta());
1175 fHistograms->FillHistogram("MC_ESD_Eta0_Y", eta0->Rapidity());
1176 fHistograms->FillHistogram("MC_ESD_Eta0_Phi", eta0->Phi());
1177 fHistograms->FillHistogram("MC_ESD_Eta0_Pt_vs_Y", eta0->Pt(),eta0->Rapidity());
1178 fHistograms->FillHistogram("MC_ESD_Eta0_InvMass", eta0->M());
1179 fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Pt", eta0->M(), eta0->Pt());
1180 fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Y", eta0->M(), eta0->Rapidity());
1181 fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Eta",eta0->M(),eta0->Eta());
1182 }
1183 delete eta0Dalitz;
1184
1185
1186 TClonesArray* chic0Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,0);
1187
1188 for(Int_t i=0; i < chic0Array->GetEntriesFast(); ++i)
1189 {
1190 TLorentzVector* chic0 = (TLorentzVector*) chic0Array->At(i);
1191 fHistograms->FillHistogram("MC_ESD_Chic0_InvMass", chic0->M());
1192 fHistograms->FillHistogram("MC_ESD_Chic0_InvMass_vs_Pt", chic0->M(),chic0->Pt());
1193 }
1194 delete chic0Array;
1195
1196 TClonesArray* chic1Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
1197
1198 for(Int_t i=0; i < chic1Array->GetEntriesFast(); ++i)
1199 {
1200 TLorentzVector* chic1 = (TLorentzVector*) chic1Array->At(i);
1201 fHistograms->FillHistogram("MC_ESD_Chic1_InvMass", chic1->M());
1202 fHistograms->FillHistogram("MC_ESD_Chic1_InvMass_vs_Pt", chic1->M(), chic1->Pt());
1203 }
1204 delete chic1Array;
1205
1206 TClonesArray* chic2Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,2);
1207
1208 for(Int_t i=0; i < chic2Array->GetEntriesFast(); ++i)
1209 {
1210 TLorentzVector* chic2 = (TLorentzVector*) chic2Array->At(i);
1211 fHistograms->FillHistogram("MC_ESD_Chic2_InvMass", chic2->M());
1212 fHistograms->FillHistogram("MC_ESD_Chic2_InvMass_vs_Pt", chic2->M(), chic2->Pt());
1213 }
1214 delete chic2Array;
1215
7f3c7cc6 1216
1217 // psi pair for electrons from gamma conversions assuming they came from main vertex
1218 // if(fUsePsiPairCut)
1219 for(UInt_t i=0; i < fEposCandidateIndex.size(); ++i)
1220 {
1221 AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
1222 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1223
1224 for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
1225 {
1226 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
1227 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1228
1229 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
1230
1231 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1232 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1233
1234 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
1235 }
1236 }
1237 // FIXME: eta -> e+e-gamma
1238 }
1239}
1240
1241//--------------------------------------------------------------------------
1242Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
1243{
7811d09e 1244//
1245// Get rapidity
1246//
1247 const double kEPSILON=1.e-16;
7f3c7cc6 1248
7811d09e 1249 if(p->Energy() - TMath::Abs(p->Pz()) < kEPSILON )
7f3c7cc6 1250 {
1251 return 1.e10;
1252 }
1253 return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
1254}
1255
1256//--------------------------------------------------------------------------
1257void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1258{
7811d09e 1259//
1260// Fill histogram with psipair(pos,neg)
1261//
7f3c7cc6 1262 for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1263 {
1264 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1265 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1266 {
1267 AliKFParticle* negKF = (AliKFParticle*) neg->At(j);
1268 Double_t psiPair = GetPsiPair(posKF, negKF);
1269 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negKF->GetPhi() - posKF->GetPhi());
1270 fHistograms->FillHistogram(hName, deltaPhi, psiPair);
1271 }
1272 }
1273}
1274
1275//--------------------------------------------------------------------------
1276void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1277{
7811d09e 1278//
1279// Fill histogram with angle(x,y)
1280//
7f3c7cc6 1281 for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1282 {
1283 AliKFParticle* xKF = (AliKFParticle*) x->At(i);
1284 TVector3 xMom(xKF->Px(),xKF->Py(),xKF->Pz());
1285 for( Int_t j=0; j < y->GetEntriesFast(); ++j )
1286 {
1287 AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1288 TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1289 fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1290 }
1291 }
1292}
1293
1294//--------------------------------------------------------------------------
7f3c7cc6 1295void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1296{
7811d09e 1297//
1298// Fill table with pid info
1299//
7f3c7cc6 1300 Int_t iGen=-1;
1301 switch(TMath::Abs(p->GetPdgCode()))
1302 {
1303 case ::kElectron: iGen=0; break;
1304 case ::kMuonMinus: iGen=1; break;
1305 case ::kPiPlus: iGen=2; break;
1306 case ::kKPlus: iGen=3; break;
1307 case ::kProton: iGen=4; break;
1308 default: iGen=-1;
1309 }
1310
1311 int jRec=-1;
1312 if(pid > -1 && pid < 5) jRec = pid;
1313
1314 if ((iGen > -1) && (jRec > -1))
1315 {
1316 fHistograms->FillTable("Table_PID", iGen, jRec);
1317 // sum
1318 fHistograms->FillTable("Table_PID", iGen, 5);
1319 fHistograms->FillTable("Table_PID", 5, jRec);
1320 }
1321}
1322
1323//--------------------------------------------------------------------------
1324void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1325{
7811d09e 1326//
1327// Make a copy of gamma candidates from V0reader
1328//
7f3c7cc6 1329 posIndex.clear();
1330 negIndex.clear();
1331
1332 if(gamma) delete gamma;
1333
1334 TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1335
1336 gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1337 gamma->SetOwner(kTRUE);
1338
1339 // make a copy
1340 for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1341 {
1342 AliKFParticle* gamKF = (AliKFParticle*)gammaV0->At(i);
1343 new ((*gamma)[i]) AliKFParticle(*gamKF);
1344 posIndex.push_back(fV0Reader->GetPindex(i));
1345 negIndex.push_back(fV0Reader->GetNindex(i));
1346 }
1347}
1348
1349//--------------------------------------------------------------------------
1350TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1351{
7811d09e 1352//
1353// Convert track index vector to AliKFParticle array
1354//
7f3c7cc6 1355 TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1356 indexKF->SetOwner(kTRUE);
1357
1358 for(UInt_t i = 0; i < index.size(); ++i)
1359 {
1360 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1361 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1362 }
1363
1364 return indexKF;
1365}
1366
1367//--------------------------------------------------------------------------
1368TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1369{
7811d09e 1370//
1371// Find true electrons from pi0 Dalitz decay candidates with MC
1372//
7f3c7cc6 1373 TClonesArray* elec = new TClonesArray("AliKFParticle");
1374 elec->SetOwner(kTRUE);
1375
1376 for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1377 {
1378 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1379 Int_t trackLabel = TMath::Abs(track->GetLabel());
1380
1381 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1382 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1383
1384 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1385 }
1386
1387 return elec;
1388}
1389
1390//--------------------------------------------------------------------------
1391TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1392{
7811d09e 1393//
1394// Find true gammas from pi0 Dalitz decay candidates with MC
1395//
7f3c7cc6 1396 TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1397 gammaPi0->SetOwner(kTRUE);
1398
1399 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1400 {
1401 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1402
1403 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1404 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1405
1406 if( !HaveSameMother(labelv1,labelv2) ) continue;
1407
1408 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1409
1410 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1411
1412 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1413
1414 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1415 }
1416
1417 return gammaPi0;
1418}
1419
1420//--------------------------------------------------------------------------
1421TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1422{
7811d09e 1423//
1424// Find true gammas from gamma candidates with MC
1425//
7f3c7cc6 1426 TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1427 gammaConv->SetOwner(kTRUE);
1428
1429 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1430 {
1431 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1432
1433 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1434 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1435
1436 if( !HaveSameMother(labelv1,labelv2) ) continue;
1437
1438 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1439
1440 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1441
1442 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1443 }
1444
1445 return gammaConv;
1446}
1447
1448//--------------------------------------------------------------
7f3c7cc6 1449void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1450{
7811d09e 1451//
1452// Remove repeated electron candidate tracks
1453// according to the gamma candidate array
1454//
7f3c7cc6 1455 vector<Bool_t> posTag(pos.size(),kTRUE);
1456 vector<Bool_t> negTag(neg.size(),kTRUE);
1457
1458 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1459 {
1460 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1461 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1462
1463 for( UInt_t j=0; j < pos.size(); ++j )
1464 {
1465 if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1466 }
1467
1468 for( UInt_t j=0; j < neg.size(); ++j )
1469 {
1470 if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1471 }
1472 }
1473
1474 CleanArray(pos, posTag);
1475 CleanArray(neg, negTag);
1476}
1477
1478//--------------------------------------------------------------------------
7f3c7cc6 1479void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1480{
7811d09e 1481//
1482// Remove electron candidates from gamma conversions
1483// according to the Psi pair angle
1484//
7f3c7cc6 1485 vector<Bool_t> posTag(pos.size(), kTRUE);
1486 vector<Bool_t> negTag(neg.size(), kTRUE);
1487
1488 for( UInt_t i=0; i < pos.size(); ++i )
1489 {
1490 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1491
1492 for( UInt_t j=0; j < neg.size(); ++j )
1493 {
1494 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1495
1496 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1497 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1498
1499 if(IsFromGammaConversion( psiPair, deltaPhi ))
1500 {
1501 posTag[i] = kFALSE;
1502 negTag[j] = kFALSE;
1503 }
1504 }
1505 }
1506
1507 CleanArray(pos, posTag);
1508 CleanArray(neg, negTag);
1509}
1510
1511//-----------------------------------------------------------------------------------
7f3c7cc6 1512void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1513{
7811d09e 1514//
1515// Remove electron candidates pairs
1516// which have mass not in the range (fMassCutMin,fMassCutMax)
1517//
7f3c7cc6 1518 vector<Bool_t> posTag(pos.size(), kTRUE);
1519 vector<Bool_t> negTag(neg.size(), kTRUE);
1520
1521 for( UInt_t i=0; i < pos.size(); ++i )
1522 {
1523 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1524
1525 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1526 TLorentzVector posLV;
b9295e41 1527 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1528
1529 for( UInt_t j=0; j < neg.size(); ++j )
1530 {
1531 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1532
1533 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1534 TLorentzVector negLV;
b9295e41 1535 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1536
1537 TLorentzVector posnegLV = posLV + negLV;
1538
1539 if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1540 {
1541 posTag[i] = kFALSE;
1542 negTag[j] = kFALSE;
1543 }
1544 }
1545 }
1546
1547 CleanArray(pos, posTag);
1548 CleanArray(neg, negTag);
1549}
1550
1551//-----------------------------------------------------------------------------------------------
1552void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1553{
7811d09e 1554//
1555// Clean the x array according to the tag parameter
1556//
7f3c7cc6 1557 vector<Int_t> tmp;
1558
1559 for(UInt_t i=0; i< x.size(); ++i)
1560 {
1561 if(tag[i]) tmp.push_back(x[i]);
1562 }
1563
1564 x = tmp;
1565}
1566
1567//--------------------------------------------------------------------------
7f3c7cc6 1568void 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)
1569{
7811d09e 1570//
1571// Remove gamma candidates according to
1572// the angle between the plane e+,e- and the gamma
1573//
7f3c7cc6 1574 vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1575
1576 for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1577 {
1578 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1579 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1580 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1581
1582 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1583 {
1584 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1585 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1586 TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1587
1588 // normal vector to x+y- plane
1589 TVector3 planePosNeg = xMom.Cross(yMom);
1590 for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1591 {
1592 AliKFParticle* gamKF = (AliKFParticle*)candidates->At(i);
1593 TVector3 gamMom(gamKF->Px(),gamKF->Py(),gamKF->Pz());
1594 if (planePosNeg.Angle(gamMom) < 1. || planePosNeg.Angle(gamMom) > 2.)
1595 {
1596 gammaTag[i] = kFALSE;
1597 }
1598 }
1599 }
1600 }
1601
1602 // Rebuild gamma candidates array
1603
1604 if(gamma) delete gamma;
1605 gamma = new TClonesArray("AliKFParticle");
1606 gamma->SetOwner(kTRUE);
1607
1608 posGamIdx.clear();
1609 negGamIdx.clear();
1610
1611 for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1612 {
1613 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1614 if(gammaTag[i])
1615 {
1616 new ((*gamma)[j++]) AliKFParticle(*iGamma);
1617 posGamIdx.push_back(fV0Reader->GetPindex(i));
1618 negGamIdx.push_back(fV0Reader->GetNindex(i));
1619 }
1620 }
1621}
1622
1623//--------------------------------------------------------------------------
7f3c7cc6 1624TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1625{
7811d09e 1626//
1627// Find Dalitz pair candidates
1628//
7f3c7cc6 1629 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1630 dalitz->SetOwner(kTRUE);
1631
1632 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1633 {
1634 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1635
1636 TLorentzVector posLV;
b9295e41 1637 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1638
1639 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1640 {
1641 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1642
1643 TLorentzVector negLV;
b9295e41 1644 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
7f3c7cc6 1645
1646 if(fUseAliKF)
1647 {
1648 AliKFParticle posNeg( *posKF, *negKF);
1649
1650 TLorentzVector posNegLV;
1651 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1652 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1653 }
1654 else
1655 {
1656 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1657 }
1658 }
1659 }
1660
1661 return dalitz;
1662}
1663
1664//--------------------------------------------------------------------------
2eedd4ed 1665TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleDalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma,Int_t opc)
7f3c7cc6 1666{
7811d09e 1667//
1668// Find pi0 Dalitz decay candidates
1669//
7f3c7cc6 1670 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1671 pi0->SetOwner(kTRUE);
1672
1673 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1674 {
1675 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1676
1677 TLorentzVector posLV;
b9295e41 1678 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1679
1680 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1681 {
1682 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1683
1684 TLorentzVector negLV;
b9295e41 1685 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
2eedd4ed 1686
1687 AliKFParticle posNegKF(*posKF,*negKF);
1688
7f3c7cc6 1689
1690 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1691 {
1692 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
2eedd4ed 1693 AliKFParticle posNegGam( *posKF, *negKF, *gamKF );
1694
1695 Double_t lDiffMass = posNegGam.GetMass() - posNegKF.GetMass();
1696
1697 if( opc == 1 )
1698 {
1699 fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_Diff",lDiffMass );
1700 fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1701 }
1702 else if ( opc == 2 )
1703 {
1704 fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_Diff",lDiffMass );
1705 fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1706 }
1707 else if ( opc == 3 )
1708 {
1709 fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_Diff",lDiffMass );
1710 fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1711 }
7f3c7cc6 1712
1713 if(fUseAliKF)
1714 {
2eedd4ed 1715
7f3c7cc6 1716 TLorentzVector posNegGamLV;
1717 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1718 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1719 }
1720 else
1721 {
1722 TLorentzVector gamLV;
1723 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1724
1725 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1726 }
1727 }
1728 }
1729 }
1730
1731 return pi0;
1732}
1733
1734//--------------------------------------------------------------------------
2eedd4ed 1735TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx,Int_t motherOpc)
7f3c7cc6 1736{
7811d09e 1737//
1738// Find true Dalitz pairs from Dalitz pair candidats with MC
1739//
7f3c7cc6 1740 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1741 dalitz->SetOwner(kTRUE);
1742
1743 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1744 {
1745 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1746 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1747 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1748
1749 TLorentzVector posLV;
b9295e41 1750 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1751
1752 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1753
1754 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1755 {
1756 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1757 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1758
2eedd4ed 1759 if(!IsDalitzPair(posLabel,negLabel,motherOpc)) continue;
7f3c7cc6 1760
1761 if(fUseAliKF)
1762 {
1763 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1764 AliKFParticle posNeg( posKF, negKF);
1765
1766 TLorentzVector posNegLV;
1767 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1768
1769 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1770 }
1771 else // TLorentzVector
1772 {
1773 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1774
1775 TLorentzVector negLV;
b9295e41 1776 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1777
1778 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1779 }
1780 }
1781 }
1782
1783 return dalitz;
1784}
1785
2eedd4ed 1786TClonesArray* AliAnalysisTaskGammaConvDalitz::FindJpsi(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx,Int_t motherOpc)
1787{
1788//
1789// Find true Jpsi's
1790// If mother
1791// -1: from the all sources
1792// 0: from the Chic_0
1793// 1: from the Chic_1
1794// 2: from the Chic_2
1795
1796 TClonesArray* jPsi = new TClonesArray("TLorentzVector");
1797 jPsi->SetOwner(kTRUE);
1798
1799 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1800 {
1801 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1802 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1803 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1804
1805 if( fStack->Particle(posLabel)->GetPdgCode() != ::kPositron ) continue;
1806
1807 TLorentzVector posLV;
1808 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1809
1810 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1811
1812 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1813 {
1814 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1815 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1816
1817 if( fStack->Particle(negLabel)->GetPdgCode() != ::kElectron ) continue;
1818
1819 if( !HaveSameMother(posLabel,negLabel) ) continue;
1820
1821 Int_t motherLabel = fStack->Particle(negLabel)->GetMother(0);
1822 TParticle *motherJpsi = fStack->Particle(motherLabel);
1823
1824 if( motherJpsi->GetPdgCode() != 443 ){
1825 continue;
1826 }
1827
1828
1829
1830 if( motherOpc > -1 )
1831 {
1832 if( motherJpsi->GetMother(0) < 0 ) continue;
1833
1834 TParticle *gmotherChic = fStack->Particle(motherJpsi->GetMother(0));
1835 Int_t pdgCode = gmotherChic->GetPdgCode();
1836
1837 Bool_t lson = kTRUE;
1838
1839 switch(motherOpc){
1840
1841 case 0: if ( pdgCode != 10441 )
1842 lson = kFALSE;
1843 break;
1844 case 1: if ( pdgCode != 20443 )
1845 lson = kFALSE;
1846 break;
1847 case 2: if ( pdgCode != 445 )
1848 lson = kFALSE;
1849 break;
1850 }
1851
1852 if( lson == kFALSE ) continue;
1853 }
1854
1855
1856 if(fUseAliKF)
1857 {
1858 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1859 AliKFParticle posNeg( posKF, negKF);
1860
1861 TLorentzVector posNegLV;
1862 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1863
1864 new ((*jPsi)[j++]) TLorentzVector(posNegLV);
1865 }
1866 else // TLorentzVector
1867 {
1868 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1869
1870 TLorentzVector negLV;
1871 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1872
1873 new ((*jPsi)[j++]) TLorentzVector(posLV + negLV);
1874 }
1875 }
1876 }
1877
1878 return jPsi;
1879}
1880
1881
1882
7f3c7cc6 1883//--------------------------------------------------------------------------
2eedd4ed 1884TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleDalitz(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam,Int_t motherOpc)
7f3c7cc6 1885{
7811d09e 1886//
1887// Find true pi0 Dalitz decay from pi0 candidates with MC
1888//
7f3c7cc6 1889 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1890 pi0->SetOwner(kTRUE);
1891
1892 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1893 {
1894 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1895 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1896 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1897
1898 TLorentzVector posLV;
b9295e41 1899 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1900
1901 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1902
1903 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1904 {
1905 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1906 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
2eedd4ed 1907
1908
1909 if( !HaveSameMother(posLabel,negLabel) ) continue; //Check if both particles have same mother
1910 if(!IsDalitzPair(posLabel,negLabel,motherOpc)) continue; //check if mohter is eta0 or pi0
1911
7f3c7cc6 1912 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1913
1914 TLorentzVector negLV;
b9295e41 1915 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1916
1917 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1918
1919 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1920 {
1921 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1922
1923 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1924 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1925
1926 if( !HaveSameMother(labelv1,labelv2) ) continue;
1927
1928 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1929
1930 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
2eedd4ed 1931
1932
7f3c7cc6 1933 if( !HaveSameMother(labelGamma, posLabel) ) continue;
2eedd4ed 1934
7f3c7cc6 1935
1936 if(fUseAliKF)
1937 {
1938 AliKFParticle posNegGam( posKF, negKF, *gamKF );
1939 TLorentzVector posNegGamLV;
1940 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1941 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1942 }
1943 else // TLorentzVector
1944 {
1945 TLorentzVector gamLV;
1946 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1947
1948 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1949 }
1950 }
1951 }
1952 }
1953
1954 return pi0;
1955}
2eedd4ed 1956TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleChic(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam,Int_t motherOpc)
1957{
1958//
1959// Find true pi0 Dalitz decay from pi0 candidates with MC
1960//
1961 TClonesArray* chic = new TClonesArray("TLorentzVector");
1962 chic->SetOwner(kTRUE);
1963
1964 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1965 {
1966 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1967 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1968 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1969
1970 if( fStack->Particle(posLabel)->GetPdgCode() != ::kPositron ) continue;
1971
1972 TLorentzVector posLV;
1973 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1974
1975 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1976
1977 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1978 {
1979 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1980 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1981
1982
1983 if( fStack->Particle(negLabel)->GetPdgCode() != ::kElectron ) continue;
1984
1985
1986 if( !HaveSameMother(posLabel,negLabel) ) continue;
1987
1988 Int_t jpsiLabel = fStack->Particle(negLabel)->GetMother(0);
1989
1990 if( fStack->Particle(jpsiLabel)->GetPdgCode() != 443 ) continue;
1991
1992 TParticle *jpsiParticle = fStack->Particle(jpsiLabel);
1993
1994 if ( jpsiParticle->GetMother(0) < 0 ) continue;
1995
1996
1997 Int_t chicLabel = jpsiParticle->GetMother(0);
1998
1999 Int_t pdgCode = fStack->Particle(chicLabel)->GetPdgCode();
2000
2001
2002 Bool_t lSon = kTRUE;
2003
2004 switch(motherOpc){
2005
2006 case 0: if ( pdgCode != 10441 )
2007 lSon = kFALSE;
2008 break;
2009 case 1: if ( pdgCode != 20443 )
2010 lSon = kFALSE;
2011 break;
2012 case 2: if ( pdgCode != 445 )
2013 lSon = kFALSE;
2014 break;
2015 }
2016
2017
2018 if( lSon == kFALSE ) continue;
2019
2020
2021
2022
2023
2024 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
2025
2026 TLorentzVector negLV;
2027 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
2028
2029 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
2030
2031 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
2032 {
2033 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
2034
2035 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
2036 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
2037
2038 if( !HaveSameMother(labelv1,labelv2) ) continue;
2039
2040 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
2041
2042 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
2043
2044
2045 if( !HaveSameMother(labelGamma, jpsiLabel) ) continue;
2046
2047
2048 if(fUseAliKF)
2049 {
2050 AliKFParticle posNegGam( posKF, negKF, *gamKF );
2051 TLorentzVector posNegGamLV;
2052 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
2053 new ((*chic)[j++]) TLorentzVector(posNegGamLV);
2054 }
2055 else // TLorentzVector
2056 {
2057 TLorentzVector gamLV;
2058 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
2059
2060 new ((*chic)[j++]) TLorentzVector(posLV + negLV + gamLV);
2061 }
2062 }
2063 }
2064 }
2065
2066 return chic;
2067}
7f3c7cc6 2068
2069//-----------------------------------------------------------------------------------------------
2070void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
2071{
7811d09e 2072//
2073// Update gamma event pool for background computation
2074//
7f3c7cc6 2075 if( fDebug ) AliInfo("=> UpdateGammaPool");
2076
2077 // cycle
2078 for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
2079 {
2080 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
2081 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
2082 ++fGamPoolPos;
2083 if(fGamPoolPos == fPoolMaxSize)
2084 {
2085 fGamPoolPos = 0;
2086 }
2087 }
2088}
2089
2090void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
2091{
7811d09e 2092//
2093// Update electron event pool for background computation
2094//
2095 Int_t multiplicity = fV0Reader->CountESDTracks();
2096
2097
2098 fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
7f3c7cc6 2099}
2100
2101//-----------------------------------------------------------------------------------------------
7f3c7cc6 2102TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
2103{
7811d09e 2104//
2105// Gamma copy from events with same multiplicity and Z
2106//
7f3c7cc6 2107 if( fDebug ) AliInfo("=> GammasFromBGHandler");
7811d09e 2108
2109 Int_t zbin = fBGEventHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2110 Int_t mbin = fBGEventHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2111
2112
7f3c7cc6 2113
2114 TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
2115 gammaPool->SetOwner(kTRUE);
2116
2117 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
2118 {
7811d09e 2119 AliGammaConversionKFVector* gammaV0s = fBGEventHandler->GetBGGoodV0s(zbin,mbin,iEventBG);
7f3c7cc6 2120 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
2121 {
2122 new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
2123 }
2124 }
2125
2126 return gammaPool;
2127}
2128
2129//-----------------------------------------------------------------------------------------------
7f3c7cc6 2130TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
2131{
7811d09e 2132//
2133// Electron copy from events with same multiplicity and Z
2134//
7f3c7cc6 2135 if( fDebug ) AliInfo("=> ElectronFromBGHandler");
2136
2137 TClonesArray* electronPool = new TClonesArray("AliKFParticle");
2138 electronPool->SetOwner(kTRUE);
7811d09e 2139
2140 Int_t multiplicity = fV0Reader->CountESDTracks();
2141
2142
2143
7f3c7cc6 2144
2145 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
2146 {
7811d09e 2147 AliGammaConversionKFVector* electronNeg = fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
7f3c7cc6 2148 for (UInt_t i = 0; i < electronNeg->size(); ++i )
2149 {
2150 new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
2151 }
2152 }
2153
2154 return electronPool;
2155}
2156
2157//-----------------------------------------------------------------------------------------------
2158Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
2159{
7811d09e 2160//
2161// Get track pid according to MC
2162//
7f3c7cc6 2163 Int_t label = TMath::Abs(t->GetLabel());
2164 Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
2165
2166 switch(pdgCode)
2167 {
2168 case ::kElectron: return AliPID::kElectron;
2169 case ::kMuonMinus: return AliPID::kMuon;
2170 case ::kPiPlus: return AliPID::kPion;
2171 case ::kKPlus: return AliPID::kKaon;
2172 case ::kProton: return AliPID::kProton;
2173 }
2174
2175 return -1;
2176}
2177
2178//-----------------------------------------------------------------------------------------------
7811d09e 2179//FIXME PID ITS
2180// NOTE prior should be estimated from data
2181// NOTE: move to config
2182
7f3c7cc6 2183Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
2184{
7811d09e 2185//
2186// Get track pid according to Bayes' formula
2187//
7f3c7cc6 2188 double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
2189 Double_t detectoProb[AliPID::kSPECIES];
2190
2191 if( trackType == kITSsaTrack ) // ITS standalone pid
2192 {
2193 t->GetITSpid( detectoProb );
2194 }
2195 else // global track
2196 {
2197 t->GetESDpid( detectoProb );
2198 }
2199
2200 AliPID bayesPID( detectoProb );
2201 return bayesPID.GetMostProbable( priors );
2202}
2203
2204//-----------------------------------------------------------------------------------------------
2205Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
2206{
7811d09e 2207//
2208// Get track pid according to a n-sigma cut around ITS and/or TPC signals
2209//
7f3c7cc6 2210 if( trackType == kITSsaTrack) // ITS standalone tracks
2211 {
2212 Double_t mom = t->GetP();
2213
2214
2215 // ITS Number of sigmas (FIXME: add new fESDpidCuts)
2216 // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
2217 Double_t nElecSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kElectron );
2218 Double_t nPionSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kPion );
2219
2220 if( nElecSigma < 4. && nElecSigma > -3. && mom < .2 && nPionSigma < -1.5 )
2221 {
2222 return AliPID::kElectron;
2223 }
2224 }
2225 else // global track
2226 {
2227 Double_t nElecSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kElectron );
2228 Double_t nPionSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kPion );
2229 Double_t nKaonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kKaon ));
2230 Double_t nProtonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kProton));
2231 if( nElecSigma > fNSigmaBelowElecTPCbethe && nElecSigma < fNSigmaAboveElecTPCbethe &&
2232 nPionSigma > fNSigmaAbovePionTPCbethe && //NOTE mom > 0.5
2233 nKaonSigma > fNSigmaAboveKaonTPCbethe &&
2234 nProtonSigma > fNSigmaAboveProtonTPCbethe )
2235 {
2236 return AliPID::kElectron;
2237 }
2238 // NOTE: add other particle types
2239 }
2240
2241 return -1;
2242}
2243
2244//-----------------------------------------------------------------------------------------------
2eedd4ed 2245Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel,Int_t motherOpc ) const
7f3c7cc6 2246{
7811d09e 2247//
2248// Returns true if the two particles is a Dalitz pair
2249//
2eedd4ed 2250//motherOpc 1: for Pi0Dalitz
2251//motherOpc 2: for EtaDalitz
2252
7f3c7cc6 2253 if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
2254
2255 TParticle* pos = fStack->Particle( posLabel );
2256 TParticle* neg = fStack->Particle( negLabel );
2257
2258 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
2259 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
2260
2261 //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
2262 //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
2263
2264 Int_t motherLabel = pos->GetMother(0);
2265 if( motherLabel < 0 ) return kFALSE;
2266
2267 TParticle* mother = fStack->Particle( motherLabel );
2eedd4ed 2268
2269 if( mother->GetNDaughters() != 3) return kFALSE;
7f3c7cc6 2270
2eedd4ed 2271 if( motherOpc == 1 ){ //Pi0Dalitz
2272 if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
2273 }
2274 else if(motherOpc == 2){
2275 if( mother->GetPdgCode() != 221 ) return kFALSE;
2276 }
2277 else {
2278 return kFALSE;
2279 }
7f3c7cc6 2280 // NOTE: one of them must be a photon
2281
2282 return kTRUE;
2283}
2284
2285//-----------------------------------------------------------------------------------------------
7f3c7cc6 2286Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
2287{
7811d09e 2288//
2289// Returns true if the particle comes from Pi0 -> e+ e- gamma
2290//
7f3c7cc6 2291 Bool_t ePlusFlag = kFALSE;
2292 Bool_t eMinusFlag = kFALSE;
2293 Bool_t gammaFlag = kFALSE;
2294
2295 Int_t motherLabel = fStack->Particle( label )->GetMother(0);
2296
2297 if( motherLabel < 0 ) return kFALSE;
2298
2299 TParticle* mother = fStack->Particle( motherLabel );
2300
2301 if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
2302
2303 if ( mother->GetNDaughters() != 3 ) return kFALSE;
2304
2305 for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
2306 {
2307 switch( fStack->Particle(idx)->GetPdgCode())
2308 {
2309 case ::kPositron:
2310 ePlusFlag = kTRUE;
2311 break;
2312 case ::kElectron:
2313 eMinusFlag = kTRUE;
2314 break;
2315 case ::kGamma:
2316 gammaFlag = kTRUE;
2317 break;
2318 }
2319 }
2320
2321 return ( ePlusFlag && eMinusFlag && gammaFlag );
2322}
2323
2324//--------------------------------------------------------------------------
2325Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
2326{
7811d09e 2327//
2328// Returns true if it is a gamma conversion according to psi pair value
2329//
7f3c7cc6 2330 return ( (deltaPhi > fDeltaPhiCutMin && deltaPhi < fDeltaPhiCutMax) &&
2331 TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
2332}
2333
2334//--------------------------------------------------------------------------
2335Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
2336{
7811d09e 2337//
2338// Returns true if it is a gamma conversion according to MC
2339//
7f3c7cc6 2340 if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
2341
2342 TParticle* pos = fStack->Particle( posLabel );
2343 TParticle* neg = fStack->Particle( negLabel );
2344
2345 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
2346 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
2347
2348 if( pos->GetUniqueID() != kPPair ) return kFALSE;
2349
2350 Int_t motherLabel = pos->GetMother(0);
2351 if( motherLabel < 0 ) return kFALSE;
2352
2353 TParticle* mother = fStack->Particle( motherLabel );
2354
2355 if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
2356
2357 return kTRUE;
2358}
2359
2360//-----------------------------------------------------------------------------------------------
2361Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
2362{
7811d09e 2363//
2364// Returns true if the two particle have the same mother
2365//
7f3c7cc6 2366 if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
2367 return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
2368}
2369
2370//-----------------------------------------------------------------------------------------------
2371Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
2372{
7811d09e 2373//
2374// This angle is a measure for the contribution of the opening in polar
2375// direction Δ0 to the opening angle ξ Pair
2376//
2377// Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2378// Master Thesis. Thorsten Dahms. 2005
2379// https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2380//
b9295e41 2381 Double_t momPos[3];
2382 Double_t momNeg[3];
2383 if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
2384 if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
7f3c7cc6 2385
b9295e41 2386 TVector3 posDaughter;
2387 TVector3 negDaughter;
7f3c7cc6 2388
b9295e41 2389 posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
2390 negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
7f3c7cc6 2391
b9295e41 2392 Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
2393 Double_t openingAngle = posDaughter.Angle( negDaughter ); //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
7f3c7cc6 2394 if( openingAngle < 1e-20 ) return 0.;
2395 Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
2396
2397 return psiAngle;
2398}
2399
2400//-----------------------------------------------------------------------------------------------
7811d09e 2401Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const AliKFParticle* xPos, const AliKFParticle* yNeg ) const
7f3c7cc6 2402{
7811d09e 2403//
2404// Get psi pair value
2405//
b9295e41 2406 TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
2407 TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
7f3c7cc6 2408
2409 Double_t deltaTheta = neg.Theta() - pos.Theta();
2410 Double_t openingAngle = pos.Angle( neg );
2411
2412 if( openingAngle < 1e-20 ) return 0.;
2413
2414 return TMath::ASin( deltaTheta/openingAngle );
2415}
2416
2417//-----------------------------------------------------------------------------------------------
2418Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
2419{
7811d09e 2420//
2421// Get psi pair value
2422//
7f3c7cc6 2423 Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
2424 Double_t openingAngle = xPos->Angle( yNeg->Vect() );
2425
2426 if( openingAngle < 1e-20 ) return 0.;
2427
2428 return TMath::ASin( deltaTheta/openingAngle );;
2429}
2430
2431// tmp NOTE: Should go to AliV0Reader
2432//-----------------------------------------------------------------------------------------------
2433/*
2434Double_t AliAnalysisTaskGammaConvDalitz::GetPsiAngleV0(
2435 Double_t radiusVO, // radius at XY of VO vertex
2436 const AliExternalTrackParam* trackPos, // pos. track parm. at V0 vertex
2437 const AliExternalTrackParam* trackNeg // neg. track parm. at V0 vertex
2438 )
2439{
2440 // This angle is a measure for the contribution of the opening in polar
2441 // direction Δ0 to the opening angle ξ Pair
2442
2443 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2444 // Master Thesis. Thorsten Dahms. 2005
2445 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2446
2447 const Double_t kForceDeltaPhi = 1.2;
2448 static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
2449
2450 Double_t psiAngle = 0.; // Default value
2451
2452 static Double_t MagnField = fESDEvent->GetMagneticField();
2453 static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
2454
2455 if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
2456
2457 // compute propagation radius for a fixed angle
2458 Double_t Rpos = radiusVO + trackPos->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2459 Double_t Rneg = radiusVO + trackNeg->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2460
7811d09e 2461 Double_t MomPos[3];
2462 Double_t MomNeg[3];
2463 if( trackPos->GetPxPyPzAt( Rpos, MagnField, MomPos ) == 0 ) trackPos->GetPxPyPz( MomPos );
2464 if( trackNeg->GetPxPyPzAt( Rneg, MagnField, MomNeg ) == 0 ) trackNeg->GetPxPyPz( MomNeg );
7f3c7cc6 2465
7811d09e 2466 TVector3 PosDaughter;
2467 TVector3 NegDaughter;
2468 PosDaughter.SetXYZ( MomPos[0], MomPos[1], MomPos[2] );
2469 NegDaughter.SetXYZ( MomNeg[0], MomNeg[1], MomNeg[2] );
7f3c7cc6 2470
7811d09e 2471 Double_t deltaTheta = NegDaughter.Theta() - PosDaughter.Theta();
2472 Double_t chiPar = PosDaughter.Angle( NegDaughter ); //TMath::ACos( PosDaughter.Dot(NegDaughter) / (NegDaughter.Mag()*PosDaughter.Mag()) );
7f3c7cc6 2473 if( chiPar < 1e-20 ) return psiAngle;
2474
2475 psiAngle = TMath::ASin( deltaTheta / chiPar );
2476
2477 return psiAngle;
2478}
2479*/