]> 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
Changes to fast embedding and jetresponse: QA plot for delta AOD, Jet eta phi selecti...
[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//
7f3c7cc6 238 if( fDebug ) AliInfo("=> UserExec");
239
240 if( fV0Reader == 0 )
241 {
242 AliFatal("no pointer to AliV0Reader");
243 return;
244 }
245
7f3c7cc6 246 // Create list of gamma candidates in standalone mode
247 // otherwise use the created ones by AliAnalysisTaskGammaConversion
248 if( fStandalone )
249 {
7811d09e 250
251 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
252 AliESDInputHandler *esdHandler=0;
253 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
254 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
255 } else {
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()){
260 if (fMCEvent ) {
261 AliV0Reader::InitESDpid();
262 } else {
263 AliV0Reader::InitESDpid(1);
264 }
265 }
266 }
267
268
269 if (fMCEvent ) {
270
271 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
272 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
273
274 if (!mcHandler){
275 AliError("Could not retrive MC event handler!");
276 return;
277 }
278
279 if (!mcHandler->InitOk() ){
280
281 return;
282 }
283 if (!mcHandler->TreeK() ){
284 return;
285 }
286 if (!mcHandler->TreeTR() ) {
287 return;
288 }
289 }
290
291
292
293
294 fV0Reader->SetInputAndMCEvent( InputEvent(), MCEvent() );
7f3c7cc6 295 fV0Reader->Initialize();
296 }
297
298 if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
299 {
300 if( fDebug ) AliInfo("no contributors to primary vertex");
301 return;
302 }
303
7811d09e 304 if( fV0Reader->CheckForPrimaryVertexZ() == kFALSE )
305 {
306
307 if( fDebug ) AliInfo("z vertex out of range");
308 return;
309 }
310
7f3c7cc6 311 // Get Pointers
312 fBGEventHandler = fV0Reader->GetBGHandler();
313 fESDpid = fV0Reader->GetESDpid();
314 fESDEvent = fV0Reader->GetESDEvent();
315 if(fDoMC)
316 {
7811d09e 317 fStack= MCEvent()->Stack();
318 fGCMCEvent=MCEvent();
7f3c7cc6 319 }
320
321 // Read the magnetic field sign from ESD
322 if ( fReadMagFieldSign == kTRUE )
323 {
324 fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
325 }
326
327 // Process MC information
328 if(fDoMC)
329 {
330 ProcessMCData();
331 }
332
333 if( fStandalone ){
334 while(fV0Reader->NextV0()){}; //SelectGammas
335 fV0Reader->ResetV0IndexNumber();
336 }
7811d09e 337
7f3c7cc6 338 CreateListOfDalitzPairCandidates();
339 ProcessGammaElectronsForDalitzAnalysis();
7811d09e 340
341 if ( fStandalone ){
342
343 fV0Reader->UpdateEventByEventData();
344
345 }
7f3c7cc6 346
347 PostData( 1, fOutputContainer );
348}
349
350
351void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
352{
7811d09e 353//
7f3c7cc6 354 if( fDebug ) AliInfo("Not to do anything in Terminate");
355}
356
357//-----------------------------------------------------------------------------------------------
7f3c7cc6 358void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
359{
7811d09e 360//
361// set user ITSsa track cuts
362//
7f3c7cc6 363 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
364
365 if( esdCuts )
366 {
367 fITSsaTrackCuts = esdCuts;
368 }
369 else
370 {
7811d09e 371 // default cuts
7f3c7cc6 372 fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
373
374 fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
375 fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
376
377 fITSsaTrackCuts->SetMinNClustersITS(2);
378 fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
379 fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
380
381 fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
382 fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
383
384 fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
385 }
386}
387
388//-----------------------------------------------------------------------------------------------
7f3c7cc6 389void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
390{
7811d09e 391//
392// set user global track cuts
393//
7f3c7cc6 394 if( fESDtrackCuts ) delete fESDtrackCuts;
395
396 if( esdCuts )
397 {
398 fESDtrackCuts = esdCuts;
399 }
400 else
401 {
7811d09e 402 //default cuts
7f3c7cc6 403 fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
404
405 fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
406 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
407
408 fESDtrackCuts->SetMinNClustersITS(2);
409 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
410 fESDtrackCuts->SetRequireITSRefit(kTRUE);
411
412 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
413 fESDtrackCuts->SetMaxNsigmaToVertex(3);
414
415 fESDtrackCuts->SetMinNClustersTPC(80);
416 fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
417 fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
418 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
419 }
420}
421
422//-----------------------------------------------------------------------------------------------
423void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
424{
7811d09e 425//
426// set user pid cuts
427//
7f3c7cc6 428 if( fESDpidCuts ) delete fESDpidCuts;
429 if( esdPIDCuts )
430 {
431 fESDpidCuts = esdPIDCuts;
432 }
7811d09e 433 else // default cuts
7f3c7cc6 434 {
435 fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
436 // fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
437 fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
438 }
439}
440
441//-----------------------------------------------------------------------------------------------
7f3c7cc6 442void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
443{
7811d09e 444//
445// Process generation
446//
7f3c7cc6 447 if( fDebug ) AliInfo("=> ProcessMCData");
448
449 fHistograms->FillTable("Table_Generation", 0); //number of events
450
451 for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
452 {
453 TParticle* iParticle = fStack->Particle( i );
454 if( !iParticle ) continue;
455
456 if ( i >= fStack->GetNprimary() ) continue; // is primary?
457 if ( iParticle->GetPdgCode() != ::kPi0 ) continue; // is Pi0?
458
459 if( iParticle->GetNDaughters() == 2 &&
460 fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
461 fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
462 {
463 fHistograms->FillTable("Table_Generation", 1); // pi0 -> gg
464 }
465
466 if ( iParticle->GetNDaughters() != 3 ) continue; // Num == 3 (e+,e-,gamma)
467
468 // Check for Pi0 Dalitz decay
469 TParticle* eposPi0 = 0;
470 TParticle* enegPi0 = 0;
471 TParticle* gammaPi0 = 0;
472
473 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
474 {
475 switch(fStack->Particle(idxPi0)->GetPdgCode())
476 {
477 case ::kPositron:
478 eposPi0 = fStack->Particle(idxPi0);
479 break;
480 case ::kElectron:
481 enegPi0 = fStack->Particle(idxPi0);
482 break;
483 case ::kGamma:
484 gammaPi0 = fStack->Particle(idxPi0);
485 break;
486 }
487 }
488
489 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
490
491 // found a Pi0 Dalitz decay
492
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());
504
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);
511
512 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
513
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());
517
518
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;
524
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 )
528 {
529 for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
530 {
531 TParticle* tmpDaughter = fStack->Particle(tIndex);
7811d09e 532
7f3c7cc6 533 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
534 if( tmpDaughter->GetPdgCode() == ::kElectron )
535 { // e+
536 daugGammaElectronAll = kTRUE;
7811d09e 537
7f3c7cc6 538 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
539 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
540 (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
541 {
542 daugGammaElectron = kTRUE;
543 }
544 }
545 else if( tmpDaughter->GetPdgCode() == ::kPositron )
546 {
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() ) )
551 {
552 daugGammaPositron = kTRUE;
553 }
554 }
555 }
556 }
557
558
559 if( daugGammaElectronAll && daugGammaPositronAll )
560 {
561 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
562 fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
563 }
564
565 Float_t etaMin, etaMax;
566 fESDtrackCuts->GetEtaRange( etaMin, etaMax );
567
568 // e+e- pair within acceptance
569 if ( TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
570 {
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());
576 }
577
578 // Pi0 (e+e-gamma) within acceptance
579 //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
580
581 if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
582 TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
583 {
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());
7811d09e 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());
7f3c7cc6 594 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
595 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
596
597 // Pi0 within acceptance with gamma converted
598
599 if ( daugGammaElectron && daugGammaPositron )
600 {
601 fHistograms->FillTable("Table_Generation",4); //
602
603 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
604 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
7811d09e 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());
7f3c7cc6 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());
616 }
617 }
618 }
619}
620
621//-----------------------------------------------------------------------------------------------
7f3c7cc6 622void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
623{
7811d09e 624//
625// Dalitz pair candidates
626//
7f3c7cc6 627 if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
628
629 fEposCandidateIndex.clear();
630 fEnegCandidateIndex.clear();
631
632 fHistograms->FillTable("Table_Cuts", 0);
633
634 for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
635 {
636 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
637 if ( !iTrack ) continue;
638
639
640 Double_t p[3];
641
642 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
643
644 TVector3 iMom(p[0],p[1],p[2]);
645
646 //
647 // Check track cuts and find track type
648 //
649
650 Bool_t isTrackAccepted = 0;
651 Int_t trackType = -1;
652 switch(fTrkSelectionCriteria)
653 {
654 case kITSsaTrack:
655 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
656 trackType = kITSsaTrack;
657 break;
658
659 case kGlobalTrack:
660 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
661 trackType = kGlobalTrack;
662 break;
663
664 case kITSsaGlobalTrack:
665 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
666 {
667 isTrackAccepted = kTRUE;
668 if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
669 else trackType = kGlobalTrack;
670 }
671 break;
672 }
673
674 if(!isTrackAccepted) continue;
675
676 //
677 // PID
678 //
679
680 Int_t pid=-1;
681 Int_t pidMC=-1;
682
683 if(fUseBayesPID)
684 {
685 pid = GetBayesPid(iTrack,trackType);
686 }
687 else
688 {
689 pid = GetNSigmaPid(iTrack,trackType);
690 }
691
692 if( fDoMC )
693 {
694 pidMC = GetMonteCarloPid(iTrack);
695 // pid table
696 Int_t iLabel = TMath::Abs(iTrack->GetLabel());
697 TParticle* iParticle = fStack->Particle(iLabel);
698 FillPidTable(iParticle, pid);
699 }
7811d09e 700
7f3c7cc6 701 // ITS standalone tracks
702 if( trackType == kITSsaTrack)
703 {
704 Double_t mom = iTrack->GetP();
705 Double_t signal = iTrack->GetITSsignal();
706
707 fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
708
709 if( pid == AliPID::kElectron )
710 {
711
712 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
713 if(fDoMC && pid == pidMC)
714 {
715 fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
716 }
717 }
718
719 if( fDoMC && pidMC == AliPID::kElectron)
720 {
721 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
722 }
723 }
724
725 else // global tracks
726 {
727 const AliExternalTrackParam *in = iTrack->GetInnerParam();
728 Double_t mom = in->GetP();
729 Double_t signal = iTrack->GetTPCsignal();
730
731 fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
732
733 if( fDoMC && pidMC == AliPID::kElectron )
734 {
735 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
736 }
737
738 if( pid == AliPID::kElectron )
739 {
740 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
741 if(fDoMC && pid == pidMC)
742 {
743 fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
744 }
745 }
746 }
747
748 if( AliPID::kElectron != pid) continue;
749
750 // electron track candidates from here
751
752 if( iTrack->GetSign() > 0 )
753 {
754 fEposCandidateIndex.push_back(i);
755 }
756 else
757 {
758 fEnegCandidateIndex.push_back(i);
759 }
760 }
761
762 // gamma candidates
763 GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
764
765 if(fDoMC)
766 {
767 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
768 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
769 delete pi0Dalitz;
770 }
771
772 if(fUseTrackIndexCut) // remove repeated tracks
773 {
774 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
775
776 if(fDoMC)
777 {
778 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
779 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
780 delete pi0Dalitz;
781 }
782 }
783
784 if(fUsePsiPairCut) // remove electrons from gamma conversions
785 {
786 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
787
788 if(fDoMC)
789 {
790 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
791 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
792 delete pi0Dalitz;
793 }
794 }
795
796 if( fUseMassCut )
797 {
798 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
799
800 if(fDoMC)
801 {
802 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
803 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
804 delete pi0Dalitz;
805 }
806 }
807
808 if(fUseGammaCut)
809 {
810 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
811
812 if(fDoMC)
813 {
814 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
815 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
816 delete pi0Dalitz;
817 }
818 }
819}
820
821//-----------------------------------------------------------------------------------------------
822void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
823{
7811d09e 824//
825// Process gamma and electrons for pi0 Dalitz decay
826//
7f3c7cc6 827 if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
828
829 fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
830
831 TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
832
833 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
834 {
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);
839 }
840
841 TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
842
843 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
844 {
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);
849 }
850
851 TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
852 for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
853 {
854 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
855
856 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
857 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Mass", dalitz->M());
858 }
859
860 // gamma candidates
861 for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
862 {
863 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
864 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
865 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
866 }
867
868 // psi pair for all candidates
869 //if(fUsePsiPairCut)
870 FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
871
872 // Angle epos,eneg gamma
873 FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
874 FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
875
876 TClonesArray* pi0Candidates = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaCandidates);
877 for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
878 {
879 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
880
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());
891 }
892
893 delete dalitzPairCandidates;
894 delete pi0Candidates;
895
896 if(fComputeBkg)
897 {
898 // 1) e+e- with with gammas used in the signal
899 TClonesArray* pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaPool);
900
901 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
902 {
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());
906 }
907
908 if(ePosCandidates->GetEntriesFast() > 0 &&
909 eNegCandidates->GetEntriesFast() > 0 &&
910 fGammaCandidates->GetEntriesFast() > 0)
911 {
912 UpdateGammaPool(fGammaCandidates);
913 }
914
915 delete pi0Bkg;
916
917 // 2) e+e- with gammas from a pool of events
7811d09e 918 TClonesArray* gammaBGHandler = GammasFromBGHandler();
7f3c7cc6 919 pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, gammaBGHandler);
920
921 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
922 {
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());
926 }
927
928 delete pi0Bkg;
929
930 // 3) e+ with e-, gamma from a pool of events
931 TClonesArray* elecBGHandler = ElectronFromBGHandler();
932 pi0Bkg = FindPi0Dalitz(ePosCandidates, elecBGHandler, gammaBGHandler);
933
934 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
935 {
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());
939 }
940
941 if(eNegCandidates->GetEntriesFast() > 0)
942 {
943 UpdateElectronPool(eNegCandidates);
944 }
945
946 delete gammaBGHandler;
947 delete elecBGHandler;
948 delete pi0Bkg;
7811d09e 949
7f3c7cc6 950 }
951
952 delete ePosCandidates;
953 delete eNegCandidates;
954
955 if(fDoMC)
956 {
957 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
958 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
959 {
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);
964 }
965
966 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
967 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
968 {
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);
973 }
974
975 TClonesArray* dalitzPair = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex);
976 for(Int_t i=0; i < dalitzPair->GetEntriesFast(); ++i)
977 {
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 );
982 }
983
984 // psi pair for dalitz pairs
985 //if(fUsePsiPairCut)
986 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_DalitzPair_PsiPair_vs_DPhi");
987
988 delete ePosPi0Dalitz;
989 delete eNegPi0Dalitz;
990 delete dalitzPair;
991
992 // all gammas
993 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
994 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
995 {
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());
999 }
1000
1001 delete gamma;
1002
1003 // gamma from pi0 dalitz
1004 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1005 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
1006 {
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);
1011 }
1012
1013 delete gammaPi0Dalitz;
1014
1015 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1016
1017 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
1018 {
1019 TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
1020
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);
1032 }
1033 delete pi0Dalitz;
1034
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)
1038 {
1039 AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
1040 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1041
1042 for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
1043 {
1044 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
1045 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1046
1047 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
1048
1049 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1050 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1051
1052 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
1053 }
1054 }
1055 // FIXME: eta -> e+e-gamma
1056 }
1057}
1058
1059//--------------------------------------------------------------------------
1060Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
1061{
7811d09e 1062//
1063// Get rapidity
1064//
1065 const double kEPSILON=1.e-16;
7f3c7cc6 1066
7811d09e 1067 if(p->Energy() - TMath::Abs(p->Pz()) < kEPSILON )
7f3c7cc6 1068 {
1069 return 1.e10;
1070 }
1071 return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
1072}
1073
1074//--------------------------------------------------------------------------
1075void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1076{
7811d09e 1077//
1078// Fill histogram with psipair(pos,neg)
1079//
7f3c7cc6 1080 for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1081 {
1082 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1083 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1084 {
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);
1089 }
1090 }
1091}
1092
1093//--------------------------------------------------------------------------
1094void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1095{
7811d09e 1096//
1097// Fill histogram with angle(x,y)
1098//
7f3c7cc6 1099 for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1100 {
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 )
1104 {
1105 AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1106 TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1107 fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1108 }
1109 }
1110}
1111
1112//--------------------------------------------------------------------------
7f3c7cc6 1113void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1114{
7811d09e 1115//
1116// Fill table with pid info
1117//
7f3c7cc6 1118 Int_t iGen=-1;
1119 switch(TMath::Abs(p->GetPdgCode()))
1120 {
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;
1126 default: iGen=-1;
1127 }
1128
1129 int jRec=-1;
1130 if(pid > -1 && pid < 5) jRec = pid;
1131
1132 if ((iGen > -1) && (jRec > -1))
1133 {
1134 fHistograms->FillTable("Table_PID", iGen, jRec);
1135 // sum
1136 fHistograms->FillTable("Table_PID", iGen, 5);
1137 fHistograms->FillTable("Table_PID", 5, jRec);
1138 }
1139}
1140
1141//--------------------------------------------------------------------------
1142void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1143{
7811d09e 1144//
1145// Make a copy of gamma candidates from V0reader
1146//
7f3c7cc6 1147 posIndex.clear();
1148 negIndex.clear();
1149
1150 if(gamma) delete gamma;
1151
1152 TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1153
1154 gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1155 gamma->SetOwner(kTRUE);
1156
1157 // make a copy
1158 for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1159 {
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));
1164 }
1165}
1166
1167//--------------------------------------------------------------------------
1168TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1169{
7811d09e 1170//
1171// Convert track index vector to AliKFParticle array
1172//
7f3c7cc6 1173 TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1174 indexKF->SetOwner(kTRUE);
1175
1176 for(UInt_t i = 0; i < index.size(); ++i)
1177 {
1178 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1179 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1180 }
1181
1182 return indexKF;
1183}
1184
1185//--------------------------------------------------------------------------
1186TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1187{
7811d09e 1188//
1189// Find true electrons from pi0 Dalitz decay candidates with MC
1190//
7f3c7cc6 1191 TClonesArray* elec = new TClonesArray("AliKFParticle");
1192 elec->SetOwner(kTRUE);
1193
1194 for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1195 {
1196 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1197 Int_t trackLabel = TMath::Abs(track->GetLabel());
1198
1199 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1200 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1201
1202 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1203 }
1204
1205 return elec;
1206}
1207
1208//--------------------------------------------------------------------------
1209TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1210{
7811d09e 1211//
1212// Find true gammas from pi0 Dalitz decay candidates with MC
1213//
7f3c7cc6 1214 TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1215 gammaPi0->SetOwner(kTRUE);
1216
1217 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1218 {
1219 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1220
1221 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1222 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1223
1224 if( !HaveSameMother(labelv1,labelv2) ) continue;
1225
1226 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1227
1228 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1229
1230 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1231
1232 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1233 }
1234
1235 return gammaPi0;
1236}
1237
1238//--------------------------------------------------------------------------
1239TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1240{
7811d09e 1241//
1242// Find true gammas from gamma candidates with MC
1243//
7f3c7cc6 1244 TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1245 gammaConv->SetOwner(kTRUE);
1246
1247 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1248 {
1249 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1250
1251 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1252 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1253
1254 if( !HaveSameMother(labelv1,labelv2) ) continue;
1255
1256 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1257
1258 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1259
1260 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1261 }
1262
1263 return gammaConv;
1264}
1265
1266//--------------------------------------------------------------
7f3c7cc6 1267void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1268{
7811d09e 1269//
1270// Remove repeated electron candidate tracks
1271// according to the gamma candidate array
1272//
7f3c7cc6 1273 vector<Bool_t> posTag(pos.size(),kTRUE);
1274 vector<Bool_t> negTag(neg.size(),kTRUE);
1275
1276 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1277 {
1278 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1279 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1280
1281 for( UInt_t j=0; j < pos.size(); ++j )
1282 {
1283 if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1284 }
1285
1286 for( UInt_t j=0; j < neg.size(); ++j )
1287 {
1288 if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1289 }
1290 }
1291
1292 CleanArray(pos, posTag);
1293 CleanArray(neg, negTag);
1294}
1295
1296//--------------------------------------------------------------------------
7f3c7cc6 1297void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1298{
7811d09e 1299//
1300// Remove electron candidates from gamma conversions
1301// according to the Psi pair angle
1302//
7f3c7cc6 1303 vector<Bool_t> posTag(pos.size(), kTRUE);
1304 vector<Bool_t> negTag(neg.size(), kTRUE);
1305
1306 for( UInt_t i=0; i < pos.size(); ++i )
1307 {
1308 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1309
1310 for( UInt_t j=0; j < neg.size(); ++j )
1311 {
1312 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1313
1314 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1315 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1316
1317 if(IsFromGammaConversion( psiPair, deltaPhi ))
1318 {
1319 posTag[i] = kFALSE;
1320 negTag[j] = kFALSE;
1321 }
1322 }
1323 }
1324
1325 CleanArray(pos, posTag);
1326 CleanArray(neg, negTag);
1327}
1328
1329//-----------------------------------------------------------------------------------
7f3c7cc6 1330void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1331{
7811d09e 1332//
1333// Remove electron candidates pairs
1334// which have mass not in the range (fMassCutMin,fMassCutMax)
1335//
7f3c7cc6 1336 vector<Bool_t> posTag(pos.size(), kTRUE);
1337 vector<Bool_t> negTag(neg.size(), kTRUE);
1338
1339 for( UInt_t i=0; i < pos.size(); ++i )
1340 {
1341 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1342
1343 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1344 TLorentzVector posLV;
b9295e41 1345 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1346
1347 for( UInt_t j=0; j < neg.size(); ++j )
1348 {
1349 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1350
1351 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1352 TLorentzVector negLV;
b9295e41 1353 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1354
1355 TLorentzVector posnegLV = posLV + negLV;
1356
1357 if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1358 {
1359 posTag[i] = kFALSE;
1360 negTag[j] = kFALSE;
1361 }
1362 }
1363 }
1364
1365 CleanArray(pos, posTag);
1366 CleanArray(neg, negTag);
1367}
1368
1369//-----------------------------------------------------------------------------------------------
1370void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1371{
7811d09e 1372//
1373// Clean the x array according to the tag parameter
1374//
7f3c7cc6 1375 vector<Int_t> tmp;
1376
1377 for(UInt_t i=0; i< x.size(); ++i)
1378 {
1379 if(tag[i]) tmp.push_back(x[i]);
1380 }
1381
1382 x = tmp;
1383}
1384
1385//--------------------------------------------------------------------------
7f3c7cc6 1386void 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)
1387{
7811d09e 1388//
1389// Remove gamma candidates according to
1390// the angle between the plane e+,e- and the gamma
1391//
7f3c7cc6 1392 vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1393
1394 for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1395 {
1396 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1397 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1398 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1399
1400 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1401 {
1402 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1403 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1404 TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1405
1406 // normal vector to x+y- plane
1407 TVector3 planePosNeg = xMom.Cross(yMom);
1408 for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1409 {
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.)
1413 {
1414 gammaTag[i] = kFALSE;
1415 }
1416 }
1417 }
1418 }
1419
1420 // Rebuild gamma candidates array
1421
1422 if(gamma) delete gamma;
1423 gamma = new TClonesArray("AliKFParticle");
1424 gamma->SetOwner(kTRUE);
1425
1426 posGamIdx.clear();
1427 negGamIdx.clear();
1428
1429 for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1430 {
1431 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1432 if(gammaTag[i])
1433 {
1434 new ((*gamma)[j++]) AliKFParticle(*iGamma);
1435 posGamIdx.push_back(fV0Reader->GetPindex(i));
1436 negGamIdx.push_back(fV0Reader->GetNindex(i));
1437 }
1438 }
1439}
1440
1441//--------------------------------------------------------------------------
7f3c7cc6 1442TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1443{
7811d09e 1444//
1445// Find Dalitz pair candidates
1446//
7f3c7cc6 1447 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1448 dalitz->SetOwner(kTRUE);
1449
1450 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1451 {
1452 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1453
1454 TLorentzVector posLV;
b9295e41 1455 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1456
1457 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1458 {
1459 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1460
1461 TLorentzVector negLV;
b9295e41 1462 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
7f3c7cc6 1463
1464 if(fUseAliKF)
1465 {
1466 AliKFParticle posNeg( *posKF, *negKF);
1467
1468 TLorentzVector posNegLV;
1469 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1470 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1471 }
1472 else
1473 {
1474 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1475 }
1476 }
1477 }
1478
1479 return dalitz;
1480}
1481
1482//--------------------------------------------------------------------------
1483TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma)
1484{
7811d09e 1485//
1486// Find pi0 Dalitz decay candidates
1487//
7f3c7cc6 1488 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1489 pi0->SetOwner(kTRUE);
1490
1491 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1492 {
1493 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1494
1495 TLorentzVector posLV;
b9295e41 1496 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1497
1498 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1499 {
1500 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1501
1502 TLorentzVector negLV;
b9295e41 1503 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
7f3c7cc6 1504
1505 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1506 {
1507 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1508
1509 if(fUseAliKF)
1510 {
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);
1515 }
1516 else
1517 {
1518 TLorentzVector gamLV;
1519 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1520
1521 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1522 }
1523 }
1524 }
1525 }
1526
1527 return pi0;
1528}
1529
1530//--------------------------------------------------------------------------
7f3c7cc6 1531TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1532{
7811d09e 1533//
1534// Find true Dalitz pairs from Dalitz pair candidats with MC
1535//
7f3c7cc6 1536 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1537 dalitz->SetOwner(kTRUE);
1538
1539 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1540 {
1541 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1542 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1543 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1544
1545 TLorentzVector posLV;
b9295e41 1546 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1547
1548 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1549
1550 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1551 {
1552 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1553 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1554
1555 if(!IsDalitzPair(posLabel,negLabel)) continue;
1556
1557 if(fUseAliKF)
1558 {
1559 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1560 AliKFParticle posNeg( posKF, negKF);
1561
1562 TLorentzVector posNegLV;
1563 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1564
1565 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1566 }
1567 else // TLorentzVector
1568 {
1569 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1570
1571 TLorentzVector negLV;
b9295e41 1572 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1573
1574 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1575 }
1576 }
1577 }
1578
1579 return dalitz;
1580}
1581
1582//--------------------------------------------------------------------------
7f3c7cc6 1583TClonesArray* 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)
1584{
7811d09e 1585//
1586// Find true pi0 Dalitz decay from pi0 candidates with MC
1587//
7f3c7cc6 1588 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1589 pi0->SetOwner(kTRUE);
1590
1591 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1592 {
1593 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1594 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1595 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1596
1597 TLorentzVector posLV;
b9295e41 1598 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1599
1600 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1601
1602 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1603 {
1604 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1605 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1606
1607 if(!IsDalitzPair(posLabel,negLabel)) continue;
1608
1609 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1610
1611 TLorentzVector negLV;
b9295e41 1612 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1613
1614 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1615
1616 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1617 {
1618 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1619
1620 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1621 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1622
1623 if( !HaveSameMother(labelv1,labelv2) ) continue;
1624
1625 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1626
1627 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1628 if( !HaveSameMother(labelGamma, posLabel) ) continue;
1629
1630 if(fUseAliKF)
1631 {
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);
1636 }
1637 else // TLorentzVector
1638 {
1639 TLorentzVector gamLV;
1640 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1641
1642 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1643 }
1644 }
1645 }
1646 }
1647
1648 return pi0;
1649}
1650
1651//-----------------------------------------------------------------------------------------------
1652void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
1653{
7811d09e 1654//
1655// Update gamma event pool for background computation
1656//
7f3c7cc6 1657 if( fDebug ) AliInfo("=> UpdateGammaPool");
1658
1659 // cycle
1660 for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
1661 {
1662 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
1663 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
1664 ++fGamPoolPos;
1665 if(fGamPoolPos == fPoolMaxSize)
1666 {
1667 fGamPoolPos = 0;
1668 }
1669 }
1670}
1671
1672void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
1673{
7811d09e 1674//
1675// Update electron event pool for background computation
1676//
1677 Int_t multiplicity = fV0Reader->CountESDTracks();
1678
1679
1680 fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
7f3c7cc6 1681}
1682
1683//-----------------------------------------------------------------------------------------------
7f3c7cc6 1684TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
1685{
7811d09e 1686//
1687// Gamma copy from events with same multiplicity and Z
1688//
7f3c7cc6 1689 if( fDebug ) AliInfo("=> GammasFromBGHandler");
7811d09e 1690
1691 Int_t zbin = fBGEventHandler->GetZBinIndex(fV0Reader->GetVertexZ());
1692 Int_t mbin = fBGEventHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1693
1694
7f3c7cc6 1695
1696 TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
1697 gammaPool->SetOwner(kTRUE);
1698
1699 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1700 {
7811d09e 1701 AliGammaConversionKFVector* gammaV0s = fBGEventHandler->GetBGGoodV0s(zbin,mbin,iEventBG);
7f3c7cc6 1702 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
1703 {
1704 new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
1705 }
1706 }
1707
1708 return gammaPool;
1709}
1710
1711//-----------------------------------------------------------------------------------------------
7f3c7cc6 1712TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
1713{
7811d09e 1714//
1715// Electron copy from events with same multiplicity and Z
1716//
7f3c7cc6 1717 if( fDebug ) AliInfo("=> ElectronFromBGHandler");
1718
1719 TClonesArray* electronPool = new TClonesArray("AliKFParticle");
1720 electronPool->SetOwner(kTRUE);
7811d09e 1721
1722 Int_t multiplicity = fV0Reader->CountESDTracks();
1723
1724
1725
7f3c7cc6 1726
1727 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1728 {
7811d09e 1729 AliGammaConversionKFVector* electronNeg = fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
7f3c7cc6 1730 for (UInt_t i = 0; i < electronNeg->size(); ++i )
1731 {
1732 new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
1733 }
1734 }
1735
1736 return electronPool;
1737}
1738
1739//-----------------------------------------------------------------------------------------------
1740Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
1741{
7811d09e 1742//
1743// Get track pid according to MC
1744//
7f3c7cc6 1745 Int_t label = TMath::Abs(t->GetLabel());
1746 Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
1747
1748 switch(pdgCode)
1749 {
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;
1755 }
1756
1757 return -1;
1758}
1759
1760//-----------------------------------------------------------------------------------------------
7811d09e 1761//FIXME PID ITS
1762// NOTE prior should be estimated from data
1763// NOTE: move to config
1764
7f3c7cc6 1765Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
1766{
7811d09e 1767//
1768// Get track pid according to Bayes' formula
1769//
7f3c7cc6 1770 double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
1771 Double_t detectoProb[AliPID::kSPECIES];
1772
1773 if( trackType == kITSsaTrack ) // ITS standalone pid
1774 {
1775 t->GetITSpid( detectoProb );
1776 }
1777 else // global track
1778 {
1779 t->GetESDpid( detectoProb );
1780 }
1781
1782 AliPID bayesPID( detectoProb );
1783 return bayesPID.GetMostProbable( priors );
1784}
1785
1786//-----------------------------------------------------------------------------------------------
1787Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
1788{
7811d09e 1789//
1790// Get track pid according to a n-sigma cut around ITS and/or TPC signals
1791//
7f3c7cc6 1792 if( trackType == kITSsaTrack) // ITS standalone tracks
1793 {
1794 Double_t mom = t->GetP();
1795
1796
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 );
1801
1802 if( nElecSigma < 4. && nElecSigma > -3. && mom < .2 && nPionSigma < -1.5 )
1803 {
1804 return AliPID::kElectron;
1805 }
1806 }
1807 else // global track
1808 {
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 )
1817 {
1818 return AliPID::kElectron;
1819 }
1820 // NOTE: add other particle types
1821 }
1822
1823 return -1;
1824}
1825
1826//-----------------------------------------------------------------------------------------------
7f3c7cc6 1827Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel ) const
1828{
7811d09e 1829//
1830// Returns true if the two particles is a Dalitz pair
1831//
7f3c7cc6 1832 if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
1833
1834 TParticle* pos = fStack->Particle( posLabel );
1835 TParticle* neg = fStack->Particle( negLabel );
1836
1837 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1838 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1839
1840 //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
1841 //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
1842
1843 Int_t motherLabel = pos->GetMother(0);
1844 if( motherLabel < 0 ) return kFALSE;
1845
1846 TParticle* mother = fStack->Particle( motherLabel );
1847
1848 if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1849 if( mother->GetNDaughters() != 3) return kFALSE;
1850 // NOTE: one of them must be a photon
1851
1852 return kTRUE;
1853}
1854
1855//-----------------------------------------------------------------------------------------------
7f3c7cc6 1856Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
1857{
7811d09e 1858//
1859// Returns true if the particle comes from Pi0 -> e+ e- gamma
1860//
7f3c7cc6 1861 Bool_t ePlusFlag = kFALSE;
1862 Bool_t eMinusFlag = kFALSE;
1863 Bool_t gammaFlag = kFALSE;
1864
1865 Int_t motherLabel = fStack->Particle( label )->GetMother(0);
1866
1867 if( motherLabel < 0 ) return kFALSE;
1868
1869 TParticle* mother = fStack->Particle( motherLabel );
1870
1871 if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1872
1873 if ( mother->GetNDaughters() != 3 ) return kFALSE;
1874
1875 for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
1876 {
1877 switch( fStack->Particle(idx)->GetPdgCode())
1878 {
1879 case ::kPositron:
1880 ePlusFlag = kTRUE;
1881 break;
1882 case ::kElectron:
1883 eMinusFlag = kTRUE;
1884 break;
1885 case ::kGamma:
1886 gammaFlag = kTRUE;
1887 break;
1888 }
1889 }
1890
1891 return ( ePlusFlag && eMinusFlag && gammaFlag );
1892}
1893
1894//--------------------------------------------------------------------------
1895Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
1896{
7811d09e 1897//
1898// Returns true if it is a gamma conversion according to psi pair value
1899//
7f3c7cc6 1900 return ( (deltaPhi > fDeltaPhiCutMin && deltaPhi < fDeltaPhiCutMax) &&
1901 TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
1902}
1903
1904//--------------------------------------------------------------------------
1905Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
1906{
7811d09e 1907//
1908// Returns true if it is a gamma conversion according to MC
1909//
7f3c7cc6 1910 if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
1911
1912 TParticle* pos = fStack->Particle( posLabel );
1913 TParticle* neg = fStack->Particle( negLabel );
1914
1915 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1916 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1917
1918 if( pos->GetUniqueID() != kPPair ) return kFALSE;
1919
1920 Int_t motherLabel = pos->GetMother(0);
1921 if( motherLabel < 0 ) return kFALSE;
1922
1923 TParticle* mother = fStack->Particle( motherLabel );
1924
1925 if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
1926
1927 return kTRUE;
1928}
1929
1930//-----------------------------------------------------------------------------------------------
1931Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
1932{
7811d09e 1933//
1934// Returns true if the two particle have the same mother
1935//
7f3c7cc6 1936 if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
1937 return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
1938}
1939
1940//-----------------------------------------------------------------------------------------------
1941Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
1942{
7811d09e 1943//
1944// This angle is a measure for the contribution of the opening in polar
1945// direction Δ0 to the opening angle ξ Pair
1946//
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
1950//
b9295e41 1951 Double_t momPos[3];
1952 Double_t momNeg[3];
1953 if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
1954 if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
7f3c7cc6 1955
b9295e41 1956 TVector3 posDaughter;
1957 TVector3 negDaughter;
7f3c7cc6 1958
b9295e41 1959 posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1960 negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
7f3c7cc6 1961
b9295e41 1962 Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1963 Double_t openingAngle = posDaughter.Angle( negDaughter ); //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
7f3c7cc6 1964 if( openingAngle < 1e-20 ) return 0.;
1965 Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
1966
1967 return psiAngle;
1968}
1969
1970//-----------------------------------------------------------------------------------------------
7811d09e 1971Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const AliKFParticle* xPos, const AliKFParticle* yNeg ) const
7f3c7cc6 1972{
7811d09e 1973//
1974// Get psi pair value
1975//
b9295e41 1976 TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
1977 TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
7f3c7cc6 1978
1979 Double_t deltaTheta = neg.Theta() - pos.Theta();
1980 Double_t openingAngle = pos.Angle( neg );
1981
1982 if( openingAngle < 1e-20 ) return 0.;
1983
1984 return TMath::ASin( deltaTheta/openingAngle );
1985}
1986
1987//-----------------------------------------------------------------------------------------------
1988Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
1989{
7811d09e 1990//
1991// Get psi pair value
1992//
7f3c7cc6 1993 Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
1994 Double_t openingAngle = xPos->Angle( yNeg->Vect() );
1995
1996 if( openingAngle < 1e-20 ) return 0.;
1997
1998 return TMath::ASin( deltaTheta/openingAngle );;
1999}
2000
2001// tmp NOTE: Should go to AliV0Reader
2002//-----------------------------------------------------------------------------------------------
2003/*
2004Double_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
2008 )
2009{
2010 // This angle is a measure for the contribution of the opening in polar
2011 // direction Δ0 to the opening angle ξ Pair
2012
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
2016
2017 const Double_t kForceDeltaPhi = 1.2;
2018 static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
2019
2020 Double_t psiAngle = 0.; // Default value
2021
2022 static Double_t MagnField = fESDEvent->GetMagneticField();
2023 static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
2024
2025 if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
2026
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;
2030
7811d09e 2031 Double_t MomPos[3];
2032 Double_t MomNeg[3];
2033 if( trackPos->GetPxPyPzAt( Rpos, MagnField, MomPos ) == 0 ) trackPos->GetPxPyPz( MomPos );
2034 if( trackNeg->GetPxPyPzAt( Rneg, MagnField, MomNeg ) == 0 ) trackNeg->GetPxPyPz( MomNeg );
7f3c7cc6 2035
7811d09e 2036 TVector3 PosDaughter;
2037 TVector3 NegDaughter;
2038 PosDaughter.SetXYZ( MomPos[0], MomPos[1], MomPos[2] );
2039 NegDaughter.SetXYZ( MomNeg[0], MomNeg[1], MomNeg[2] );
7f3c7cc6 2040
7811d09e 2041 Double_t deltaTheta = NegDaughter.Theta() - PosDaughter.Theta();
2042 Double_t chiPar = PosDaughter.Angle( NegDaughter ); //TMath::ACos( PosDaughter.Dot(NegDaughter) / (NegDaughter.Mag()*PosDaughter.Mag()) );
7f3c7cc6 2043 if( chiPar < 1e-20 ) return psiAngle;
2044
2045 psiAngle = TMath::ASin( deltaTheta / chiPar );
2046
2047 return psiAngle;
2048}
2049*/