]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Macro to read and display SSD common mode values (Panos)
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConversion.cxx
CommitLineData
d7d7e825 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
5 * Version 1.1 *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////
17//---------------------------------------------
18// Class used to do analysis on conversion pairs
19//---------------------------------------------
20////////////////////////////////////////////////
21
22// root
23#include <TChain.h>
24
25// analysis
26#include "AliAnalysisTaskGammaConversion.h"
27#include "AliStack.h"
28#include "AliLog.h"
29#include "AliESDtrackCuts.h"
30#include "TNtuple.h"
31
32class AliKFVertex;
33class AliAODHandler;
34class AliAODEvent;
35class ALiESDEvent;
36class AliMCEvent;
37class AliMCEventHandler;
38class AliESDInputHandler;
39class AliAnalysisManager;
40class Riostream;
41class TFile;
42class TInterpreter;
43class TSystem;
44class TROOT;
45
46ClassImp(AliAnalysisTaskGammaConversion)
47
48
49AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
50AliAnalysisTaskSE(),
51 fV0Reader(NULL),
52 fStack(NULL),
53 fESDEvent(NULL),
54 fOutputContainer(NULL),
55 fHistograms(NULL),
56 fDoMCTruth(kFALSE),
57 fDoNeutralMeson(kFALSE),
58 fDoJet(kFALSE),
59 fDoChic(kFALSE),
6c84d371 60 fKFReconstructedGammasTClone(NULL),
61 fCurrentEventPosElectronTClone(NULL),
62 fCurrentEventNegElectronTClone(NULL),
63 fKFReconstructedGammasCutTClone(NULL),
64 fPreviousEventTLVNegElectronTClone(NULL),
65 fPreviousEventTLVPosElectronTClone(NULL),
66// fKFReconstructedGammas(),
d7d7e825 67 fElectronv1(),
68 fElectronv2(),
6c84d371 69// fCurrentEventPosElectron(),
70// fCurrentEventNegElectron(),
71// fKFReconstructedGammasCut(),
72// fPreviousEventTLVNegElectron(),
73// fPreviousEventTLVPosElectron(),
d7d7e825 74 fElectronMass(-1),
75 fGammaMass(-1),
76 fPi0Mass(-1),
77 fEtaMass(-1),
78 fGammaWidth(-1),
79 fPi0Width(-1),
80 fEtaWidth(-1),
81 fMinOpeningAngleGhostCut(0.),
82 fEsdTrackCuts(0),
83 fCalculateBackground(kFALSE),
84 fWriteNtuple(kFALSE),
85 fGammaNtuple(NULL),
86 fNeutralMesonNtuple(NULL),
87 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 88 fChargedParticles(NULL),
d7d7e825 89 fChargedParticlesId(),
90 fGammaPtHighest(0.),
91 fMinPtForGammaJet(1.),
92 fMinIsoConeSize(0.2),
93 fMinPtIsoCone(0.7),
94 fMinPtGamChargedCorr(0.5),
95 fMinPtJetCone(0.5),
96 fLeadingChargedIndex(-1),
97 fAODBranch(NULL),
6c84d371 98 fAODBranchName("GammaConv")//,
99// fAODObjects(NULL)
d7d7e825 100{
101 // Default constructor
102 // Common I/O in slot 0
103 DefineInput (0, TChain::Class());
104 DefineOutput(0, TTree::Class());
105
106 // Your private output
107 DefineOutput(1, TList::Class());
108
109 // Define standard ESD track cuts for Gamma-hadron correlation
110 SetESDtrackCuts();
111}
112
113AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
114 AliAnalysisTaskSE(name),
115 fV0Reader(NULL),
116 fStack(NULL),
117 fESDEvent(NULL),
118 fOutputContainer(0x0),
119 fHistograms(NULL),
120 fDoMCTruth(kFALSE),
121 fDoNeutralMeson(kFALSE),
122 fDoJet(kFALSE),
123 fDoChic(kFALSE),
6c84d371 124 fKFReconstructedGammasTClone(NULL),
125 fCurrentEventPosElectronTClone(NULL),
126 fCurrentEventNegElectronTClone(NULL),
127 fKFReconstructedGammasCutTClone(NULL),
128 fPreviousEventTLVNegElectronTClone(NULL),
129 fPreviousEventTLVPosElectronTClone(NULL),
130 // fKFReconstructedGammas(),
d7d7e825 131 fElectronv1(),
132 fElectronv2(),
6c84d371 133 // fCurrentEventPosElectron(),
134 // fCurrentEventNegElectron(),
135 // fKFReconstructedGammasCut(),
136 // fPreviousEventTLVNegElectron(),
137 // fPreviousEventTLVPosElectron(),
d7d7e825 138 fElectronMass(-1),
139 fGammaMass(-1),
140 fPi0Mass(-1),
141 fEtaMass(-1),
142 fGammaWidth(-1),
143 fPi0Width(-1),
144 fEtaWidth(-1),
145 fMinOpeningAngleGhostCut(0.),
146 fEsdTrackCuts(0),
147 fCalculateBackground(kFALSE),
148 fWriteNtuple(kFALSE),
149 fGammaNtuple(NULL),
150 fNeutralMesonNtuple(NULL),
151 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 152 fChargedParticles(NULL),
d7d7e825 153 fChargedParticlesId(),
154 fGammaPtHighest(0.),
155 fMinPtForGammaJet(1.),
156 fMinIsoConeSize(0.2),
157 fMinPtIsoCone(0.7),
158 fMinPtGamChargedCorr(0.5),
159 fMinPtJetCone(0.5),
160 fLeadingChargedIndex(-1),
161 fAODBranch(NULL),
6c84d371 162 fAODBranchName("GammaConv")//,
163 // fAODObjects(NULL)
d7d7e825 164{
165 // Common I/O in slot 0
166 DefineInput (0, TChain::Class());
167 DefineOutput(0, TTree::Class());
168
169 // Your private output
170 DefineOutput(1, TList::Class());
171
172 // Define standard ESD track cuts for Gamma-hadron correlation
173 SetESDtrackCuts();
174}
175
176AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
177{
178 // Remove all pointers
179
180 if(fOutputContainer){
181 fOutputContainer->Clear() ;
182 delete fOutputContainer ;
183 }
184 if(fHistograms){
185 delete fHistograms;
186 }
187 if(fV0Reader){
188 delete fV0Reader;
189 }
190 if (fAODBranch) {
191 fAODBranch->Clear();
192 delete fAODBranch ;
193 }
194}
195
196
197void AliAnalysisTaskGammaConversion::Init()
198{
199 // Initialization
200 AliLog::SetGlobalLogLevel(AliLog::kError);
201}
202void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
203{
204 // SetESDtrackCuts
205
206 fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
207//standard cuts from:
208//http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
209//fEsdTrackCuts->SetMinNClustersTPC(50);
210//fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
211//fEsdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
212 fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
213 fEsdTrackCuts->SetRequireITSRefit(kTRUE);
214 fEsdTrackCuts->SetMaxNsigmaToVertex(3);
215 fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
216 // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
217
218}
219
220void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)
221{
222 // Execute analysis for current event
223
224 ConnectInputData("");
225
226 //Each event needs an empty branch
227 fAODBranch->Clear();
228
6c84d371 229 if(fKFReconstructedGammasTClone == NULL){
230 fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
231 }
232 if(fCurrentEventPosElectronTClone == NULL){
233 fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
234 }
235 if(fCurrentEventNegElectronTClone == NULL){
236 fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
237 }
238 if(fKFReconstructedGammasCutTClone == NULL){
239 fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
240 }
241 if(fPreviousEventTLVNegElectronTClone == NULL){
242 fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
243 }
244 if(fPreviousEventTLVPosElectronTClone == NULL){
245 fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
246 }
247 if(fChargedParticles == NULL){
248 fChargedParticles = new TClonesArray("AliESDtrack",0);
249 }
250
251 //clear TClones
252 fKFReconstructedGammasTClone->Clear();
253 fCurrentEventPosElectronTClone->Clear();
254 fCurrentEventNegElectronTClone->Clear();
255 fKFReconstructedGammasCutTClone->Clear();
256 fPreviousEventTLVNegElectronTClone->Clear();
257 fPreviousEventTLVPosElectronTClone->Clear();
258
d7d7e825 259 //clear vectors
6c84d371 260 // fKFReconstructedGammas.clear();
d7d7e825 261 fElectronv1.clear();
262 fElectronv2.clear();
6c84d371 263 // fCurrentEventPosElectron.clear();
264 // fCurrentEventNegElectron.clear();
265 // fKFReconstructedGammasCut.clear();
d7d7e825 266
6c84d371 267 fChargedParticles->Clear();
d7d7e825 268 fChargedParticlesId.clear();
269
270 //Clear the data in the v0Reader
271 fV0Reader->UpdateEventByEventData();
272
273
274 // Process the MC information
275 if(fDoMCTruth){
276 ProcessMCData();
277 }
278
279 //Process the v0 information with no cuts
280 ProcessV0sNoCut();
281
282 // Process the v0 information
283 ProcessV0s();
284
285 //Fill Gamma AOD
286 FillAODWithConversionGammas() ;
287
288 //calculate background if flag is set
289 if(fCalculateBackground){
290 CalculateBackground();
291 }
292
293
294
295
296 // Process reconstructed gammas
297 if(fDoNeutralMeson == kTRUE){
298 ProcessGammasForNeutralMesonAnalysis();
299 }
300
301 CheckV0Efficiency();
302
303 //Process reconstructed gammas electrons for Chi_c Analysis
6c84d371 304 if(fDoChic == kTRUE){
d7d7e825 305 ProcessGammaElectronsForChicAnalysis();
306 }
307 // Process reconstructed gammas for gamma Jet/hadron correlations
308 if(fDoJet == kTRUE){
309 ProcessGammasForGammaJetAnalysis();
310 }
311
312 PostData(1, fOutputContainer);
313
314}
315
316void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){
317 // see header file for documentation
318
319 if(fV0Reader == NULL){
320 // Write warning here cuts and so on are default if this ever happens
321 }
322 fV0Reader->Initialize();
323}
324
325
326
327void AliAnalysisTaskGammaConversion::ProcessMCData(){
328 // see header file for documentation
329
330 fStack = fV0Reader->GetMCStack();
331
332 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
333 return; // aborts if the primary vertex does not have contributors.
334 }
335
336 for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
337 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
338
339 if (!particle) {
340 //print warning here
341 continue;
342 }
343
344 ///////////////////////Begin Chic Analysis/////////////////////////////
345
346
347 if(particle->GetPdgCode() == 443){//Is JPsi
348
349 if(particle->GetNDaughters()==2){
350 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
351 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
352 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
353 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
354 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
355 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
356
357 if( TMath::Abs(daug0->Eta()) < 0.9){
358 if(daug0->GetPdgCode() == -11)
359 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
360 else
361 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
362
363 }
364 if(TMath::Abs(daug1->Eta()) < 0.9){
365 if(daug1->GetPdgCode() == -11)
366 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
367 else
368 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
369 }
370 }
371 }
372 }
373 // const int CHI_C0 = 10441;
374 // const int CHI_C1 = 20443;
375 // const int CHI_C2 = 445
376 if(particle->GetPdgCode() == 22){//gamma from JPsi
377 if(particle->GetMother(0) > -1){
378 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
379 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
380 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
381 if(TMath::Abs(particle->Eta()) < 1.2)
382 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
383 }
384 }
385 }
386 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
387 if( particle->GetNDaughters() == 2){
388 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
389 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
390
391 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
392 if( daug0->GetPdgCode() == 443){
393 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
394 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
395 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
396 fHistograms->FillTable("Table_Electrons",18);
397
398 }//if
399 else if (daug1->GetPdgCode() == 443){
400 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
401 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
402 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
403 fHistograms->FillTable("Table_Electrons",18);
404 }//else if
405 }//gamma o Jpsi
406 }//GetNDaughters
407 }
408
409
410 /////////////////////End Chic Analysis////////////////////////////
411
412
413 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
414
415 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
416
417 Double_t tmpPhi=particle->Phi();
418
419 if(particle->Phi()> TMath::Pi()){
420 tmpPhi = particle->Phi()-(2*TMath::Pi());
421 }
422
423 Double_t rapidity;
424 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
425 rapidity=0;
426 }
427 else{
428 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
429 }
430
431 //process the gammas
432 if (particle->GetPdgCode() == 22){
433
434 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
435 continue; // no photon as mothers!
436 }
437
438 if(particle->GetMother(0) >= fStack->GetNprimary()){
439 continue; // the gamma has a mother, and it is not a primary particle
440 }
441
d7d7e825 442 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
443 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
444 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
445 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
446 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
447
448
449 if(particle->GetMother(0) < 0){ // direct gamma
450 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
451 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
452 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
453 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
454 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
455 }
456
457
458 // looking for conversion (electron + positron from pairbuilding (= 5) )
459 TParticle* ePos = NULL;
460 TParticle* eNeg = NULL;
461
462 if(particle->GetNDaughters() >= 2){
463 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
464 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
465 if(tmpDaughter->GetUniqueID() == 5){
466 if(tmpDaughter->GetPdgCode() == 11){
467 eNeg = tmpDaughter;
468 }
469 else if(tmpDaughter->GetPdgCode() == -11){
470 ePos = tmpDaughter;
471 }
472 }
473 }
474 }
475
476
477 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
478 continue;
479 }
480
481
482 Double_t ePosPhi = ePos->Phi();
483 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
484
485 Double_t eNegPhi = eNeg->Phi();
486 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
487
488
489 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
490 continue; // no reconstruction below the Pt cut
491 }
492
493 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
494 continue;
495 }
496
497 if(ePos->R()>fV0Reader->GetMaxRCut()){
498 continue; // cuts on distance from collision point
499 }
500
501 if(TMath::Abs(ePos->Vz()) > 240){
502 continue; // outside material
503 }
504
505 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
506 continue; // line cut to exclude regions where we do not reconstruct
507 }
508
509 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
510 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
511 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
512 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
513 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
514 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
515
516 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
517 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
518 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
519 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
520
521 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
522 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
523 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
524 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
525
526
527
528 //cout << "filled histos for converted gamma, ePos, eNeg" << endl;
529
530 // begin Mapping
531 Int_t rBin = fHistograms->GetRBin(ePos->R());
532 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
533
534 TString nameMCMappingPhiR="";
535 nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);
536 fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
537
538 TString nameMCMappingPhi="";
539 nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);
540 fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
541
542 TString nameMCMappingR="";
543 nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);
544 fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
545
546 TString nameMCMappingPhiInR="";
547 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);
548 fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
549 //end mapping
550
551 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
552 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
553 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
554 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
555
556 //cout << "mapping is done" << endl;
557
558
559 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
560 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
561 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
562 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
563 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
564 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
565
566 } // end direct gamma
567 else{ // mother exits
6c84d371 568 /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
d7d7e825 569 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
570 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
571 ){
572 fMCGammaChic.push_back(particle);
573 }
6c84d371 574 */
d7d7e825 575 } // end if mother exits
576 } // end if particle is a photon
577
578
579
580 // process motherparticles (2 gammas as daughters)
581 // the motherparticle had already to pass the R and the eta cut, but no line cut.
582 // the line cut is just valid for the conversions!
583
584 if(particle->GetNDaughters() == 2){
585
586 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
587 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
588
589 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
590
591 // Check the acceptance for both gammas
01b7fdcc 592 Bool_t gammaEtaCut = kTRUE;
593 if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
594 Bool_t gammaRCut = kTRUE;
595 if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
d7d7e825 596
597
598
599 // check for conversions now -> have to pass eta, R and line cut!
600 Bool_t daughter0Electron = kFALSE;
601 Bool_t daughter0Positron = kFALSE;
602 Bool_t daughter1Electron = kFALSE;
603 Bool_t daughter1Positron = kFALSE;
604
605 if(daughter0->GetNDaughters() >= 2){ // first gamma
606 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
607 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
608 if(tmpDaughter->GetUniqueID() == 5){
609 if(tmpDaughter->GetPdgCode() == 11){
610 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
611 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
612 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
613 daughter0Electron = kTRUE;
614 }
615 }
616
617 }
618 }
619 else if(tmpDaughter->GetPdgCode() == -11){
620 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
621 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
622 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
623 daughter0Positron = kTRUE;
624 }
625 }
626
627 }
628
629 }
630 }
631 }
632 }
633
634
635 if(daughter1->GetNDaughters() >= 2){ // second gamma
636 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
637 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
638 if(tmpDaughter->GetUniqueID() == 5){
639 if(tmpDaughter->GetPdgCode() == 11){
640 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
641 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
642 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
643 daughter1Electron = kTRUE;
644 }
645 }
646
647 }
648 }
649 else if(tmpDaughter->GetPdgCode() == -11){
650 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
651 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
652 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
653 daughter1Positron = kTRUE;
654 }
655 }
656
657 }
658
659 }
660 }
661 }
662 }
663
664
665
666
667 if(particle->GetPdgCode()==111){ //Pi0
668 if( iTracks >= fStack->GetNprimary()){
669 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
670 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
671 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
672 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
673 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
674 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
675 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
676 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
677 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
678
01b7fdcc 679 if(gammaEtaCut && gammaRCut){
d7d7e825 680 //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
681 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
682 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
683 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
684 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
685 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
686 }
687 }
688 }
689 else{
690 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
691 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
692 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
693 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
694 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
695 fHistograms->FillHistogram("MC_Pi0_R", particle->R());
696 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
697 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
698 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
699
01b7fdcc 700 if(gammaEtaCut && gammaRCut){
d7d7e825 701 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
702 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
703 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
704 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
705 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
706 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
707 fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
708 }
709 }
710 }
711 }
712
713 if(particle->GetPdgCode()==221){ //Eta
714 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
715 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
716 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
717 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
718 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
719 fHistograms->FillHistogram("MC_Eta_R", particle->R());
720 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
721 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
722 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
723
01b7fdcc 724 if(gammaEtaCut && gammaRCut){
d7d7e825 725 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
726 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
727 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
728 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
729 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
730 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
731 fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
732 }
733
734 }
735
736 }
737
738 // all motherparticles with 2 gammas as daughters
739 fHistograms->FillHistogram("MC_Mother_R", particle->R());
740 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
741 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
742 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
743 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
744 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
745 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
746 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
747 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
748 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
749 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
750
01b7fdcc 751 if(gammaEtaCut && gammaRCut){
d7d7e825 752 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
753 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
754 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
755 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
756 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
757 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
758 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
759 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
760
761 }
762
763
764 } // end passed R and eta cut
765
766 } // end if(particle->GetNDaughters() == 2)
767
768 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
769
770 //cout << "right before the end of processMCdata" << endl;
771
772} // end ProcessMCData
773
774
775
776void AliAnalysisTaskGammaConversion::FillNtuple(){
777 //Fills the ntuple with the different values
778
779 if(fGammaNtuple == NULL){
780 return;
781 }
782 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
783 for(Int_t i=0;i<numberOfV0s;i++){
784 Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
785 AliESDv0 * cV0 = fV0Reader->GetV0(i);
786 Double_t negPID=0;
787 Double_t posPID=0;
788 fV0Reader->GetPIDProbability(negPID,posPID);
789 values[0]=cV0->GetOnFlyStatus();
790 values[1]=fV0Reader->CheckForPrimaryVertex();
791 values[2]=negPID;
792 values[3]=posPID;
793 values[4]=fV0Reader->GetX();
794 values[5]=fV0Reader->GetY();
795 values[6]=fV0Reader->GetZ();
796 values[7]=fV0Reader->GetXYRadius();
797 values[8]=fV0Reader->GetMotherCandidateNDF();
798 values[9]=fV0Reader->GetMotherCandidateChi2();
799 values[10]=fV0Reader->GetMotherCandidateEnergy();
800 values[11]=fV0Reader->GetMotherCandidateEta();
801 values[12]=fV0Reader->GetMotherCandidatePt();
802 values[13]=fV0Reader->GetMotherCandidateMass();
803 values[14]=fV0Reader->GetMotherCandidateWidth();
804 // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
805 values[16]=fV0Reader->GetOpeningAngle();
806 values[17]=fV0Reader->GetNegativeTrackEnergy();
807 values[18]=fV0Reader->GetNegativeTrackPt();
808 values[19]=fV0Reader->GetNegativeTrackEta();
809 values[20]=fV0Reader->GetNegativeTrackPhi();
810 values[21]=fV0Reader->GetPositiveTrackEnergy();
811 values[22]=fV0Reader->GetPositiveTrackPt();
812 values[23]=fV0Reader->GetPositiveTrackEta();
813 values[24]=fV0Reader->GetPositiveTrackPhi();
814 values[25]=fV0Reader->HasSameMCMother();
815 if(values[25] != 0){
816 values[26]=fV0Reader->GetMotherMCParticlePDGCode();
817 values[15]=fV0Reader->GetMotherMCParticle()->Pt();
818 }
819 fTotalNumberOfAddedNtupleEntries++;
820 fGammaNtuple->Fill(values);
821 }
822 fV0Reader->ResetV0IndexNumber();
823
824}
825
826void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
827 // Process all the V0's without applying any cuts to it
828
829 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
830 for(Int_t i=0;i<numberOfV0s;i++){
831 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
832
833 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
834 return;
835 }
836
837 if(fDoMCTruth){
838
839 if(fV0Reader->HasSameMCMother() == kFALSE){
840 continue;
841 }
842
843 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
844 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
845
846 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
847 continue;
848 }
849 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
850 continue;
851 }
852
853 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
854
855 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
856 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
857 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
858 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
859 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
860 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
861 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
862 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
863 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
864 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
865
866 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
867 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
868
869 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
870 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
871 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
872 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
873
874 //store MCTruth properties
875 fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
876 fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
877 fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
878
879
880 }
881 }
882 }
883 fV0Reader->ResetV0IndexNumber();
884}
885
886void AliAnalysisTaskGammaConversion::ProcessV0s(){
887 // see header file for documentation
888
889 if(fWriteNtuple == kTRUE){
890 FillNtuple();
891 }
892
893 Int_t nSurvivingV0s=0;
894 while(fV0Reader->NextV0()){
895 nSurvivingV0s++;
896
897
898 //-------------------------- filling v0 information -------------------------------------
899 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
900 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
901 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
902 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
903
904 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
905 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
906 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
907 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
908
909 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
910 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
911 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
912 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
913
914 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
915 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
916 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
917 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
918 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
919 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
920 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
921 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
922 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
923 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
924
925 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
926 fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
927
928
929
930 // begin mapping
931 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
932 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
933 Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
934
935 TString nameESDMappingPhiR="";
936 nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);
937 fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
938
939 TString nameESDMappingPhi="";
940 nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);
941 fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
942
943 TString nameESDMappingR="";
944 nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);
945 fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
946
947 TString nameESDMappingPhiInR="";
948 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);
949 fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
950 // end mapping
951
6c84d371 952 new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
953
954 // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
d7d7e825 955 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
956 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
957
958
959 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
960 if(fDoMCTruth){
961
d7d7e825 962 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
6c84d371 963 // TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle(); // not used any longer
d7d7e825 964
d7d7e825 965 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
d7d7e825 966
967 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
968 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
969 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
970 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
971 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
972 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
973 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
974 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
975 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
976 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
977 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
978 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
979 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
980 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
981
982 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
983 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
984
985
986 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
987 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
988 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
989 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
990
991 //store MCTruth properties
992 fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
993 fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
994 fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
995
996 //resolution
997 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
998 Double_t esdpt = fV0Reader->GetMotherCandidatePt();
999 Double_t resdPt = 0;
1000 if(mcpt != 0){
1001 resdPt = ((esdpt - mcpt)/mcpt)*100;
1002 }
1003
1004 fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);
1005 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1006 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1007
1008 Double_t resdZ = 0;
1009 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1010 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
1011 }
1012
1013 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1014 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1015 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1016
1017 Double_t resdR = 0;
1018 if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1019 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
1020 }
1021
1022 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1023 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1024 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1025 fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);
1026 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
d7d7e825 1027 }//if(fDoMCTruth)
1028 }//while(fV0Reader->NextV0)
1029 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1030 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1031
1032 //cout << "nearly at the end of doMCTruth" << endl;
1033
1034}
1035
1036void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1037 // Fill AOD with reconstructed Gamma
1038
6c84d371 1039 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1040 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
d7d7e825 1041 //Create AOD particle object from AliKFParticle
1042
1043 /* AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1044 //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1045 //but this means that I have to work a little bit more in my side.
1046 //AODPWG4Particle objects are simpler and lighter, I think
1047 AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
1048 gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1049 gamma.SetCaloLabel(-1,-1); //How to get the MC label of the 2 electrons that form the gamma?
1050 gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
1051 gamma.SetPdg(AliCaloPID::kPhotonConv); //photon id
1052 gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
1053
1054 //Add it to the aod list
1055 Int_t i = fAODBranch->GetEntriesFast();
1056 new((*fAODBranch)[i]) AliAODPWG4Particle(gamma);
1057 */
6c84d371 1058 // AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1059 AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
d7d7e825 1060 AliGammaConversionAODObject aodObject;
1061 aodObject.SetPx(gammakf->GetPx());
1062 aodObject.SetPy(gammakf->GetPy());
1063 aodObject.SetPz(gammakf->GetPz());
1064 aodObject.SetLabel1(fElectronv1[gammaIndex]);
1065 aodObject.SetLabel2(fElectronv2[gammaIndex]);
1066 Int_t i = fAODBranch->GetEntriesFast();
6c84d371 1067 new((*fAODBranch)[i]) AliGammaConversionAODObject(aodObject);
d7d7e825 1068 }
1069
1070}
1071
1072
1073void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1074 // see header file for documentation
1075
6c84d371 1076 // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1077 // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
1078 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1079 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
d7d7e825 1080
6c84d371 1081 // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1082 // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
1083
1084 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1085 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
d7d7e825 1086
1087 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1088 continue;
1089 }
1090 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1091 continue;
1092 }
1093
d7d7e825 1094 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1095
1096 Double_t massTwoGammaCandidate = 0.;
1097 Double_t widthTwoGammaCandidate = 0.;
1098 Double_t chi2TwoGammaCandidate =10000.;
1099 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
1100 if(twoGammaCandidate->GetNDF()>0){
1101 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1102
6c84d371 1103 if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
d7d7e825 1104
1105 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1106 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1107
1108 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
1109 Double_t rapidity;
1110 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){
1111 rapidity=0;
1112 }
1113 else{
1114 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1115 }
1116
1117 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue; // minimum opening angle to avoid using ghosttracks
1118
1119 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1120 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1121 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1122 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1123 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
1124 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1125 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
1126 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
1127 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1128 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1129 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1130 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
1131 }
1132 }
1133 delete twoGammaCandidate;
1134
1135 //cout << "nearly at the end of processgamma for neutral meson ..." << endl;
1136
1137
1138 }
1139 }
1140}
1141
1142void AliAnalysisTaskGammaConversion::CalculateBackground(){
1143 // see header file for documentation
1144
1145 vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();
1146 vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();
6c84d371 1147
d7d7e825 1148 for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){
1149 AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);
1150 for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){
1151 AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);
1152
1153 AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);
1154
1155 Double_t massBG =0.;
1156 Double_t widthBG = 0.;
1157 Double_t chi2BG =10000.;
1158 backgroundCandidate->GetMass(massBG,widthBG);
1159 if(backgroundCandidate->GetNDF()>0){
1160 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
1161 if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
1162
01b7fdcc 1163 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
1164 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
d7d7e825 1165
1166 Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);
1167
1168 Double_t rapidity;
1169 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
1170 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
1171
1172
1173
1174
1175 if(openingAngleBG < fMinOpeningAngleGhostCut ) continue; // minimum opening angle to avoid using ghosttracks
1176
1177
1178 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
1179 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
01b7fdcc 1180 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
1181 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
d7d7e825 1182 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
01b7fdcc 1183 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
d7d7e825 1184 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
01b7fdcc 1185 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
1186 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
d7d7e825 1187 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
01b7fdcc 1188 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
d7d7e825 1189 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
1190 }
1191 }
1192 delete backgroundCandidate;
1193 //cout << "nearly at the end of background" << endl;
1194
1195 }
1196 }
1197}
1198
1199
1200
1201void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
1202 //ProcessGammasForGammaJetAnalysis
1203
1204 Double_t distIsoMin;
1205
1206 CreateListOfChargedParticles();
1207
1208
6c84d371 1209 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1210 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1211 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 1212 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
d7d7e825 1213
01b7fdcc 1214 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
d7d7e825 1215 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
1216
1217 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
01b7fdcc 1218 CalculateJetCone(gammaIndex);
d7d7e825 1219 }
1220 }
1221 }
1222}
1223
1224void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
1225 // CreateListOfChargedParticles
1226
1227 fESDEvent = fV0Reader->GetESDEvent();
1228 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1229 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1230
1231 if(!curTrack){
1232 continue;
1233 }
1234
1235 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
6c84d371 1236 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
1237 // fChargedParticles.push_back(curTrack);
d7d7e825 1238 fChargedParticlesId.push_back(iTracks);
1239 }
1240 }
1241}
01b7fdcc 1242void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
6c84d371 1243 // CaculateJetCone
d7d7e825 1244
1245 Double_t cone;
1246 Double_t coneSize=0.3;
1247 Double_t ptJet=0;
1248
6c84d371 1249 // AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
1250 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 1251 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
d7d7e825 1252
6c84d371 1253 AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
d7d7e825 1254 Double_t momLeadingCharged[3];
1255 leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
1256
01b7fdcc 1257 TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
d7d7e825 1258
01b7fdcc 1259 Double_t phi1=momentumVectorLeadingCharged.Phi();
1260 Double_t eta1=momentumVectorLeadingCharged.Eta();
1261 Double_t phi3=momentumVectorCurrentGamma.Phi();
d7d7e825 1262
6c84d371 1263 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1264 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 1265 Int_t chId = fChargedParticlesId[iCh];
1266 if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
1267 Double_t mom[3];
1268 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 1269 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1270 Double_t phi2=momentumVectorChargedParticle.Phi();
1271 Double_t eta2=momentumVectorChargedParticle.Eta();
d7d7e825 1272
1273
1274 cone=100.;
1275 if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
1276 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
1277 }else{
1278 if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
1279 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
1280 }
1281 if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
1282 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
1283 }
1284 }
1285
01b7fdcc 1286 if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
1287 ptJet+= momentumVectorChargedParticle.Pt();
1288 Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
1289 Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
d7d7e825 1290 fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
1291 fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
1292
1293 }
1294
1295 Double_t dphiHdrGam=phi3-phi2;
1296 if ( dphiHdrGam < (-TMath::PiOver2())){
1297 dphiHdrGam+=(TMath::TwoPi());
1298 }
1299
1300 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1301 dphiHdrGam-=(TMath::TwoPi());
1302 }
1303
01b7fdcc 1304 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 1305 fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
1306 }
1307 }//track loop
1308
1309
1310}
1311
1312Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
1313 // GetMinimumDistanceToCharge
1314
1315 Double_t fIsoMin=100.;
1316 Double_t ptLeadingCharged=-1.;
1317
6c84d371 1318 AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
01b7fdcc 1319 TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
d7d7e825 1320
01b7fdcc 1321 Double_t phi1=momentumVectorgammaHighestPt.Phi();
1322 Double_t eta1=momentumVectorgammaHighestPt.Eta();
d7d7e825 1323
6c84d371 1324 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1325 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 1326 Int_t chId = fChargedParticlesId[iCh];
1327 if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
1328 Double_t mom[3];
1329 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 1330 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1331 Double_t phi2=momentumVectorChargedParticle.Phi();
1332 Double_t eta2=momentumVectorChargedParticle.Eta();
d7d7e825 1333 Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
1334
01b7fdcc 1335 if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
d7d7e825 1336 if (iso<fIsoMin){
1337 fIsoMin=iso;
1338 }
1339 }
1340
1341 Double_t dphiHdrGam=phi1-phi2;
1342 if ( dphiHdrGam < (-TMath::PiOver2())){
1343 dphiHdrGam+=(TMath::TwoPi());
1344 }
1345
1346 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
1347 dphiHdrGam-=(TMath::TwoPi());
1348 }
01b7fdcc 1349 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 1350 fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
1351 }
1352
1353 if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
01b7fdcc 1354 if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
1355 ptLeadingCharged=momentumVectorChargedParticle.Pt();
d7d7e825 1356 fLeadingChargedIndex=iCh;
1357 }
1358 }
1359
1360 }//track loop
1361 fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
1362 return fIsoMin;
1363
1364}
1365
1366 Int_t AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
1367 //GetIndexHighestPtGamma
1368
1369 Int_t indexHighestPtGamma=-1;
01b7fdcc 1370 //Double_t
1371 fGammaPtHighest = -100.;
d7d7e825 1372
6c84d371 1373 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1374 AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
01b7fdcc 1375 TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
1376 if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
1377 fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
d7d7e825 1378 //gammaHighestPt = gammaHighestPtCandidate;
1379 indexHighestPtGamma=firstGammaIndex;
1380 }
1381 }
1382
1383 return indexHighestPtGamma;
1384
1385}
1386
1387
1388void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
1389{
1390 // Terminate analysis
1391 //
1392 AliDebug(1,"Do nothing in Terminate");
1393}
1394
1395void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
1396{
1397 //AOD
1398 fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
1399 fAODBranch->SetName(fAODBranchName);
1400 AddAODBranch("TClonesArray", &fAODBranch);
1401
1402 // Create the output container
1403 if(fOutputContainer != NULL){
1404 delete fOutputContainer;
1405 fOutputContainer = NULL;
1406 }
1407 if(fOutputContainer == NULL){
1408 fOutputContainer = new TList();
1409 }
1410
1411 //Adding the histograms to the output container
1412 fHistograms->GetOutputContainer(fOutputContainer);
1413
1414
1415 if(fWriteNtuple){
1416 if(fGammaNtuple == NULL){
1417 fGammaNtuple = new TNtuple("V0ntuple","V0ntuple","OnTheFly:HasVertex:NegPIDProb:PosPIDProb:X:Y:Z:R:MotherCandidateNDF:MotherCandidateChi2:MotherCandidateEnergy:MotherCandidateEta:MotherCandidatePt:MotherCandidateMass:MotherCandidateWidth:MCMotherCandidatePT:EPOpeningAngle:ElectronEnergy:ElectronPt:ElectronEta:ElectronPhi:PositronEnergy:PositronPt:PositronEta:PositronPhi:HasSameMCMother:MotherMCParticlePIDCode",50000);
1418 }
1419 if(fNeutralMesonNtuple == NULL){
1420 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
1421 }
1422 TList * ntupleTList = new TList();
1423 ntupleTList->SetName("Ntuple");
1424 ntupleTList->Add((TNtuple*)fGammaNtuple);
1425 fOutputContainer->Add(ntupleTList);
1426 }
1427
1428 fOutputContainer->SetName(GetName());
1429}
1430
1431Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
1432 //helper function
1433 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
1434 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
1435 return v3D0.Angle(v3D1);
1436}
1437
1438void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
1439 // see header file for documentation
1440
1441 vector<Int_t> indexOfGammaParticle;
1442
1443 fStack = fV0Reader->GetMCStack();
1444
1445 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1446 return; // aborts if the primary vertex does not have contributors.
1447 }
1448
1449 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
1450 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1451 if(particle->GetPdgCode()==22){ //Gamma
1452 if(particle->GetNDaughters() >= 2){
1453 TParticle* electron=NULL;
1454 TParticle* positron=NULL;
1455 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1456 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
1457 if(tmpDaughter->GetUniqueID() == 5){
1458 if(tmpDaughter->GetPdgCode() == 11){
1459 electron = tmpDaughter;
1460 }
1461 else if(tmpDaughter->GetPdgCode() == -11){
1462 positron = tmpDaughter;
1463 }
1464 }
1465 }
1466 if(electron!=NULL && positron!=0){
1467 if(electron->R()<160){
1468 indexOfGammaParticle.push_back(iTracks);
1469 }
1470 }
1471 }
1472 }
1473 }
1474
1475 Int_t nFoundGammas=0;
1476 Int_t nNotFoundGammas=0;
1477
1478 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1479 for(Int_t i=0;i<numberOfV0s;i++){
1480 fV0Reader->GetV0(i);
1481
1482 if(fV0Reader->HasSameMCMother() == kFALSE){
1483 continue;
1484 }
1485
1486 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1487 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1488
1489 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1490 continue;
1491 }
1492 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1493 continue;
1494 }
1495
1496 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1497 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
1498 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
1499 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
1500 nFoundGammas++;
1501 }
1502 else{
1503 nNotFoundGammas++;
1504 }
1505 }
1506 }
1507 }
1508 // cout<<"Found: "<<nFoundGammas<<" of: "<<indexOfGammaParticle.size()<<endl;
1509}
1510
1511
1512void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
1513 // see header file for documantation
1514
1515 fESDEvent = fV0Reader->GetESDEvent();
1516
1517
6c84d371 1518 TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
1519 TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
1520 TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
1521 TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
1522 TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
1523 TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
d7d7e825 1524
6c84d371 1525 /*
1526 vector <AliESDtrack*> vESDeNegTemp(0);
1527 vector <AliESDtrack*> vESDePosTemp(0);
1528 vector <AliESDtrack*> vESDxNegTemp(0);
1529 vector <AliESDtrack*> vESDxPosTemp(0);
1530 vector <AliESDtrack*> vESDeNegNoJPsi(0);
1531 vector <AliESDtrack*> vESDePosNoJPsi(0);
1532 */
d7d7e825 1533
1534
1535 fHistograms->FillTable("Table_Electrons",0);//Count number of Events
1536
1537 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1538 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1539
1540 if(!curTrack){
1541 //print warning here
1542 continue;
1543 }
1544
1545 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
1546 double r[3];curTrack->GetConstrainedXYZ(r);
1547
1548 TVector3 rXYZ(r);
1549
1550 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
1551
1552 Bool_t flagKink = kTRUE;
1553 Bool_t flagTPCrefit = kTRUE;
1554 Bool_t flagTRDrefit = kTRUE;
1555 Bool_t flagITSrefit = kTRUE;
1556 Bool_t flagTRDout = kTRUE;
1557 Bool_t flagVertex = kTRUE;
1558
1559
1560 //Cuts ---------------------------------------------------------------
1561
1562 if(curTrack->GetKinkIndex(0) > 0){
1563 fHistograms->FillHistogram("Table_Electrons",5);//Count kink
1564 flagKink = kFALSE;
1565 }
1566
1567 ULong_t trkStatus = curTrack->GetStatus();
1568
1569 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
1570
1571 if(!tpcRefit){
1572 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
1573 flagTPCrefit = kFALSE;
1574 }
1575
1576 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
1577 if(!itsRefit){
1578 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
1579 flagITSrefit = kFALSE;
1580 }
1581
1582 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
1583
1584 if(!trdRefit){
1585 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
1586 flagTRDrefit = kFALSE;
1587 }
1588
1589 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
1590
1591 if(!trdOut) {
1592 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
1593 flagTRDout = kFALSE;
1594 }
1595
1596 double nsigmaToVxt = GetSigmaToVertex(curTrack);
1597
1598 if(nsigmaToVxt > 3){
1599 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
1600 flagVertex = kFALSE;
1601 }
1602
1603 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
1604 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
1605
1606
1607 Stat_t pid, weight;
1608 GetPID(curTrack, pid, weight);
1609
1610 if(pid!=0){
1611 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
1612 }
1613
1614 if(pid == 0){
1615 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
1616 }
1617
1618
1619
1620
1621 Int_t labelMC = TMath::Abs(curTrack->GetLabel());
1622 TParticle* curParticle = fStack->Particle(labelMC);
1623
1624
1625
1626
1627 TLorentzVector curElec;
1628 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
1629
1630
1631
1632
1633 if(curTrack->GetSign() > 0){
1634
6c84d371 1635 // vESDxPosTemp.push_back(curTrack);
1636 new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
d7d7e825 1637
1638 if( pid == 0){
1639
1640 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1641 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
1642 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1643 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1644 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
6c84d371 1645 // vESDePosTemp.push_back(curTrack);
1646 new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
d7d7e825 1647
1648 }
1649
1650 }
1651 else {
6c84d371 1652 // vESDxNegTemp.push_back(curTrack);
1653 if(vESDxNegTemp == NULL){
1654 cout<<"TCloes is zero god damn it"<<endl;
1655 }
1656 if(curTrack == NULL){
1657 cout<<"curTrack is zero god damn it"<<endl;
1658 }
1659
1660 new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
1661
d7d7e825 1662 if( pid == 0){
1663
1664 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1665 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
1666 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1667 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1668 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
6c84d371 1669 //vESDeNegTemp.push_back(curTrack);
1670 new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
d7d7e825 1671
1672 }
6c84d371 1673
d7d7e825 1674 }
1675
1676 }
1677
1678
1679 Bool_t ePosJPsi = kFALSE;
1680 Bool_t eNegJPsi = kFALSE;
1681 Bool_t ePosPi0 = kFALSE;
1682 Bool_t eNegPi0 = kFALSE;
1683
1684 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
1685
6c84d371 1686 for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
1687 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
1688 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
1689 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
d7d7e825 1690 TParticle* partMother = fStack ->Particle(labelMother);
1691 if (partMother->GetPdgCode() == 111){
1692 ieNegPi0 = iNeg;
1693 eNegPi0 = kTRUE;
1694 }
1695 if(partMother->GetPdgCode() == 443){ //Mother JPsi
1696 fHistograms->FillTable("Table_Electrons",14);
1697 ieNegJPsi = iNeg;
1698 eNegJPsi = kTRUE;
1699 }
1700 else{
6c84d371 1701 // vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
1702 new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
d7d7e825 1703 // cout<<"ESD No Positivo JPsi "<<endl;
1704 }
1705
1706 }
1707 }
1708
6c84d371 1709 for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
1710 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
1711 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
1712 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
d7d7e825 1713 TParticle* partMother = fStack ->Particle(labelMother);
1714 if (partMother->GetPdgCode() == 111){
1715 iePosPi0 = iPos;
1716 ePosPi0 = kTRUE;
1717 }
1718 if(partMother->GetPdgCode() == 443){ //Mother JPsi
1719 fHistograms->FillTable("Table_Electrons",15);
1720 iePosJPsi = iPos;
1721 ePosJPsi = kTRUE;
1722 }
1723 else{
6c84d371 1724 // vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
1725 new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
d7d7e825 1726 // cout<<"ESD No Negativo JPsi "<<endl;
1727 }
1728
1729 }
1730 }
1731
1732 if( eNegJPsi && ePosJPsi ){
1733 TVector3 tempeNegV,tempePosV;
6c84d371 1734 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
1735 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
d7d7e825 1736 fHistograms->FillTable("Table_Electrons",16);
1737 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
6c84d371 1738 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
1739 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
d7d7e825 1740 }
1741
1742 if( eNegPi0 && ePosPi0 ){
1743 TVector3 tempeNegV,tempePosV;
6c84d371 1744 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
1745 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
d7d7e825 1746 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
6c84d371 1747 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
1748 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
d7d7e825 1749 }
1750
1751
1752 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
1753
6c84d371 1754 CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
d7d7e825 1755
6c84d371 1756 // vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
1757 // vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
1758
1759 TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
1760 TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
d7d7e825 1761
1762
1763 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
1764
1765
1766
1767
1768 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
1769
1770
1771 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
1772 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
1773
1774
1775
6c84d371 1776 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
d7d7e825 1777
1778 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
6c84d371 1779 *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
d7d7e825 1780
1781 //BackGround
1782
1783 //Like Sign e+e-
1784 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
1785 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
1786 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
1787 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
1788
1789 // Like Sign e+e- no JPsi
1790 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
1791 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
1792
1793 //Mixed Event
1794
6c84d371 1795 if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
d7d7e825 1796 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
6c84d371 1797 *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
1798 *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
1799 *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
d7d7e825 1800
1801 }
1802
1803 /*
1804 //Photons P
1805 Double_t vtx[3];
1806 vtx[0]=0;vtx[1]=0;vtx[2]=0;
1807 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
1808
1809 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
1810
1811
1812
1813 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
1814 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
1815
1816 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
1817
1818 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
1819
1820 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
1821
1822
1823 }
1824
1825
1826 */
1827
1828
1829}
1830
6c84d371 1831/*
1832 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
d7d7e825 1833 //see header file for documentation
1834 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
6c84d371 1835 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
1836 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
1837 }
1838 }
1839 }
1840*/
1841void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
1842 //see header file for documentation
1843 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
1844 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
1845 fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
d7d7e825 1846 }
1847 }
1848}
6c84d371 1849void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
d7d7e825 1850 //see header file for documentation
6c84d371 1851 for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
1852 TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
1853 for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
1854 TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
d7d7e825 1855 TLorentzVector np = ep + en;
1856 fHistograms->FillHistogram(histoName.Data(),np.M());
1857 }
1858 }
d7d7e825 1859}
1860
6c84d371 1861void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
1862 TClonesArray const tlVeNeg,TClonesArray const tlVePos)
d7d7e825 1863{
1864 //see header file for documentation
1865
6c84d371 1866 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
d7d7e825 1867
6c84d371 1868 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
d7d7e825 1869
6c84d371 1870 TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
d7d7e825 1871
6c84d371 1872 for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
d7d7e825 1873
6c84d371 1874 // AliKFParticle * gammaCandidate = &fKFGammas[iGam];
1875 AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
d7d7e825 1876 TLorentzVector g;
1877
1878 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
1879 TLorentzVector xyg = xy + g;
1880 fHistograms->FillHistogram(histoMass.Data(),xyg.M());
1881 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
1882 }
1883 }
1884 }
1885
1886}
6c84d371 1887void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
d7d7e825 1888{
1889 // see header file for documentation
6c84d371 1890 for(Int_t i=0; i < e.GetEntriesFast(); i++)
d7d7e825 1891 {
6c84d371 1892 for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
d7d7e825 1893 {
6c84d371 1894 TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
d7d7e825 1895
1896 fHistograms->FillHistogram(hBg.Data(),ee.M());
1897 }
1898 }
1899}
1900
1901
6c84d371 1902void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
1903 TClonesArray const positiveElectrons,
1904 TClonesArray const gammas){
d7d7e825 1905 // see header file for documentation
1906
6c84d371 1907 UInt_t sizeN = negativeElectrons.GetEntriesFast();
1908 UInt_t sizeP = positiveElectrons.GetEntriesFast();
1909 UInt_t sizeG = gammas.GetEntriesFast();
d7d7e825 1910
1911
1912
1913 vector <Bool_t> xNegBand(sizeN);
1914 vector <Bool_t> xPosBand(sizeP);
1915 vector <Bool_t> gammaBand(sizeG);
1916
1917
1918 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
1919 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
1920 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
1921
1922
1923 for(UInt_t iPos=0; iPos < sizeP; iPos++){
1924
6c84d371 1925 Double_t aP[3];
1926 ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
d7d7e825 1927
1928 TVector3 ePosV(aP[0],aP[1],aP[2]);
1929
1930 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
1931
6c84d371 1932 Double_t aN[3];
1933 ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
d7d7e825 1934 TVector3 eNegV(aN[0],aN[1],aN[2]);
1935
1936 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
1937 xPosBand[iPos]=kFALSE;
1938 xNegBand[iNeg]=kFALSE;
1939 }
1940
1941 for(UInt_t iGam=0; iGam < sizeG; iGam++){
6c84d371 1942 AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
d7d7e825 1943 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
1944 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
1945 gammaBand[iGam]=kFALSE;
1946 }
1947 }
1948 }
1949
1950
1951
1952
1953 for(UInt_t iPos=0; iPos < sizeP; iPos++){
1954 if(xPosBand[iPos]){
6c84d371 1955 new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
1956 // fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
d7d7e825 1957 }
1958 }
1959 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
1960 if(xNegBand[iNeg]){
6c84d371 1961 new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
1962 // fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
d7d7e825 1963 }
1964 }
1965 for(UInt_t iGam=0; iGam < sizeG; iGam++){
1966 if(gammaBand[iGam]){
6c84d371 1967 new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
1968 //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
d7d7e825 1969 }
1970 }
1971}
1972
1973
1974void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
1975{
1976 // see header file for documentation
1977 pid = -1;
1978 weight = -1;
1979
1980 double wpart[5];
1981 double wpartbayes[5];
1982
1983 //get probability of the diffenrent particle types
1984 track->GetESDpid(wpart);
1985
1986 // Tentative particle type "concentrations"
1987 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
1988
1989 //Bayes' formula
1990 double rcc = 0.;
1991 for (int i = 0; i < 5; i++)
1992 {
1993 rcc+=(c[i] * wpart[i]);
1994 }
1995
1996
1997
1998 for (int i=0; i<5; i++) {
1999 if( rcc!=0){
2000 wpartbayes[i] = c[i] * wpart[i] / rcc;
2001 }
2002 }
2003
2004
2005
2006 Float_t max=0.;
2007 int ipid=-1;
2008 //find most probable particle in ESD pid
2009 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
2010 for (int i = 0; i < 5; i++)
2011 {
2012 if (wpartbayes[i] > max)
2013 {
2014 ipid = i;
2015 max = wpartbayes[i];
2016 }
2017 }
2018
2019 pid = ipid;
2020 weight = max;
2021}
2022double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
2023{
2024 // Calculates the number of sigma to the vertex.
2025
2026 Float_t b[2];
2027 Float_t bRes[2];
2028 Float_t bCov[3];
2029 t->GetImpactParameters(b,bCov);
2030 if (bCov[0]<=0 || bCov[2]<=0) {
2031 AliDebug(1, "Estimated b resolution lower or equal zero!");
2032 bCov[0]=0; bCov[2]=0;
2033 }
2034 bRes[0] = TMath::Sqrt(bCov[0]);
2035 bRes[1] = TMath::Sqrt(bCov[2]);
2036
2037 // -----------------------------------
2038 // How to get to a n-sigma cut?
2039 //
2040 // The accumulated statistics from 0 to d is
2041 //
2042 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
2043 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
2044 //
2045 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
2046 // Can this be expressed in a different way?
2047
2048 if (bRes[0] == 0 || bRes[1] ==0)
2049 return -1;
2050
2051 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
2052
2053 // stupid rounding problem screws up everything:
2054 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
2055 if (TMath::Exp(-d * d / 2) < 1e-10)
2056 return 1000;
2057
2058
2059 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
2060 return d;
2061}
6c84d371 2062
2063//vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
2064TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
01b7fdcc 2065 //Return TLoresntz vector of track?
6c84d371 2066 // vector <TLorentzVector> tlVtrack(0);
2067 TClonesArray array("TLorentzVector",0);
d7d7e825 2068
6c84d371 2069 for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
2070 double p[3];
2071 //esdTrack[itrack]->GetConstrainedPxPyPz(p);
2072 ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
d7d7e825 2073 TLorentzVector currentTrack;
01b7fdcc 2074 currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
6c84d371 2075 new((array)[array.GetEntriesFast()]) TLorentzVector(currentTrack);
2076 // tlVtrack.push_back(currentTrack);
d7d7e825 2077 }
6c84d371 2078
2079 return array;
d7d7e825 2080
6c84d371 2081 // return tlVtrack;
d7d7e825 2082}
2083