]> 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
- remove unused variable nr_contributorsOff
[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),
b9295e41 59// fMCEvent(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{
b9295e41 102 // constructor
103
7f3c7cc6 104 AdoptITSsaTrackCuts();
105 AdoptESDtrackCuts();
106 AdoptESDpidCuts();
107
108 fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
109 fGammaPool->SetOwner(kTRUE);
110}
111
112//-----------------------------------------------------------------------------------------------
113AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz( const char* name ):
114 AliAnalysisTaskSE( name ),
115 fStack(0),
b9295e41 116 // fMCEvent(0),
7f3c7cc6 117 fESDEvent(0),
118 fEposCandidateIndex(),
119 fEnegCandidateIndex(),
120 fGammaCandidatePosIndex(),
121 fGammaCandidateNegIndex(),
122 fGammaCandidates(0),
123 fGammaPool(0),
124 fPoolMaxSize(10),
125 fGamPoolPos(0),
126 fBGEventHandler(0),
127 fOutputContainer(0),
128 fMCTruth(0),
129 fV0Reader(0),
130 fESDpid(0),
131 fESDtrackCuts(0),
132 fITSsaTrackCuts(0),
133 fESDpidCuts(0),
134 fHistograms(0),
135 fStandalone(kFALSE),
136 fDoMC(kFALSE),
137 fComputeBkg(kTRUE),
138 fUseBayesPID(kFALSE),
139 fUseTrackIndexCut(kTRUE),
140 fUsePsiPairCut(kTRUE),
141 fUseMassCut(kFALSE),
142 fUseGammaCut(kFALSE),
143 fReadMagFieldSign(kTRUE),
144 fUseAliKF(kFALSE),
145 fMagFieldSign(1),
b9295e41 146 fkElectronMass(0.00051099891),
7f3c7cc6 147 fPsiPairCut(0.45),
148 fDeltaPhiCutMin(0.),
149 fDeltaPhiCutMax(0.12),
150 fMassCutMin(0.),
151 fMassCutMax(0.1),
152 fNSigmaBelowElecTPCbethe(-2.),
153 fNSigmaAboveElecTPCbethe(3.),
154 fNSigmaAbovePionTPCbethe(3.),
155 fNSigmaAboveKaonTPCbethe(3.),
156 fNSigmaAboveProtonTPCbethe(3.),
157 fTrkSelectionCriteria(kGlobalTrack)
158{
b9295e41 159 //constructor
7f3c7cc6 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{
b9295e41 179 //destructor
7f3c7cc6 180 if( fOutputContainer ) delete fOutputContainer;
181 if( fHistograms ) delete fHistograms;
182 if( fStandalone && fV0Reader ) delete fV0Reader;
183 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
184 if( fESDtrackCuts ) delete fESDtrackCuts;
185 if( fESDpidCuts ) delete fESDpidCuts;
186 if( fGammaCandidates) delete fGammaCandidates;
187 if( fGammaPool ) delete fGammaPool;
188}
189
190//-----------------------------------------------------------------------------------------------
191void AliAnalysisTaskGammaConvDalitz::ConnectInputData(Option_t *option)
192{
193 // Connect Input Data
194 if( fDebug ) AliInfo("=> ConnectInputData");
195
196 AliAnalysisTaskSE::ConnectInputData(option);
197
198 if( fV0Reader == 0 )
199 {
200 AliFatal("There is not pointer to AliV0Reader object!!!");
201 }
202}
203
204//-----------------------------------------------------------------------------------------------
205void AliAnalysisTaskGammaConvDalitz::UserCreateOutputObjects()
206{
b9295e41 207 //usercreateoutputobjects
7f3c7cc6 208 if( fDebug ) AliInfo("=> UserCreateOutputObjects");
209
210 // Create the output container
211 if( fOutputContainer != 0 )
212 {
213 delete fOutputContainer;
214 }
215
216 fOutputContainer = new TList();
217
218 // Add the histograms to the output container
219 fHistograms->GetOutputContainer( fOutputContainer );
220 fOutputContainer->SetOwner(kTRUE);
221
222 PostData( 1, fOutputContainer );
223}
224
225//-----------------------------------------------------------------------------------------------
226// Execute analysis for current event
227void AliAnalysisTaskGammaConvDalitz::UserExec(Option_t */*option*/)
228{
b9295e41 229 //userexec
7f3c7cc6 230 if( fDebug ) AliInfo("=> UserExec");
231
232 if( fV0Reader == 0 )
233 {
234 AliFatal("no pointer to AliV0Reader");
235 return;
236 }
237
238
239
240 // Create list of gamma candidates in standalone mode
241 // otherwise use the created ones by AliAnalysisTaskGammaConversion
242 if( fStandalone )
243 {
244 fV0Reader->UpdateEventByEventData();
245 fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
246 fV0Reader->Initialize();
247 }
248
249 if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
250 {
251 if( fDebug ) AliInfo("no contributors to primary vertex");
252 return;
253 }
254
255 // Get Pointers
256 fBGEventHandler = fV0Reader->GetBGHandler();
257 fESDpid = fV0Reader->GetESDpid();
258 fESDEvent = fV0Reader->GetESDEvent();
259 if(fDoMC)
260 {
261 fStack = fV0Reader->GetMCStack();
262 }
263
264 // Read the magnetic field sign from ESD
265 if ( fReadMagFieldSign == kTRUE )
266 {
267 fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
268 }
269
270 // Process MC information
271 if(fDoMC)
272 {
273 ProcessMCData();
274 }
275
276 if( fStandalone ){
277 while(fV0Reader->NextV0()){}; //SelectGammas
278 fV0Reader->ResetV0IndexNumber();
279 }
280
281 CreateListOfDalitzPairCandidates();
282 ProcessGammaElectronsForDalitzAnalysis();
283
284 PostData( 1, fOutputContainer );
285}
286
287
288void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
289{
290 if( fDebug ) AliInfo("Not to do anything in Terminate");
291}
292
293//-----------------------------------------------------------------------------------------------
294//
295void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
296{
b9295e41 297 //AdoptITSsaTrackCuts
7f3c7cc6 298 if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
299
300 if( esdCuts )
301 {
302 fITSsaTrackCuts = esdCuts;
303 }
304 else
305 {
306 // set defaults cuts
307 fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
308
309 fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
310 fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
311
312 fITSsaTrackCuts->SetMinNClustersITS(2);
313 fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
314 fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
315
316 fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
317 fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
318
319 fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
320 }
321}
322
323//-----------------------------------------------------------------------------------------------
324//
325void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
326{
b9295e41 327 //AdoptESDtrackCuts
7f3c7cc6 328 if( fESDtrackCuts ) delete fESDtrackCuts;
329
330 if( esdCuts )
331 {
332 fESDtrackCuts = esdCuts;
333 }
334 else
335 {
336 // set defaults cuts
337 fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
338
339 fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
340 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
341
342 fESDtrackCuts->SetMinNClustersITS(2);
343 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
344 fESDtrackCuts->SetRequireITSRefit(kTRUE);
345
346 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
347 fESDtrackCuts->SetMaxNsigmaToVertex(3);
348
349 fESDtrackCuts->SetMinNClustersTPC(80);
350 fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
351 fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
352 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
353 }
354}
355
356//-----------------------------------------------------------------------------------------------
357void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
358{
b9295e41 359 //AdoptESDpidCuts
7f3c7cc6 360 if( fESDpidCuts ) delete fESDpidCuts;
361 if( esdPIDCuts )
362 {
363 fESDpidCuts = esdPIDCuts;
364 }
365 else // set defaults cuts
366 {
367 fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
368 // fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
369 fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
370 }
371}
372
373//-----------------------------------------------------------------------------------------------
374// Process generation
375void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
376{
b9295e41 377 //ProcessMCData
7f3c7cc6 378 if( fDebug ) AliInfo("=> ProcessMCData");
379
380 fHistograms->FillTable("Table_Generation", 0); //number of events
381
382 for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
383 {
384 TParticle* iParticle = fStack->Particle( i );
385 if( !iParticle ) continue;
386
387 if ( i >= fStack->GetNprimary() ) continue; // is primary?
388 if ( iParticle->GetPdgCode() != ::kPi0 ) continue; // is Pi0?
389
390 if( iParticle->GetNDaughters() == 2 &&
391 fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
392 fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
393 {
394 fHistograms->FillTable("Table_Generation", 1); // pi0 -> gg
395 }
396
397 if ( iParticle->GetNDaughters() != 3 ) continue; // Num == 3 (e+,e-,gamma)
398
399 // Check for Pi0 Dalitz decay
400 TParticle* eposPi0 = 0;
401 TParticle* enegPi0 = 0;
402 TParticle* gammaPi0 = 0;
403
404 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
405 {
406 switch(fStack->Particle(idxPi0)->GetPdgCode())
407 {
408 case ::kPositron:
409 eposPi0 = fStack->Particle(idxPi0);
410 break;
411 case ::kElectron:
412 enegPi0 = fStack->Particle(idxPi0);
413 break;
414 case ::kGamma:
415 gammaPi0 = fStack->Particle(idxPi0);
416 break;
417 }
418 }
419
420 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
421
422 // found a Pi0 Dalitz decay
423
424 fHistograms->FillTable("Table_Generation", 2);
425 fHistograms->FillHistogram("MC_Pi0Dalitz_P", iParticle->P());
426 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt", iParticle->Pt());
427 fHistograms->FillHistogram("MC_Pi0Dalitz_Eta", iParticle->Eta());
428 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
429 fHistograms->FillHistogram("MC_EposDalitz_Pt", eposPi0->Pt());
430 fHistograms->FillHistogram("MC_EposDalitz_Eta", eposPi0->Eta());
431 fHistograms->FillHistogram("MC_EnegDalitz_Pt", enegPi0->Pt());
432 fHistograms->FillHistogram("MC_EnegDalitz_Eta", enegPi0->Eta());
433 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Pt", gammaPi0->Pt());
434 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Eta", gammaPi0->Eta());
435
436 // Angle between the gamma and the plane e+e-
437 TVector3 ePosMom( eposPi0->Px(), eposPi0->Py(), eposPi0->Pz() );
438 TVector3 eNegMom( enegPi0->Px(), enegPi0->Py(), enegPi0->Pz() );
439 TVector3 gamMom( gammaPi0->Px(), gammaPi0->Py() , gammaPi0->Pz() );
440 TVector3 planeEposEneg = eNegMom.Cross( ePosMom );
441 Double_t anglePlaneGamma = planeEposEneg.Angle(gamMom);
442
443 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
444
445 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle", anglePlaneGamma);
446 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_P", anglePlaneGamma, gammaPi0->P());
447 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_Pt", anglePlaneGamma, gammaPi0->Pt());
448
449
450 // check for gamma conversion
451 Bool_t daugGammaElectron = kFALSE;
452 Bool_t daugGammaPositron = kFALSE; // acceptance
453 Bool_t daugGammaElectronAll = kFALSE;
454 Bool_t daugGammaPositronAll = kFALSE;
455
456 // is the gamma converted? -> has 2 daughter e+e-
457 // are e+ e- from gamma in the acceptance for the V0s
458 if( gammaPi0->GetNDaughters() >= 2 )
459 {
460 for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
461 {
462 TParticle* tmpDaughter = fStack->Particle(tIndex);
463
464 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
465 if( tmpDaughter->GetPdgCode() == ::kElectron )
466 { // e+
467 daugGammaElectronAll = kTRUE;
468
469 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
470 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
471 (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
472 {
473 daugGammaElectron = kTRUE;
474 }
475 }
476 else if( tmpDaughter->GetPdgCode() == ::kPositron )
477 {
478 daugGammaPositronAll = kTRUE;
479 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
480 ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
481 (tmpDaughter->R() < fV0Reader->GetMaxRCut() ) )
482 {
483 daugGammaPositron = kTRUE;
484 }
485 }
486 }
487 }
488
489
490 if( daugGammaElectronAll && daugGammaPositronAll )
491 {
492 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
493 fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
494 }
495
496 Float_t etaMin, etaMax;
497 fESDtrackCuts->GetEtaRange( etaMin, etaMax );
498
499 // e+e- pair within acceptance
500 if ( TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
501 {
502 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt", eposPi0->Pt());
503 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta", eposPi0->Eta());
504 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt", enegPi0->Pt());
505 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta", enegPi0->Eta());
506 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_EposPt_vs_EnegPt", eposPi0->Pt(), enegPi0->Pt());
507 }
508
509 // Pi0 (e+e-gamma) within acceptance
510 //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
511
512 if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
513 TMath::Abs( eposPi0->Eta() ) < etaMax && TMath::Abs( enegPi0->Eta() ) < etaMax )
514 {
515 fHistograms->FillTable("Table_Generation",3); //
516 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt",iParticle->Pt());
517 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Eta",iParticle->Eta());
518 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Radius",iParticle->R());
519 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Pt",gammaPi0->Pt());
520 fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Eta",gammaPi0->Eta());
521 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt",eposPi0->Pt());
522 fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta",eposPi0->Eta());
523 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt",enegPi0->Pt());
524 fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta",enegPi0->Eta());
525 fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
526 fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
527
528 // Pi0 within acceptance with gamma converted
529
530 if ( daugGammaElectron && daugGammaPositron )
531 {
532 fHistograms->FillTable("Table_Generation",4); //
533
534 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
535 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
536 fHistograms->FillHistogram("MC_Acceptance_GC_EposDalitz_Pt",eposPi0->Pt());
537 fHistograms->FillHistogram("MC_Acceptance_GC_EposDalitz_Eta",eposPi0->Eta());
538 fHistograms->FillHistogram("MC_Acceptance_GC_EnegDalitz_Pt",enegPi0->Pt());
539 fHistograms->FillHistogram("MC_Acceptance_GC_EnegDalitz_Eta",enegPi0->Eta());
540 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Pt",gammaPi0->Pt());
541 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Eta",gammaPi0->Eta());
542 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle",anglePlaneGamma);
543 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle_vs_Pt",anglePlaneGamma,gammaPi0->Pt());
544 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
545 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
546 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
547 }
548 }
549 }
550}
551
552//-----------------------------------------------------------------------------------------------
553// Dalitz pair candidates
554void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
555{
b9295e41 556 //CreateListOfDalitzPairCandidates
7f3c7cc6 557 if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
558
559 fEposCandidateIndex.clear();
560 fEnegCandidateIndex.clear();
561
562 fHistograms->FillTable("Table_Cuts", 0);
563
564 for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
565 {
566 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
567 if ( !iTrack ) continue;
568
569
570 Double_t p[3];
571
572 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
573
574 TVector3 iMom(p[0],p[1],p[2]);
575
576 //
577 // Check track cuts and find track type
578 //
579
580 Bool_t isTrackAccepted = 0;
581 Int_t trackType = -1;
582 switch(fTrkSelectionCriteria)
583 {
584 case kITSsaTrack:
585 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
586 trackType = kITSsaTrack;
587 break;
588
589 case kGlobalTrack:
590 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
591 trackType = kGlobalTrack;
592 break;
593
594 case kITSsaGlobalTrack:
595 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
596 {
597 isTrackAccepted = kTRUE;
598 if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
599 else trackType = kGlobalTrack;
600 }
601 break;
602 }
603
604 if(!isTrackAccepted) continue;
605
606 //
607 // PID
608 //
609
610 Int_t pid=-1;
611 Int_t pidMC=-1;
612
613 if(fUseBayesPID)
614 {
615 pid = GetBayesPid(iTrack,trackType);
616 }
617 else
618 {
619 pid = GetNSigmaPid(iTrack,trackType);
620 }
621
622 if( fDoMC )
623 {
624 pidMC = GetMonteCarloPid(iTrack);
625 // pid table
626 Int_t iLabel = TMath::Abs(iTrack->GetLabel());
627 TParticle* iParticle = fStack->Particle(iLabel);
628 FillPidTable(iParticle, pid);
629 }
630
631 // ITS standalone tracks
632 if( trackType == kITSsaTrack)
633 {
634 Double_t mom = iTrack->GetP();
635 Double_t signal = iTrack->GetITSsignal();
636
637 fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
638
639 if( pid == AliPID::kElectron )
640 {
641
642 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
643 if(fDoMC && pid == pidMC)
644 {
645 fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
646 }
647 }
648
649 if( fDoMC && pidMC == AliPID::kElectron)
650 {
651 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
652 }
653 }
654
655 else // global tracks
656 {
657 const AliExternalTrackParam *in = iTrack->GetInnerParam();
658 Double_t mom = in->GetP();
659 Double_t signal = iTrack->GetTPCsignal();
660
661 fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
662
663 if( fDoMC && pidMC == AliPID::kElectron )
664 {
665 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
666 }
667
668 if( pid == AliPID::kElectron )
669 {
670 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
671 if(fDoMC && pid == pidMC)
672 {
673 fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
674 }
675 }
676 }
677
678 if( AliPID::kElectron != pid) continue;
679
680 // electron track candidates from here
681
682 if( iTrack->GetSign() > 0 )
683 {
684 fEposCandidateIndex.push_back(i);
685 }
686 else
687 {
688 fEnegCandidateIndex.push_back(i);
689 }
690 }
691
692 // gamma candidates
693 GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
694
695 if(fDoMC)
696 {
697 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
698 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
699 delete pi0Dalitz;
700 }
701
702 if(fUseTrackIndexCut) // remove repeated tracks
703 {
704 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
705
706 if(fDoMC)
707 {
708 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
709 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
710 delete pi0Dalitz;
711 }
712 }
713
714 if(fUsePsiPairCut) // remove electrons from gamma conversions
715 {
716 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
717
718 if(fDoMC)
719 {
720 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
721 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
722 delete pi0Dalitz;
723 }
724 }
725
726 if( fUseMassCut )
727 {
728 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
729
730 if(fDoMC)
731 {
732 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
733 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
734 delete pi0Dalitz;
735 }
736 }
737
738 if(fUseGammaCut)
739 {
740 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
741
742 if(fDoMC)
743 {
744 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
745 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
746 delete pi0Dalitz;
747 }
748 }
749}
750
751//-----------------------------------------------------------------------------------------------
752void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
753{
b9295e41 754 //ProcessGammaElectronsForDalitzAnalysis
7f3c7cc6 755 if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
756
757 fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
758
759 TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
760
761 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
762 {
763 AliKFParticle* epos = (AliKFParticle*) ePosCandidates->At(i);
764 fHistograms->FillHistogram("ESD_EposCandidates_Pt", epos->GetPt());
765 fHistograms->FillHistogram("ESD_EposCandidates_Eta", epos->GetEta());
766 fHistograms->FillTable( "Table_Reconstruction", 1);
767 }
768
769 TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
770
771 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
772 {
773 AliKFParticle* eneg = (AliKFParticle*) eNegCandidates->At(i);
774 fHistograms->FillHistogram("ESD_EnegCandidates_Pt", eneg->GetPt());
775 fHistograms->FillHistogram("ESD_EnegCandidates_Eta", eneg->GetEta());
776 fHistograms->FillTable( "Table_Reconstruction", 2);
777 }
778
779 TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
780 for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
781 {
782 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
783
784 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
785 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Mass", dalitz->M());
786 }
787
788 // gamma candidates
789 for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
790 {
791 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
792 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
793 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
794 }
795
796 // psi pair for all candidates
797 //if(fUsePsiPairCut)
798 FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
799
800 // Angle epos,eneg gamma
801 FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
802 FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
803
804 TClonesArray* pi0Candidates = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaCandidates);
805 for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
806 {
807 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
808
809 fHistograms->FillHistogram("ESD_Pi0_P", pi0->P());
810 fHistograms->FillHistogram("ESD_Pi0_Pt", pi0->Pt());
811 fHistograms->FillHistogram("ESD_Pi0_Eta", pi0->Eta());
812 fHistograms->FillHistogram("ESD_Pi0_Y", pi0->Rapidity());
813 fHistograms->FillHistogram("ESD_Pi0_Phi", pi0->Phi());
814 fHistograms->FillHistogram("ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
815 fHistograms->FillHistogram("ESD_Pi0_Mass", pi0->M());
816 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
817 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Y", pi0->Rapidity(),pi0->M());
818 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
819 }
820
821 delete dalitzPairCandidates;
822 delete pi0Candidates;
823
824 if(fComputeBkg)
825 {
826 // 1) e+e- with with gammas used in the signal
827 TClonesArray* pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaPool);
828
829 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
830 {
831 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
832 fHistograms->FillHistogram("ESD_BKG_PrevGamma", pi0->M());
833 fHistograms->FillHistogram("ESD_BKG_PrevGamma_vs_Pt", pi0->Pt(), pi0->M());
834 }
835
836 if(ePosCandidates->GetEntriesFast() > 0 &&
837 eNegCandidates->GetEntriesFast() > 0 &&
838 fGammaCandidates->GetEntriesFast() > 0)
839 {
840 UpdateGammaPool(fGammaCandidates);
841 }
842
843 delete pi0Bkg;
844
845 // 2) e+e- with gammas from a pool of events
3c048fe4 846/* TClonesArray* gammaBGHandler = GammasFromBGHandler();
7f3c7cc6 847 pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, gammaBGHandler);
848
849 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
850 {
851 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
852 fHistograms->FillHistogram("ESD_BKG_BGHandler", pi0->M());
853 fHistograms->FillHistogram("ESD_BKG_BGHandler_vs_Pt", pi0->Pt(), pi0->M());
854 }
855
856 delete pi0Bkg;
857
858 // 3) e+ with e-, gamma from a pool of events
859 TClonesArray* elecBGHandler = ElectronFromBGHandler();
860 pi0Bkg = FindPi0Dalitz(ePosCandidates, elecBGHandler, gammaBGHandler);
861
862 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
863 {
864 TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
865 fHistograms->FillHistogram("ESD_BKG_Electron", pi0->M());
866 fHistograms->FillHistogram("ESD_BKG_Electron_vs_Pt", pi0->Pt(), pi0->M());
867 }
868
869 if(eNegCandidates->GetEntriesFast() > 0)
870 {
871 UpdateElectronPool(eNegCandidates);
872 }
873
874 delete gammaBGHandler;
875 delete elecBGHandler;
876 delete pi0Bkg;
3c048fe4 877*/
7f3c7cc6 878 }
879
880 delete ePosCandidates;
881 delete eNegCandidates;
882
883 if(fDoMC)
884 {
885 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
886 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
887 {
888 AliKFParticle* epos = (AliKFParticle*) ePosPi0Dalitz->At(i);
889 fHistograms->FillHistogram("MC_ESD_EposDalitz_Pt", epos->GetPt());
890 fHistograms->FillHistogram("MC_ESD_EposDalitz_Eta", epos->GetEta());
891 fHistograms->FillTable( "Table_Reconstruction", 3);
892 }
893
894 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
895 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
896 {
897 AliKFParticle* eneg = (AliKFParticle*) eNegPi0Dalitz->At(i);
898 fHistograms->FillHistogram("MC_ESD_EnegDalitz_Pt", eneg->GetPt());
899 fHistograms->FillHistogram("MC_ESD_EnegDalitz_Eta", eneg->GetEta());
900 fHistograms->FillTable( "Table_Reconstruction", 4);
901 }
902
903 TClonesArray* dalitzPair = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex);
904 for(Int_t i=0; i < dalitzPair->GetEntriesFast(); ++i)
905 {
906 TLorentzVector* dalitz = (TLorentzVector*) dalitzPair->At(i);
907 fHistograms->FillHistogram("MC_ESD_DalitzPair_Pt", dalitz->Pt());
908 fHistograms->FillHistogram("MC_ESD_DalitzPair_Mass", dalitz->M());
909 fHistograms->FillHistogram( "Table_Reconstruction", 5 );
910 }
911
912 // psi pair for dalitz pairs
913 //if(fUsePsiPairCut)
914 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_DalitzPair_PsiPair_vs_DPhi");
915
916 delete ePosPi0Dalitz;
917 delete eNegPi0Dalitz;
918 delete dalitzPair;
919
920 // all gammas
921 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
922 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
923 {
924 AliKFParticle* iGamma = (AliKFParticle*) gamma->At(i);
925 fHistograms->FillHistogram("MC_ESD_Gamma_Pt", iGamma->GetPt());
926 fHistograms->FillHistogram("MC_ESD_Gamma_Eta", iGamma->GetEta());
927 }
928
929 delete gamma;
930
931 // gamma from pi0 dalitz
932 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
933 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
934 {
935 AliKFParticle* iGamma = (AliKFParticle*) gammaPi0Dalitz->At(i);
936 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Pt", iGamma->GetPt());
937 fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Eta", iGamma->GetEta());
938 fHistograms->FillTable( "Table_Reconstruction", 6);
939 }
940
941 delete gammaPi0Dalitz;
942
943 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
944
945 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
946 {
947 TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
948
949 fHistograms->FillHistogram("MC_ESD_Pi0_P", pi0->P());
950 fHistograms->FillHistogram("MC_ESD_Pi0_Pt", pi0->Pt());
951 fHistograms->FillHistogram("MC_ESD_Pi0_Eta", pi0->Eta());
952 fHistograms->FillHistogram("MC_ESD_Pi0_Y", pi0->Rapidity());
953 fHistograms->FillHistogram("MC_ESD_Pi0_Phi", pi0->Phi());
954 fHistograms->FillHistogram("MC_ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
955 fHistograms->FillHistogram("MC_ESD_Pi0_Mass", pi0->M());
956 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
957 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Y", pi0->Rapidity(),pi0->M());
958 fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
959 fHistograms->FillHistogram( "Table_Reconstruction", 7);
960 }
961 delete pi0Dalitz;
962
963 // psi pair for electrons from gamma conversions assuming they came from main vertex
964 // if(fUsePsiPairCut)
965 for(UInt_t i=0; i < fEposCandidateIndex.size(); ++i)
966 {
967 AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
968 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
969
970 for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
971 {
972 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
973 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
974
975 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
976
977 Double_t psiPair = GetPsiPair(posTrack, negTrack);
978 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
979
980 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
981 }
982 }
983 // FIXME: eta -> e+e-gamma
984 }
985}
986
987//--------------------------------------------------------------------------
988Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
989{
b9295e41 990 //Rapidity
991 const double epsilon=1.e-16;
7f3c7cc6 992
b9295e41 993 if(p->Energy() - TMath::Abs(p->Pz()) < epsilon )
7f3c7cc6 994 {
995 return 1.e10;
996 }
997 return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
998}
999
1000//--------------------------------------------------------------------------
1001void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1002{
b9295e41 1003 //FillPsiPair
7f3c7cc6 1004 for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1005 {
1006 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1007 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1008 {
1009 AliKFParticle* negKF = (AliKFParticle*) neg->At(j);
1010 Double_t psiPair = GetPsiPair(posKF, negKF);
1011 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negKF->GetPhi() - posKF->GetPhi());
1012 fHistograms->FillHistogram(hName, deltaPhi, psiPair);
1013 }
1014 }
1015}
1016
1017//--------------------------------------------------------------------------
1018void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1019{
b9295e41 1020 //FillAngle
7f3c7cc6 1021 for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1022 {
1023 AliKFParticle* xKF = (AliKFParticle*) x->At(i);
1024 TVector3 xMom(xKF->Px(),xKF->Py(),xKF->Pz());
1025 for( Int_t j=0; j < y->GetEntriesFast(); ++j )
1026 {
1027 AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1028 TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1029 fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1030 }
1031 }
1032}
1033
1034//--------------------------------------------------------------------------
1035
1036void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1037{
b9295e41 1038 //FillPidTable
7f3c7cc6 1039 Int_t iGen=-1;
1040 switch(TMath::Abs(p->GetPdgCode()))
1041 {
1042 case ::kElectron: iGen=0; break;
1043 case ::kMuonMinus: iGen=1; break;
1044 case ::kPiPlus: iGen=2; break;
1045 case ::kKPlus: iGen=3; break;
1046 case ::kProton: iGen=4; break;
1047 default: iGen=-1;
1048 }
1049
1050 int jRec=-1;
1051 if(pid > -1 && pid < 5) jRec = pid;
1052
1053 if ((iGen > -1) && (jRec > -1))
1054 {
1055 fHistograms->FillTable("Table_PID", iGen, jRec);
1056 // sum
1057 fHistograms->FillTable("Table_PID", iGen, 5);
1058 fHistograms->FillTable("Table_PID", 5, jRec);
1059 }
1060}
1061
1062//--------------------------------------------------------------------------
1063void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1064{
b9295e41 1065 //GetGammaCandidates
7f3c7cc6 1066 posIndex.clear();
1067 negIndex.clear();
1068
1069 if(gamma) delete gamma;
1070
1071 TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1072
1073 gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1074 gamma->SetOwner(kTRUE);
1075
1076 // make a copy
1077 for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1078 {
1079 AliKFParticle* gamKF = (AliKFParticle*)gammaV0->At(i);
1080 new ((*gamma)[i]) AliKFParticle(*gamKF);
1081 posIndex.push_back(fV0Reader->GetPindex(i));
1082 negIndex.push_back(fV0Reader->GetNindex(i));
1083 }
1084}
1085
1086//--------------------------------------------------------------------------
1087TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1088{
b9295e41 1089 //IndexToAliKFParticle
7f3c7cc6 1090 TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1091 indexKF->SetOwner(kTRUE);
1092
1093 for(UInt_t i = 0; i < index.size(); ++i)
1094 {
1095 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1096 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1097 }
1098
1099 return indexKF;
1100}
1101
1102//--------------------------------------------------------------------------
1103TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1104{
b9295e41 1105 //FindElectronFromPi0Dalitz
7f3c7cc6 1106 TClonesArray* elec = new TClonesArray("AliKFParticle");
1107 elec->SetOwner(kTRUE);
1108
1109 for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1110 {
1111 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1112 Int_t trackLabel = TMath::Abs(track->GetLabel());
1113
1114 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1115 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1116
1117 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1118 }
1119
1120 return elec;
1121}
1122
1123//--------------------------------------------------------------------------
1124TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1125{
b9295e41 1126 //FindGammaFromPi0Dalitz
7f3c7cc6 1127 TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1128 gammaPi0->SetOwner(kTRUE);
1129
1130 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1131 {
1132 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1133
1134 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1135 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1136
1137 if( !HaveSameMother(labelv1,labelv2) ) continue;
1138
1139 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1140
1141 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1142
1143 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1144
1145 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1146 }
1147
1148 return gammaPi0;
1149}
1150
1151//--------------------------------------------------------------------------
1152TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1153{
b9295e41 1154 //FindGamma
7f3c7cc6 1155 TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1156 gammaConv->SetOwner(kTRUE);
1157
1158 for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1159 {
1160 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1161
1162 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1163 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1164
1165 if( !HaveSameMother(labelv1,labelv2) ) continue;
1166
1167 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1168
1169 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1170
1171 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1172 }
1173
1174 return gammaConv;
1175}
1176
1177//--------------------------------------------------------------
1178// Remove repeated electron candidate tracks according to the gamma candidate array
1179void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1180{
b9295e41 1181 //ESDtrackIndexCut
7f3c7cc6 1182 vector<Bool_t> posTag(pos.size(),kTRUE);
1183 vector<Bool_t> negTag(neg.size(),kTRUE);
1184
1185 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1186 {
1187 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1188 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1189
1190 for( UInt_t j=0; j < pos.size(); ++j )
1191 {
1192 if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1193 }
1194
1195 for( UInt_t j=0; j < neg.size(); ++j )
1196 {
1197 if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1198 }
1199 }
1200
1201 CleanArray(pos, posTag);
1202 CleanArray(neg, negTag);
1203}
1204
1205//--------------------------------------------------------------------------
1206// Remove electron candidates from gamma conversions according to the Psi pair angle
1207void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1208{
b9295e41 1209 //PsiPairCut
7f3c7cc6 1210 vector<Bool_t> posTag(pos.size(), kTRUE);
1211 vector<Bool_t> negTag(neg.size(), kTRUE);
1212
1213 for( UInt_t i=0; i < pos.size(); ++i )
1214 {
1215 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1216
1217 for( UInt_t j=0; j < neg.size(); ++j )
1218 {
1219 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1220
1221 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1222 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1223
1224 if(IsFromGammaConversion( psiPair, deltaPhi ))
1225 {
1226 posTag[i] = kFALSE;
1227 negTag[j] = kFALSE;
1228 }
1229 }
1230 }
1231
1232 CleanArray(pos, posTag);
1233 CleanArray(neg, negTag);
1234}
1235
1236//-----------------------------------------------------------------------------------
1237// Remove electron candidates pairs with mass not in the range (fMassCutMin,fMassCutMax)
1238void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1239{
b9295e41 1240 //MassCut
7f3c7cc6 1241 vector<Bool_t> posTag(pos.size(), kTRUE);
1242 vector<Bool_t> negTag(neg.size(), kTRUE);
1243
1244 for( UInt_t i=0; i < pos.size(); ++i )
1245 {
1246 AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1247
1248 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1249 TLorentzVector posLV;
b9295e41 1250 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1251
1252 for( UInt_t j=0; j < neg.size(); ++j )
1253 {
1254 AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1255
1256 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1257 TLorentzVector negLV;
b9295e41 1258 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1259
1260 TLorentzVector posnegLV = posLV + negLV;
1261
1262 if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1263 {
1264 posTag[i] = kFALSE;
1265 negTag[j] = kFALSE;
1266 }
1267 }
1268 }
1269
1270 CleanArray(pos, posTag);
1271 CleanArray(neg, negTag);
1272}
1273
1274//-----------------------------------------------------------------------------------------------
1275void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1276{
b9295e41 1277 //CleanArray
7f3c7cc6 1278 vector<Int_t> tmp;
1279
1280 for(UInt_t i=0; i< x.size(); ++i)
1281 {
1282 if(tag[i]) tmp.push_back(x[i]);
1283 }
1284
1285 x = tmp;
1286}
1287
1288//--------------------------------------------------------------------------
1289// Remove gamma candidates according to the angle between the plane e+,e- and the gamma
1290void 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)
1291{
b9295e41 1292 //AngleEposEnegGammaCut
7f3c7cc6 1293 vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1294
1295 for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1296 {
1297 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1298 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1299 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1300
1301 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1302 {
1303 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1304 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1305 TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1306
1307 // normal vector to x+y- plane
1308 TVector3 planePosNeg = xMom.Cross(yMom);
1309 for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1310 {
1311 AliKFParticle* gamKF = (AliKFParticle*)candidates->At(i);
1312 TVector3 gamMom(gamKF->Px(),gamKF->Py(),gamKF->Pz());
1313 if (planePosNeg.Angle(gamMom) < 1. || planePosNeg.Angle(gamMom) > 2.)
1314 {
1315 gammaTag[i] = kFALSE;
1316 }
1317 }
1318 }
1319 }
1320
1321 // Rebuild gamma candidates array
1322
1323 if(gamma) delete gamma;
1324 gamma = new TClonesArray("AliKFParticle");
1325 gamma->SetOwner(kTRUE);
1326
1327 posGamIdx.clear();
1328 negGamIdx.clear();
1329
1330 for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1331 {
1332 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1333 if(gammaTag[i])
1334 {
1335 new ((*gamma)[j++]) AliKFParticle(*iGamma);
1336 posGamIdx.push_back(fV0Reader->GetPindex(i));
1337 negGamIdx.push_back(fV0Reader->GetNindex(i));
1338 }
1339 }
1340}
1341
1342//--------------------------------------------------------------------------
1343//
1344TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1345{
b9295e41 1346 //FindDalitzPair
7f3c7cc6 1347 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1348 dalitz->SetOwner(kTRUE);
1349
1350 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1351 {
1352 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1353
1354 TLorentzVector posLV;
b9295e41 1355 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1356
1357 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1358 {
1359 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1360
1361 TLorentzVector negLV;
b9295e41 1362 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
7f3c7cc6 1363
1364 if(fUseAliKF)
1365 {
1366 AliKFParticle posNeg( *posKF, *negKF);
1367
1368 TLorentzVector posNegLV;
1369 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1370 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1371 }
1372 else
1373 {
1374 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1375 }
1376 }
1377 }
1378
1379 return dalitz;
1380}
1381
1382//--------------------------------------------------------------------------
1383TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma)
1384{
b9295e41 1385 //FindPi0Dalitz
7f3c7cc6 1386 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1387 pi0->SetOwner(kTRUE);
1388
1389 for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1390 {
1391 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1392
1393 TLorentzVector posLV;
b9295e41 1394 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
7f3c7cc6 1395
1396 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1397 {
1398 AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1399
1400 TLorentzVector negLV;
b9295e41 1401 negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
7f3c7cc6 1402
1403 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1404 {
1405 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1406
1407 if(fUseAliKF)
1408 {
1409 AliKFParticle posNegGam( *posKF, *negKF, *gamKF );
1410 TLorentzVector posNegGamLV;
1411 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1412 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1413 }
1414 else
1415 {
1416 TLorentzVector gamLV;
1417 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1418
1419 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1420 }
1421 }
1422 }
1423 }
1424
1425 return pi0;
1426}
1427
1428//--------------------------------------------------------------------------
1429// Only with montecarlo
1430TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1431{
b9295e41 1432 //FindDalitzPair
7f3c7cc6 1433 TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1434 dalitz->SetOwner(kTRUE);
1435
1436 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1437 {
1438 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1439 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1440 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1441
1442 TLorentzVector posLV;
b9295e41 1443 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1444
1445 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1446
1447 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1448 {
1449 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1450 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1451
1452 if(!IsDalitzPair(posLabel,negLabel)) continue;
1453
1454 if(fUseAliKF)
1455 {
1456 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1457 AliKFParticle posNeg( posKF, negKF);
1458
1459 TLorentzVector posNegLV;
1460 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1461
1462 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1463 }
1464 else // TLorentzVector
1465 {
1466 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1467
1468 TLorentzVector negLV;
b9295e41 1469 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1470
1471 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1472 }
1473 }
1474 }
1475
1476 return dalitz;
1477}
1478
1479//--------------------------------------------------------------------------
1480// Only with montecarlo
1481TClonesArray* 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)
1482{
b9295e41 1483 //FindPi0Dalitz
7f3c7cc6 1484 TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1485 pi0->SetOwner(kTRUE);
1486
1487 for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1488 {
1489 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1490 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1491 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1492
1493 TLorentzVector posLV;
b9295e41 1494 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
7f3c7cc6 1495
1496 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1497
1498 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1499 {
1500 AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1501 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1502
1503 if(!IsDalitzPair(posLabel,negLabel)) continue;
1504
1505 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1506
1507 TLorentzVector negLV;
b9295e41 1508 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
7f3c7cc6 1509
1510 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1511
1512 for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1513 {
1514 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1515
1516 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1517 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1518
1519 if( !HaveSameMother(labelv1,labelv2) ) continue;
1520
1521 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1522
1523 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1524 if( !HaveSameMother(labelGamma, posLabel) ) continue;
1525
1526 if(fUseAliKF)
1527 {
1528 AliKFParticle posNegGam( posKF, negKF, *gamKF );
1529 TLorentzVector posNegGamLV;
1530 posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1531 new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1532 }
1533 else // TLorentzVector
1534 {
1535 TLorentzVector gamLV;
1536 gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1537
1538 new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1539 }
1540 }
1541 }
1542 }
1543
1544 return pi0;
1545}
1546
1547//-----------------------------------------------------------------------------------------------
1548void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
1549{
b9295e41 1550 //UpdateGammaPool
7f3c7cc6 1551 if( fDebug ) AliInfo("=> UpdateGammaPool");
1552
1553 // cycle
1554 for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
1555 {
1556 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
1557 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
1558 ++fGamPoolPos;
1559 if(fGamPoolPos == fPoolMaxSize)
1560 {
1561 fGamPoolPos = 0;
1562 }
1563 }
1564}
1565
1566void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
1567{
b9295e41 1568 //UpdateElectronPool
7f3c7cc6 1569 fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetMultiplicity());
1570}
1571
1572//-----------------------------------------------------------------------------------------------
1573// returns a copy of gammas from events with same multiplicity and Z
1574TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
1575{
b9295e41 1576 //GammasFromBGHandler
7f3c7cc6 1577 if( fDebug ) AliInfo("=> GammasFromBGHandler");
1578
1579 TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
1580 gammaPool->SetOwner(kTRUE);
1581
1582 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1583 {
1584 AliGammaConversionKFVector* gammaV0s = fV0Reader->GetBGGoodV0s( iEventBG );
1585 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
1586 {
1587 new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
1588 }
1589 }
1590
1591 return gammaPool;
1592}
1593
1594//-----------------------------------------------------------------------------------------------
1595// returns a copy of electron from events with same multiplicity and Z
1596TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
1597{
b9295e41 1598 //ElectronFromBGHandler
7f3c7cc6 1599 if( fDebug ) AliInfo("=> ElectronFromBGHandler");
1600
1601 TClonesArray* electronPool = new TClonesArray("AliKFParticle");
1602 electronPool->SetOwner(kTRUE);
1603
1604 for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1605 {
1606 AliGammaConversionKFVector* electronNeg = fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetMultiplicity());
1607 for (UInt_t i = 0; i < electronNeg->size(); ++i )
1608 {
1609 new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
1610 }
1611 }
1612
1613 return electronPool;
1614}
1615
1616//-----------------------------------------------------------------------------------------------
1617Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
1618{
b9295e41 1619 //GetMonteCarloPid
7f3c7cc6 1620 Int_t label = TMath::Abs(t->GetLabel());
1621 Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
1622
1623 switch(pdgCode)
1624 {
1625 case ::kElectron: return AliPID::kElectron;
1626 case ::kMuonMinus: return AliPID::kMuon;
1627 case ::kPiPlus: return AliPID::kPion;
1628 case ::kKPlus: return AliPID::kKaon;
1629 case ::kProton: return AliPID::kProton;
1630 }
1631
1632 return -1;
1633}
1634
1635//-----------------------------------------------------------------------------------------------
1636Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
1637{
b9295e41 1638 //GetBayesPid
7f3c7cc6 1639 //FIXME PID ITS
1640 // NOTE prior should be estimated from data
1641 // NOTE: move to config
1642 double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
1643 Double_t detectoProb[AliPID::kSPECIES];
1644
1645 if( trackType == kITSsaTrack ) // ITS standalone pid
1646 {
1647 t->GetITSpid( detectoProb );
1648 }
1649 else // global track
1650 {
1651 t->GetESDpid( detectoProb );
1652 }
1653
1654 AliPID bayesPID( detectoProb );
1655 return bayesPID.GetMostProbable( priors );
1656}
1657
1658//-----------------------------------------------------------------------------------------------
1659Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
1660{
b9295e41 1661 //GetNSigmaPid
7f3c7cc6 1662 // ITS Number of sigmas (FIXME: add new fESDpidCuts)
1663 // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
1664
1665 if( trackType == kITSsaTrack) // ITS standalone tracks
1666 {
1667 Double_t mom = t->GetP();
1668
1669
1670 // ITS Number of sigmas (FIXME: add new fESDpidCuts)
1671 // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
1672 Double_t nElecSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kElectron );
1673 Double_t nPionSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kPion );
1674
1675 if( nElecSigma < 4. && nElecSigma > -3. && mom < .2 && nPionSigma < -1.5 )
1676 {
1677 return AliPID::kElectron;
1678 }
1679 }
1680 else // global track
1681 {
1682 Double_t nElecSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kElectron );
1683 Double_t nPionSigma = fESDpid->NumberOfSigmasTPC(t, AliPID::kPion );
1684 Double_t nKaonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kKaon ));
1685 Double_t nProtonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kProton));
1686 if( nElecSigma > fNSigmaBelowElecTPCbethe && nElecSigma < fNSigmaAboveElecTPCbethe &&
1687 nPionSigma > fNSigmaAbovePionTPCbethe && //NOTE mom > 0.5
1688 nKaonSigma > fNSigmaAboveKaonTPCbethe &&
1689 nProtonSigma > fNSigmaAboveProtonTPCbethe )
1690 {
1691 return AliPID::kElectron;
1692 }
1693 // NOTE: add other particle types
1694 }
1695
1696 return -1;
1697}
1698
1699//-----------------------------------------------------------------------------------------------
1700// returns true if the two particles is a Dalitz pair
1701Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel ) const
1702{
b9295e41 1703 //IsDalitzPair
7f3c7cc6 1704 if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
1705
1706 TParticle* pos = fStack->Particle( posLabel );
1707 TParticle* neg = fStack->Particle( negLabel );
1708
1709 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1710 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1711
1712 //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
1713 //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
1714
1715 Int_t motherLabel = pos->GetMother(0);
1716 if( motherLabel < 0 ) return kFALSE;
1717
1718 TParticle* mother = fStack->Particle( motherLabel );
1719
1720 if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1721 if( mother->GetNDaughters() != 3) return kFALSE;
1722 // NOTE: one of them must be a photon
1723
1724 return kTRUE;
1725}
1726
1727//-----------------------------------------------------------------------------------------------
1728// returns true if the particle comes from Pi0 -> e+ e- gamma
1729Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
1730{
b9295e41 1731 //IsPi0DalitzDaughter
7f3c7cc6 1732 Bool_t ePlusFlag = kFALSE;
1733 Bool_t eMinusFlag = kFALSE;
1734 Bool_t gammaFlag = kFALSE;
1735
1736 Int_t motherLabel = fStack->Particle( label )->GetMother(0);
1737
1738 if( motherLabel < 0 ) return kFALSE;
1739
1740 TParticle* mother = fStack->Particle( motherLabel );
1741
1742 if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1743
1744 if ( mother->GetNDaughters() != 3 ) return kFALSE;
1745
1746 for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
1747 {
1748 switch( fStack->Particle(idx)->GetPdgCode())
1749 {
1750 case ::kPositron:
1751 ePlusFlag = kTRUE;
1752 break;
1753 case ::kElectron:
1754 eMinusFlag = kTRUE;
1755 break;
1756 case ::kGamma:
1757 gammaFlag = kTRUE;
1758 break;
1759 }
1760 }
1761
1762 return ( ePlusFlag && eMinusFlag && gammaFlag );
1763}
1764
1765//--------------------------------------------------------------------------
1766Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
1767{
1768 // triangle
1769 return ( (deltaPhi > fDeltaPhiCutMin && deltaPhi < fDeltaPhiCutMax) &&
1770 TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
1771}
1772
1773//--------------------------------------------------------------------------
1774Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
1775{
b9295e41 1776 //IsFromGammaConversion
7f3c7cc6 1777 if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
1778
1779 TParticle* pos = fStack->Particle( posLabel );
1780 TParticle* neg = fStack->Particle( negLabel );
1781
1782 if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1783 if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1784
1785 if( pos->GetUniqueID() != kPPair ) return kFALSE;
1786
1787 Int_t motherLabel = pos->GetMother(0);
1788 if( motherLabel < 0 ) return kFALSE;
1789
1790 TParticle* mother = fStack->Particle( motherLabel );
1791
1792 if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
1793
1794 return kTRUE;
1795}
1796
1797//-----------------------------------------------------------------------------------------------
1798Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
1799{
b9295e41 1800 //HaveSameMother
7f3c7cc6 1801 if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
1802 return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
1803}
1804
1805//-----------------------------------------------------------------------------------------------
1806Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
1807{
1808 // This angle is a measure for the contribution of the opening in polar
1809 // direction Δ0 to the opening angle ξ Pair
1810
1811 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
1812 // Master Thesis. Thorsten Dahms. 2005
1813 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
1814
b9295e41 1815 Double_t momPos[3];
1816 Double_t momNeg[3];
1817 if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
1818 if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
7f3c7cc6 1819
b9295e41 1820 TVector3 posDaughter;
1821 TVector3 negDaughter;
7f3c7cc6 1822
b9295e41 1823 posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1824 negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
7f3c7cc6 1825
b9295e41 1826 Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1827 Double_t openingAngle = posDaughter.Angle( negDaughter ); //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
7f3c7cc6 1828 if( openingAngle < 1e-20 ) return 0.;
1829 Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
1830
1831 return psiAngle;
1832}
1833
1834//-----------------------------------------------------------------------------------------------
b9295e41 1835Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliKFParticle* xPos, const AliKFParticle* yNeg ) const // FIXME: alikf const
7f3c7cc6 1836{
b9295e41 1837 //GetPsiPair
1838 TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
1839 TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
7f3c7cc6 1840
1841 Double_t deltaTheta = neg.Theta() - pos.Theta();
1842 Double_t openingAngle = pos.Angle( neg );
1843
1844 if( openingAngle < 1e-20 ) return 0.;
1845
1846 return TMath::ASin( deltaTheta/openingAngle );
1847}
1848
1849//-----------------------------------------------------------------------------------------------
1850Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
1851{
b9295e41 1852 //GetPsiPair
7f3c7cc6 1853 Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
1854 Double_t openingAngle = xPos->Angle( yNeg->Vect() );
1855
1856 if( openingAngle < 1e-20 ) return 0.;
1857
1858 return TMath::ASin( deltaTheta/openingAngle );;
1859}
1860
1861// tmp NOTE: Should go to AliV0Reader
1862//-----------------------------------------------------------------------------------------------
1863/*
1864Double_t AliAnalysisTaskGammaConvDalitz::GetPsiAngleV0(
1865 Double_t radiusVO, // radius at XY of VO vertex
1866 const AliExternalTrackParam* trackPos, // pos. track parm. at V0 vertex
1867 const AliExternalTrackParam* trackNeg // neg. track parm. at V0 vertex
1868 )
1869{
1870 // This angle is a measure for the contribution of the opening in polar
1871 // direction Δ0 to the opening angle ξ Pair
1872
1873 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
1874 // Master Thesis. Thorsten Dahms. 2005
1875 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
1876
1877 const Double_t kForceDeltaPhi = 1.2;
1878 static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
1879
1880 Double_t psiAngle = 0.; // Default value
1881
1882 static Double_t MagnField = fESDEvent->GetMagneticField();
1883 static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
1884
1885 if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
1886
1887 // compute propagation radius for a fixed angle
1888 Double_t Rpos = radiusVO + trackPos->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
1889 Double_t Rneg = radiusVO + trackNeg->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
1890
b9295e41 1891 Double_t momPos[3];
1892 Double_t momNeg[3];
1893 if( trackPos->GetPxPyPzAt( Rpos, MagnField, momPos ) == 0 ) trackPos->GetPxPyPz( momPos );
1894 if( trackNeg->GetPxPyPzAt( Rneg, MagnField, momNeg ) == 0 ) trackNeg->GetPxPyPz( momNeg );
7f3c7cc6 1895
b9295e41 1896 TVector3 posDaughter;
1897 TVector3 negDaughter;
1898 posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1899 negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
7f3c7cc6 1900
b9295e41 1901 Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1902 Double_t chiPar = posDaughter.Angle( negDaughter ); //TMath::ACos( posDaughter.Dot(negDaughter) / (negDaughter.Mag()*posDaughter.Mag()) );
7f3c7cc6 1903 if( chiPar < 1e-20 ) return psiAngle;
1904
1905 psiAngle = TMath::ASin( deltaTheta / chiPar );
1906
1907 return psiAngle;
1908}
1909*/