]> 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
Add histo for Event Quality and pi0 spectra for different types of process, event...
[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"
4a6157dc 31//#include "AliCFManager.h" // for CF
32//#include "AliCFContainer.h" // for CF
9c1cb6f7 33#include "AliESDInputHandler.h"
34#include "AliAnalysisManager.h"
4a6157dc 35#include "AliGammaConversionAODObject.h"
5e55d806 36#include "AliGammaConversionBGHandler.h"
6272370b 37#include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
9c1cb6f7 38#include "AliKFVertex.h"
e60f3265 39#include <AliMCEventHandler.h>
d707e3cf 40class AliESDTrackCuts;
4a6157dc 41class AliCFContainer;
42class AliCFManager;
d7d7e825 43class AliKFVertex;
44class AliAODHandler;
45class AliAODEvent;
46class ALiESDEvent;
47class AliMCEvent;
48class AliMCEventHandler;
49class AliESDInputHandler;
50class AliAnalysisManager;
51class Riostream;
52class TFile;
53class TInterpreter;
54class TSystem;
55class TROOT;
56
57ClassImp(AliAnalysisTaskGammaConversion)
58
59
60AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
61AliAnalysisTaskSE(),
62 fV0Reader(NULL),
63 fStack(NULL),
a0b94e5c 64 fMCTruth(NULL), // for CF
4a6157dc 65 fGCMCEvent(NULL), // for CF
d7d7e825 66 fESDEvent(NULL),
67 fOutputContainer(NULL),
a0b94e5c 68 fCFManager(0x0), // for CF
d7d7e825 69 fHistograms(NULL),
b5832f95 70 fTriggerCINT1B(kFALSE),
d7d7e825 71 fDoMCTruth(kFALSE),
72 fDoNeutralMeson(kFALSE),
6272370b 73 fDoOmegaMeson(kFALSE),
d7d7e825 74 fDoJet(kFALSE),
75 fDoChic(kFALSE),
9c1cb6f7 76 fRecalculateV0ForGamma(kFALSE),
6c84d371 77 fKFReconstructedGammasTClone(NULL),
6272370b 78 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 79 fKFRecalculatedGammasTClone(NULL),
6c84d371 80 fCurrentEventPosElectronTClone(NULL),
81 fCurrentEventNegElectronTClone(NULL),
82 fKFReconstructedGammasCutTClone(NULL),
83 fPreviousEventTLVNegElectronTClone(NULL),
84 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 85 fElectronv1(),
86 fElectronv2(),
6272370b 87 fGammav1(),
88 fGammav2(),
9c1cb6f7 89 fElectronRecalculatedv1(),
90 fElectronRecalculatedv2(),
d7d7e825 91 fElectronMass(-1),
92 fGammaMass(-1),
93 fPi0Mass(-1),
94 fEtaMass(-1),
95 fGammaWidth(-1),
96 fPi0Width(-1),
97 fEtaWidth(-1),
98 fMinOpeningAngleGhostCut(0.),
9640a3d1 99 fEsdTrackCuts(NULL),
d7d7e825 100 fCalculateBackground(kFALSE),
101 fWriteNtuple(kFALSE),
102 fGammaNtuple(NULL),
103 fNeutralMesonNtuple(NULL),
104 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 105 fChargedParticles(NULL),
d7d7e825 106 fChargedParticlesId(),
107 fGammaPtHighest(0.),
108 fMinPtForGammaJet(1.),
109 fMinIsoConeSize(0.2),
110 fMinPtIsoCone(0.7),
111 fMinPtGamChargedCorr(0.5),
112 fMinPtJetCone(0.5),
113 fLeadingChargedIndex(-1),
9640a3d1 114 fLowPtMapping(1.),
115 fHighPtMapping(3.),
1e7846f4 116 fDoCF(kFALSE),
04bf4381 117 fAODGamma(NULL),
118 fAODPi0(NULL),
119 fAODOmega(NULL),
037dc2db 120 fAODBranchName("GammaConv"),
d765d400 121 fKFForceAOD(kFALSE),
332f1f44 122 fKFDeltaAODFileName(""),
037dc2db 123 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 124 fUseTrackMultiplicityForBG(kTRUE),
125 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 126 fApplyChi2Cut(kFALSE),
111d75df 127 nRandomEventsForBG(15),
128 nDegreesPMBackground(15),
129 fDoRotation(kTRUE),
130 fCheckBGProbability(kTRUE),
037dc2db 131 fKFReconstructedGammasV0Index()
d7d7e825 132{
133 // Default constructor
5a34881d 134
135 /* Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
d7d7e825 136 // Common I/O in slot 0
137 DefineInput (0, TChain::Class());
138 DefineOutput(0, TTree::Class());
139
140 // Your private output
141 DefineOutput(1, TList::Class());
a0b94e5c 142
d7d7e825 143 // Define standard ESD track cuts for Gamma-hadron correlation
144 SetESDtrackCuts();
5a34881d 145 */
d7d7e825 146}
147
148AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
149 AliAnalysisTaskSE(name),
150 fV0Reader(NULL),
151 fStack(NULL),
a0b94e5c 152 fMCTruth(NULL), // for CF
4a6157dc 153 fGCMCEvent(NULL), // for CF
d7d7e825 154 fESDEvent(NULL),
155 fOutputContainer(0x0),
a0b94e5c 156 fCFManager(0x0), // for CF
d7d7e825 157 fHistograms(NULL),
b5832f95 158 fTriggerCINT1B(kFALSE),
d7d7e825 159 fDoMCTruth(kFALSE),
160 fDoNeutralMeson(kFALSE),
6272370b 161 fDoOmegaMeson(kFALSE),
d7d7e825 162 fDoJet(kFALSE),
163 fDoChic(kFALSE),
9c1cb6f7 164 fRecalculateV0ForGamma(kFALSE),
6c84d371 165 fKFReconstructedGammasTClone(NULL),
6272370b 166 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 167 fKFRecalculatedGammasTClone(NULL),
6c84d371 168 fCurrentEventPosElectronTClone(NULL),
169 fCurrentEventNegElectronTClone(NULL),
170 fKFReconstructedGammasCutTClone(NULL),
171 fPreviousEventTLVNegElectronTClone(NULL),
172 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 173 fElectronv1(),
174 fElectronv2(),
6272370b 175 fGammav1(),
176 fGammav2(),
9c1cb6f7 177 fElectronRecalculatedv1(),
178 fElectronRecalculatedv2(),
d7d7e825 179 fElectronMass(-1),
180 fGammaMass(-1),
181 fPi0Mass(-1),
182 fEtaMass(-1),
183 fGammaWidth(-1),
184 fPi0Width(-1),
185 fEtaWidth(-1),
186 fMinOpeningAngleGhostCut(0.),
9640a3d1 187 fEsdTrackCuts(NULL),
d7d7e825 188 fCalculateBackground(kFALSE),
189 fWriteNtuple(kFALSE),
190 fGammaNtuple(NULL),
191 fNeutralMesonNtuple(NULL),
192 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 193 fChargedParticles(NULL),
d7d7e825 194 fChargedParticlesId(),
195 fGammaPtHighest(0.),
196 fMinPtForGammaJet(1.),
197 fMinIsoConeSize(0.2),
198 fMinPtIsoCone(0.7),
199 fMinPtGamChargedCorr(0.5),
200 fMinPtJetCone(0.5),
201 fLeadingChargedIndex(-1),
9640a3d1 202 fLowPtMapping(1.),
203 fHighPtMapping(3.),
c59360eb 204 fDoCF(kFALSE),
205 fAODGamma(NULL),
206 fAODPi0(NULL),
207 fAODOmega(NULL),
208 fAODBranchName("GammaConv"),
209 fKFForceAOD(kFALSE),
210 fKFDeltaAODFileName(""),
211 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 212 fUseTrackMultiplicityForBG(kTRUE),
213 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 214 fApplyChi2Cut(kFALSE),
111d75df 215 nRandomEventsForBG(15),
216 nDegreesPMBackground(15),
217 fDoRotation(kTRUE),
218 fCheckBGProbability(kTRUE),
037dc2db 219 fKFReconstructedGammasV0Index()
d7d7e825 220{
221 // Common I/O in slot 0
222 DefineInput (0, TChain::Class());
223 DefineOutput(0, TTree::Class());
224
225 // Your private output
226 DefineOutput(1, TList::Class());
a0b94e5c 227 DefineOutput(2, AliCFContainer::Class()); // for CF
228
229
d7d7e825 230 // Define standard ESD track cuts for Gamma-hadron correlation
231 SetESDtrackCuts();
9c1cb6f7 232
d7d7e825 233}
234
235AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
236{
237 // Remove all pointers
238
239 if(fOutputContainer){
240 fOutputContainer->Clear() ;
241 delete fOutputContainer ;
242 }
243 if(fHistograms){
244 delete fHistograms;
245 }
246 if(fV0Reader){
247 delete fV0Reader;
248 }
a0b94e5c 249
250 // for CF
251 if(fCFManager){
252 delete fCFManager;
253 }
9640a3d1 254
255 if(fEsdTrackCuts){
256 delete fEsdTrackCuts;
257 }
258
04bf4381 259 //Delete AODs
260 if (fAODGamma) {
261 fAODGamma->Clear();
262 delete fAODGamma;
d7d7e825 263 }
04bf4381 264 fAODGamma = NULL;
265
266 if (fAODPi0) {
267 fAODPi0->Clear();
268 delete fAODPi0;
269 }
270 fAODPi0 = NULL;
d7d7e825 271
04bf4381 272 if (fAODOmega) {
273 fAODOmega->Clear();
274 delete fAODOmega;
275 }
276 fAODOmega = NULL;
277
04bf4381 278}
d7d7e825 279
12464034 280
d7d7e825 281void AliAnalysisTaskGammaConversion::Init()
282{
283 // Initialization
4a6157dc 284 // AliLog::SetGlobalLogLevel(AliLog::kError);
d7d7e825 285}
286void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
287{
288 // SetESDtrackCuts
9640a3d1 289 if (fEsdTrackCuts!=NULL){
290 delete fEsdTrackCuts;
291 }
d7d7e825 292 fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
a0b94e5c 293 //standard cuts from:
294 //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
ebcfaa7e 295
296 // Cuts used up to 3rd of March
297
298 // fEsdTrackCuts->SetMinNClustersTPC(50);
299// fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
300// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
301// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
302// fEsdTrackCuts->SetMaxNsigmaToVertex(3);
303// fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
304
305 //------- To be tested-----------
d707e3cf 306 // Cuts used up to 26th of Agost
307// Int_t minNClustersTPC = 70;
308// Double_t maxChi2PerClusterTPC = 4.0;
309// Double_t maxDCAtoVertexXY = 2.4; // cm
310// Double_t maxDCAtoVertexZ = 3.2; // cm
311// fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
312// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
313// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
314// // fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
315// fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
316// fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
317// fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
318// fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
319// fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
320// fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
321// fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
322// fEsdTrackCuts->SetPtRange(0.15);
323
324// fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
325
326
327// Using standard function for setting Cuts
328 Bool_t selectPrimaries=kTRUE;
329 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
330 fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
331 fEsdTrackCuts->SetPtRange(0.15);
ebcfaa7e 332
333 //----- From Jacek 10.03.03 ------------------/
334// minNClustersTPC = 70;
335// maxChi2PerClusterTPC = 4.0;
336// maxDCAtoVertexXY = 2.4; // cm
337// maxDCAtoVertexZ = 3.2; // cm
338
339// esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
340// esdTrackCuts->SetRequireTPCRefit(kFALSE);
341// esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
342// esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
343// esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
344// esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
345// esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
346// esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
347// esdTrackCuts->SetDCAToVertex2D(kTRUE);
348
349
350
d7d7e825 351 // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
037dc2db 352 // fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
d7d7e825 353}
354
c00009fb 355void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
d7d7e825 356{
357 // Execute analysis for current event
6272370b 358
9c1cb6f7 359 // Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
360
361 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
362 AliESDInputHandler *esdHandler=0x0;
363 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
364 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
365 } else {
366 //load esd pid bethe bloch parameters depending on the existance of the MC handler
367 // yes: MC parameters
368 // no: data parameters
369 if (!AliV0Reader::GetESDpid()){
370 if (fMCEvent ) {
371 AliV0Reader::InitESDpid();
372 } else {
373 AliV0Reader::InitESDpid(1);
374 }
375 }
376 }
377
d765d400 378 //Must set fForceAOD to true for the AOD to get filled. Should only be done when running independent chain / train.
379 if(fKFForceAOD) {
380 if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
381 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
382 }
48682642 383
384 if(fV0Reader == NULL){
385 // Write warning here cuts and so on are default if this ever happens
6272370b 386 }
48682642 387
e60f3265 388 if (fMCEvent ) {
389 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
390
391 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
392 if (!mcHandler){ AliError("Could not retrive MC event handler!"); return; }
393 if (!mcHandler->InitOk() ) return;
394 if (!mcHandler->TreeK() ) return;
395 if (!mcHandler->TreeTR() ) return;
396 }
397
48682642 398 fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
399
400 fV0Reader->Initialize();
401 fDoMCTruth = fV0Reader->GetDoMCTruth();
402
12464034 403 if(fAODGamma) fAODGamma->Delete();
404 if(fAODPi0) fAODPi0->Delete();
405 if(fAODOmega) fAODOmega->Delete();
a0b94e5c 406
6c84d371 407 if(fKFReconstructedGammasTClone == NULL){
408 fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
409 }
410 if(fCurrentEventPosElectronTClone == NULL){
411 fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
412 }
413 if(fCurrentEventNegElectronTClone == NULL){
414 fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
415 }
416 if(fKFReconstructedGammasCutTClone == NULL){
417 fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
418 }
419 if(fPreviousEventTLVNegElectronTClone == NULL){
420 fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
421 }
422 if(fPreviousEventTLVPosElectronTClone == NULL){
423 fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
424 }
425 if(fChargedParticles == NULL){
426 fChargedParticles = new TClonesArray("AliESDtrack",0);
427 }
6272370b 428
429 if(fKFReconstructedPi0sTClone == NULL){
430 fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
431 }
9c1cb6f7 432
433 if(fKFRecalculatedGammasTClone == NULL){
434 fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
435 }
436
a0b94e5c 437
6c84d371 438 //clear TClones
1e7846f4 439 fKFReconstructedGammasTClone->Delete();
440 fCurrentEventPosElectronTClone->Delete();
441 fCurrentEventNegElectronTClone->Delete();
442 fKFReconstructedGammasCutTClone->Delete();
443 fPreviousEventTLVNegElectronTClone->Delete();
444 fPreviousEventTLVPosElectronTClone->Delete();
6272370b 445 fKFReconstructedPi0sTClone->Delete();
9c1cb6f7 446 fKFRecalculatedGammasTClone->Delete();
5e55d806 447
d7d7e825 448 //clear vectors
d7d7e825 449 fElectronv1.clear();
450 fElectronv2.clear();
6272370b 451
452 fGammav1.clear();
453 fGammav2.clear();
454
9c1cb6f7 455 fElectronRecalculatedv1.clear();
456 fElectronRecalculatedv2.clear();
457
d7d7e825 458
1e7846f4 459 fChargedParticles->Delete();
5e55d806 460
d7d7e825 461 fChargedParticlesId.clear();
037dc2db 462
463 fKFReconstructedGammasV0Index.clear();
a0b94e5c 464
d7d7e825 465 //Clear the data in the v0Reader
5e55d806 466 // fV0Reader->UpdateEventByEventData();
b5832f95 467
468 //Take Only events with proper trigger
5e55d806 469 /*
b5832f95 470 if(fTriggerCINT1B){
471 if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
472 }
5e55d806 473 */
10e3319b 474
c8206114 475 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
476 // cout<< "Event not taken"<< endl;
477 return; // aborts if the primary vertex does not have contributors.
478 }
479
480
885114d1 481 if(!fV0Reader->CheckForPrimaryVertexZ() ){
482 return;
483 }
484
485
486
d7d7e825 487 // Process the MC information
488 if(fDoMCTruth){
489 ProcessMCData();
490 }
a0b94e5c 491
d7d7e825 492 //Process the v0 information with no cuts
493 ProcessV0sNoCut();
cb90a330 494
d7d7e825 495 // Process the v0 information
496 ProcessV0s();
cb90a330 497
5e55d806 498
d7d7e825 499 //Fill Gamma AOD
500 FillAODWithConversionGammas() ;
a0b94e5c 501
d7d7e825 502 // Process reconstructed gammas
503 if(fDoNeutralMeson == kTRUE){
a0b94e5c 504 ProcessGammasForNeutralMesonAnalysis();
77880bd8 505
48682642 506 }
507
77880bd8 508 if(fDoMCTruth == kTRUE){
509 CheckV0Efficiency();
510 }
d7d7e825 511 //Process reconstructed gammas electrons for Chi_c Analysis
6c84d371 512 if(fDoChic == kTRUE){
d7d7e825 513 ProcessGammaElectronsForChicAnalysis();
514 }
515 // Process reconstructed gammas for gamma Jet/hadron correlations
516 if(fDoJet == kTRUE){
517 ProcessGammasForGammaJetAnalysis();
518 }
a0b94e5c 519
5e55d806 520 //calculate background if flag is set
521 if(fCalculateBackground){
522 CalculateBackground();
523 }
48682642 524
525 if(fDoNeutralMeson == kTRUE){
526// ProcessConvPHOSGammasForNeutralMesonAnalysis();
527 if(fDoOmegaMeson == kTRUE){
528 ProcessGammasForOmegaMesonAnalysis();
529 }
530 }
531
5e55d806 532 //Clear the data in the v0Reader
533 fV0Reader->UpdateEventByEventData();
9c1cb6f7 534 if(fRecalculateV0ForGamma==kTRUE){
535 RecalculateV0ForGamma();
536 }
d7d7e825 537 PostData(1, fOutputContainer);
a0b94e5c 538 PostData(2, fCFManager->GetParticleContainer()); // for CF
d7d7e825 539
540}
541
48682642 542// void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
543// // see header file for documentation
544// // printf(" ConnectInputData %s\n", GetName());
8a685cf3 545
48682642 546// AliAnalysisTaskSE::ConnectInputData(option);
8a685cf3 547
48682642 548// if(fV0Reader == NULL){
549// // Write warning here cuts and so on are default if this ever happens
550// }
551// fV0Reader->Initialize();
552// fDoMCTruth = fV0Reader->GetDoMCTruth();
553// }
d7d7e825 554
555
556
557void AliAnalysisTaskGammaConversion::ProcessMCData(){
558 // see header file for documentation
48682642 559 //InputEvent(), MCEvent());
560 /* TestAnaMarin
d7d7e825 561 fStack = fV0Reader->GetMCStack();
a0b94e5c 562 fMCTruth = fV0Reader->GetMCTruth(); // for CF
4a6157dc 563 fGCMCEvent = fV0Reader->GetMCEvent(); // for CF
48682642 564 */
565 fStack= MCEvent()->Stack();
566 fGCMCEvent=MCEvent();
a0b94e5c 567
568 // for CF
1e7846f4 569 Double_t containerInput[3];
570 if(fDoCF){
571 if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
572 fCFManager->SetEventInfo(fGCMCEvent);
573 }
a0b94e5c 574 // end for CF
575
576
d7d7e825 577 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
578 return; // aborts if the primary vertex does not have contributors.
579 }
bd6d9fa3 580
d7d7e825 581 for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
111d75df 582 //for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
d7d7e825 583 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
bd6d9fa3 584
111d75df 585
bd6d9fa3 586
d7d7e825 587 if (!particle) {
588 //print warning here
589 continue;
590 }
a0b94e5c 591
d7d7e825 592 ///////////////////////Begin Chic Analysis/////////////////////////////
a0b94e5c 593
a68437fb 594 if(particle->GetPdgCode() == 443){//Is JPsi
d7d7e825 595 if(particle->GetNDaughters()==2){
596 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
597 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
63e16c52 598
d7d7e825 599 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
600 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
601 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
602 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
a0b94e5c 603
d7d7e825 604 if( TMath::Abs(daug0->Eta()) < 0.9){
605 if(daug0->GetPdgCode() == -11)
606 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
607 else
608 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
a0b94e5c 609
d7d7e825 610 }
611 if(TMath::Abs(daug1->Eta()) < 0.9){
612 if(daug1->GetPdgCode() == -11)
613 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
614 else
615 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
616 }
617 }
618 }
619 }
620 // const int CHI_C0 = 10441;
621 // const int CHI_C1 = 20443;
622 // const int CHI_C2 = 445
623 if(particle->GetPdgCode() == 22){//gamma from JPsi
624 if(particle->GetMother(0) > -1){
625 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
626 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
627 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
628 if(TMath::Abs(particle->Eta()) < 1.2)
629 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
630 }
631 }
632 }
633 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
634 if( particle->GetNDaughters() == 2){
635 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
636 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
a0b94e5c 637
d7d7e825 638 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
639 if( daug0->GetPdgCode() == 443){
640 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
641 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
642 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
643 fHistograms->FillTable("Table_Electrons",18);
a0b94e5c 644
d7d7e825 645 }//if
646 else if (daug1->GetPdgCode() == 443){
647 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
648 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
649 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
650 fHistograms->FillTable("Table_Electrons",18);
651 }//else if
652 }//gamma o Jpsi
653 }//GetNDaughters
654 }
a0b94e5c 655
656
d7d7e825 657 /////////////////////End Chic Analysis////////////////////////////
a0b94e5c 658
659
d7d7e825 660 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
a0b94e5c 661
d7d7e825 662 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
663
664 Double_t tmpPhi=particle->Phi();
665
666 if(particle->Phi()> TMath::Pi()){
667 tmpPhi = particle->Phi()-(2*TMath::Pi());
668 }
669
670 Double_t rapidity;
671 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
672 rapidity=0;
673 }
674 else{
675 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
676 }
677
111d75df 678
679
680 if(iTracks<=fStack->GetNprimary() ){
681 if ( particle->GetPdgCode()== -211 || particle->GetPdgCode()== 211 ||
682 particle->GetPdgCode()== 2212 || particle->GetPdgCode()==-2212 ||
683 particle->GetPdgCode()== 321 || particle->GetPdgCode()==-321 ){
684
685 // fHistograms->FillHistogram("MC_PhysicalPrimaryCharged_Pt", mcTrack->Pt());
686 }
bd6d9fa3 687 }
111d75df 688
bd6d9fa3 689
690
d7d7e825 691 //process the gammas
692 if (particle->GetPdgCode() == 22){
a607c218 693
694
d7d7e825 695 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
696 continue; // no photon as mothers!
697 }
a607c218 698
d7d7e825 699 if(particle->GetMother(0) >= fStack->GetNprimary()){
700 continue; // the gamma has a mother, and it is not a primary particle
701 }
a0b94e5c 702
d7d7e825 703 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
704 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
705 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
706 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
707 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
708
a0b94e5c 709 // for CF
1e7846f4 710 if(fDoCF){
711 containerInput[0] = particle->Pt();
712 containerInput[1] = particle->Eta();
713 if(particle->GetMother(0) >=0){
714 containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
715 }
716 else{
717 containerInput[2]=-1;
718 }
a607c218 719
1e7846f4 720 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); // generated gamma
721 }
d7d7e825 722 if(particle->GetMother(0) < 0){ // direct gamma
723 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
724 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
725 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
726 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
727 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
728 }
729
d7d7e825 730 // looking for conversion (electron + positron from pairbuilding (= 5) )
731 TParticle* ePos = NULL;
732 TParticle* eNeg = NULL;
733
734 if(particle->GetNDaughters() >= 2){
735 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
736 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
737 if(tmpDaughter->GetUniqueID() == 5){
738 if(tmpDaughter->GetPdgCode() == 11){
739 eNeg = tmpDaughter;
740 }
741 else if(tmpDaughter->GetPdgCode() == -11){
742 ePos = tmpDaughter;
743 }
744 }
745 }
746 }
747
748
749 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
750 continue;
751 }
752
753
754 Double_t ePosPhi = ePos->Phi();
755 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
756
757 Double_t eNegPhi = eNeg->Phi();
758 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
759
760
761 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
762 continue; // no reconstruction below the Pt cut
763 }
a0b94e5c 764
d7d7e825 765 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
766 continue;
767 }
a0b94e5c 768
d7d7e825 769 if(ePos->R()>fV0Reader->GetMaxRCut()){
770 continue; // cuts on distance from collision point
771 }
a0b94e5c 772
773 if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
d7d7e825 774 continue; // outside material
775 }
776
a0b94e5c 777
d7d7e825 778 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
779 continue; // line cut to exclude regions where we do not reconstruct
780 }
781
a0b94e5c 782
783 // for CF
1e7846f4 784 if(fDoCF){
785 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable); // reconstructable gamma
786 }
d7d7e825 787 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
788 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
789 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
790 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
791 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
792 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
793
794 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
795 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
796 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
797 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
798
799 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
800 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
801 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
802 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
803
804
d7d7e825 805 // begin Mapping
806 Int_t rBin = fHistograms->GetRBin(ePos->R());
9640a3d1 807 Int_t zBin = fHistograms->GetZBin(ePos->Vz());
d7d7e825 808 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
9c1cb6f7 809 Double_t rFMD=30;
9640a3d1 810
811 TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
812
d7d7e825 813 TString nameMCMappingPhiR="";
e158cbc3 814 nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 815 // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
d7d7e825 816
817 TString nameMCMappingPhi="";
e158cbc3 818 nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
9640a3d1 819 // fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
ebcfaa7e 820 //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
d7d7e825 821
822 TString nameMCMappingR="";
e158cbc3 823 nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
9640a3d1 824 // fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
ebcfaa7e 825 //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
d7d7e825 826
827 TString nameMCMappingPhiInR="";
e158cbc3 828 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 829 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
830 fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
831
832 TString nameMCMappingZInR="";
833 nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
834 fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
835
836
837 TString nameMCMappingPhiInZ="";
838 nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
839 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
840 fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
841
9c1cb6f7 842
843 if(ePos->R()<rFMD){
844 TString nameMCMappingFMDPhiInZ="";
845 nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
846 fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
847 }
848
9640a3d1 849 TString nameMCMappingRInZ="";
850 nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
851 fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
852
853 if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
854 TString nameMCMappingMidPtPhiInR="";
855 nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
856 fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
857
858 TString nameMCMappingMidPtZInR="";
859 nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
860 fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
861
862
863 TString nameMCMappingMidPtPhiInZ="";
864 nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
865 fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
9c1cb6f7 866
867
868 if(ePos->R()<rFMD){
869 TString nameMCMappingMidPtFMDPhiInZ="";
870 nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
871 fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
872 }
873
9640a3d1 874
875 TString nameMCMappingMidPtRInZ="";
876 nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
877 fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
878
879 }
880
d7d7e825 881 //end mapping
882
883 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
884 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
885 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
886 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
9640a3d1 887 fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
888 fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
889
890
d7d7e825 891 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
892 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
893 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
894 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
895 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
896 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
897
898 } // end direct gamma
899 else{ // mother exits
6c84d371 900 /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
a0b94e5c 901 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
902 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
903 ){
904 fMCGammaChic.push_back(particle);
905 }
6c84d371 906 */
d7d7e825 907 } // end if mother exits
908 } // end if particle is a photon
909
910
911
912 // process motherparticles (2 gammas as daughters)
913 // the motherparticle had already to pass the R and the eta cut, but no line cut.
914 // the line cut is just valid for the conversions!
915
916 if(particle->GetNDaughters() == 2){
917
918 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
919 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
920
921 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
922
923 // Check the acceptance for both gammas
01b7fdcc 924 Bool_t gammaEtaCut = kTRUE;
925 if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
926 Bool_t gammaRCut = kTRUE;
927 if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
d7d7e825 928
929
930
931 // check for conversions now -> have to pass eta, R and line cut!
932 Bool_t daughter0Electron = kFALSE;
933 Bool_t daughter0Positron = kFALSE;
934 Bool_t daughter1Electron = kFALSE;
935 Bool_t daughter1Positron = kFALSE;
936
937 if(daughter0->GetNDaughters() >= 2){ // first gamma
938 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
939 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
940 if(tmpDaughter->GetUniqueID() == 5){
941 if(tmpDaughter->GetPdgCode() == 11){
942 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
943 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
944 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
945 daughter0Electron = kTRUE;
946 }
947 }
948
949 }
950 }
951 else if(tmpDaughter->GetPdgCode() == -11){
952 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
953 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
954 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
955 daughter0Positron = kTRUE;
956 }
a68437fb 957 }
958 }
d7d7e825 959 }
960 }
961 }
962 }
963
964
965 if(daughter1->GetNDaughters() >= 2){ // second gamma
966 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
967 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
968 if(tmpDaughter->GetUniqueID() == 5){
969 if(tmpDaughter->GetPdgCode() == 11){
970 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
971 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
972 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
973 daughter1Electron = kTRUE;
974 }
975 }
976
977 }
978 }
979 else if(tmpDaughter->GetPdgCode() == -11){
980 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
981 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
982 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
983 daughter1Positron = kTRUE;
984 }
985 }
986
987 }
988
989 }
990 }
991 }
992 }
993
a0b94e5c 994
d7d7e825 995 if(particle->GetPdgCode()==111){ //Pi0
996 if( iTracks >= fStack->GetNprimary()){
997 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
998 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
999 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
1000 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
1001 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
1002 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
1003 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
1004 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1005 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
1006
01b7fdcc 1007 if(gammaEtaCut && gammaRCut){
d7d7e825 1008 //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1009 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1010 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1011 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1012 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1013 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1014 }
1015 }
1016 }
1017 else{
1018 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
1019 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
1020 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
1021 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
1022 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
1023 fHistograms->FillHistogram("MC_Pi0_R", particle->R());
1024 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
1025 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1026 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
63e16c52 1027 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
1028
01b7fdcc 1029 if(gammaEtaCut && gammaRCut){
d7d7e825 1030 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1031 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1032 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
63e16c52 1033 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
1034
d7d7e825 1035 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1036 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1037 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1038 fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1039 fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1040 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1041 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1042
48682642 1043 Double_t alfa=0.;
1044 if((daughter0->Energy()+daughter1->Energy())!= 0.){
1045 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1046 }
1047 fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
1048 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
1049
d7d7e825 1050 }
1051 }
1052 }
1053 }
1054
1055 if(particle->GetPdgCode()==221){ //Eta
1056 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
1057 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
1058 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
1059 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
1060 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
1061 fHistograms->FillHistogram("MC_Eta_R", particle->R());
1062 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
1063 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1064 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1065
01b7fdcc 1066 if(gammaEtaCut && gammaRCut){
a0b94e5c 1067 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1068 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1069 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1070 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1071 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1072 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1073 fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1074 fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1075 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1076 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1077
d7d7e825 1078 }
1079
1080 }
1081
1082 }
1083
1084 // all motherparticles with 2 gammas as daughters
1085 fHistograms->FillHistogram("MC_Mother_R", particle->R());
1086 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1087 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1088 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1089 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1090 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1091 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1092 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1093 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1094 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1095 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
a0b94e5c 1096
01b7fdcc 1097 if(gammaEtaCut && gammaRCut){
a0b94e5c 1098 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1099 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1100 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1101 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
1102 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1103 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1104 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1105 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
a0b94e5c 1106
d7d7e825 1107 }
1108
1109
1110 } // end passed R and eta cut
a0b94e5c 1111
d7d7e825 1112 } // end if(particle->GetNDaughters() == 2)
1113
1114 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
a0b94e5c 1115
d7d7e825 1116} // end ProcessMCData
1117
1118
1119
1120void AliAnalysisTaskGammaConversion::FillNtuple(){
1121 //Fills the ntuple with the different values
a0b94e5c 1122
d7d7e825 1123 if(fGammaNtuple == NULL){
1124 return;
1125 }
1126 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1127 for(Int_t i=0;i<numberOfV0s;i++){
1128 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};
1129 AliESDv0 * cV0 = fV0Reader->GetV0(i);
1130 Double_t negPID=0;
1131 Double_t posPID=0;
1132 fV0Reader->GetPIDProbability(negPID,posPID);
1133 values[0]=cV0->GetOnFlyStatus();
1134 values[1]=fV0Reader->CheckForPrimaryVertex();
1135 values[2]=negPID;
1136 values[3]=posPID;
1137 values[4]=fV0Reader->GetX();
1138 values[5]=fV0Reader->GetY();
1139 values[6]=fV0Reader->GetZ();
1140 values[7]=fV0Reader->GetXYRadius();
1141 values[8]=fV0Reader->GetMotherCandidateNDF();
1142 values[9]=fV0Reader->GetMotherCandidateChi2();
1143 values[10]=fV0Reader->GetMotherCandidateEnergy();
1144 values[11]=fV0Reader->GetMotherCandidateEta();
1145 values[12]=fV0Reader->GetMotherCandidatePt();
1146 values[13]=fV0Reader->GetMotherCandidateMass();
1147 values[14]=fV0Reader->GetMotherCandidateWidth();
1148 // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1149 values[16]=fV0Reader->GetOpeningAngle();
1150 values[17]=fV0Reader->GetNegativeTrackEnergy();
1151 values[18]=fV0Reader->GetNegativeTrackPt();
1152 values[19]=fV0Reader->GetNegativeTrackEta();
1153 values[20]=fV0Reader->GetNegativeTrackPhi();
1154 values[21]=fV0Reader->GetPositiveTrackEnergy();
1155 values[22]=fV0Reader->GetPositiveTrackPt();
1156 values[23]=fV0Reader->GetPositiveTrackEta();
1157 values[24]=fV0Reader->GetPositiveTrackPhi();
1158 values[25]=fV0Reader->HasSameMCMother();
1159 if(values[25] != 0){
1160 values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1161 values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1162 }
1163 fTotalNumberOfAddedNtupleEntries++;
1164 fGammaNtuple->Fill(values);
1165 }
1166 fV0Reader->ResetV0IndexNumber();
1167
1168}
1169
1170void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1171 // Process all the V0's without applying any cuts to it
a0b94e5c 1172
d7d7e825 1173 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1174 for(Int_t i=0;i<numberOfV0s;i++){
1175 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
a0b94e5c 1176
d7d7e825 1177 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
cb90a330 1178 continue;
d7d7e825 1179 }
9640a3d1 1180
77880bd8 1181 // if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1182 if( !fV0Reader->CheckV0FinderStatus(i)){
cb90a330 1183 continue;
9640a3d1 1184 }
1185
cb90a330 1186
1187 if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
1188 !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1189 continue;
9640a3d1 1190 }
cb90a330 1191
ebcfaa7e 1192 if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1193 continue;
1194 }
1195
1196 if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1197 fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1198 continue;
1199 }
9c1cb6f7 1200 if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1201 continue;
1202 }
1203 if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1204 continue;
1205 }
1206 if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1207 continue;
1208 }
1209 if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1210 continue;
1211 }
d7d7e825 1212 if(fDoMCTruth){
a0b94e5c 1213
d7d7e825 1214 if(fV0Reader->HasSameMCMother() == kFALSE){
1215 continue;
1216 }
a0b94e5c 1217
d7d7e825 1218 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1219 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a0b94e5c 1220
d7d7e825 1221 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1222 continue;
1223 }
26923b22 1224 if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
d7d7e825 1225 continue;
1226 }
a68437fb 1227
26923b22 1228 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
a68437fb 1229 continue;
1230 }
a0b94e5c 1231
d7d7e825 1232 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a0b94e5c 1233
d7d7e825 1234 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1235 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1236 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1237 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1238 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1239 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1240 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1241 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1242 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1243 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
a0b94e5c 1244
d7d7e825 1245 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1246 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1247
1248 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1249 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1250 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1251 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1252 fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1253 fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1254 fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1255 fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1256
1257 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1258 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1259 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1260 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1261
d7d7e825 1262 //store MCTruth properties
1263 fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1264 fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1265 fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
d7d7e825 1266 }
1267 }
1268 }
1269 fV0Reader->ResetV0IndexNumber();
1270}
1271
1272void AliAnalysisTaskGammaConversion::ProcessV0s(){
1273 // see header file for documentation
037dc2db 1274
1275
d7d7e825 1276 if(fWriteNtuple == kTRUE){
1277 FillNtuple();
1278 }
1279
1280 Int_t nSurvivingV0s=0;
5ce758b0 1281 fV0Reader->ResetNGoodV0s();
d7d7e825 1282 while(fV0Reader->NextV0()){
1283 nSurvivingV0s++;
1284
e60f3265 1285
1286 TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1287
d7d7e825 1288 //-------------------------- filling v0 information -------------------------------------
1289 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1290 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1291 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1292 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
e60f3265 1293
1294 // Specific histograms for beam pipe studies
1295 if( TMath::Abs(fV0Reader->GetZ()) < fV0Reader->GetLineCutZValue() ){
1296 fHistograms->FillHistogram("ESD_Conversion_XY_BeamPipe", fV0Reader->GetX(),fV0Reader->GetY());
1297 fHistograms->FillHistogram("ESD_Conversion_RPhi_BeamPipe", vtxConv.Phi(),fV0Reader->GetXYRadius());
1298 }
1299
1300
d7d7e825 1301 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1302 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1303 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1304 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
ebcfaa7e 1305 fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1306 fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
5ce758b0 1307 if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 && fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
1308 Double_t EclsToF= (Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
1309 fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),EclsToF );
1310 fHistograms->FillHistogram("ESD_E_TPCchi2", fV0Reader->GetNegativeTrackTPCchi2()/(Double_t)fV0Reader->GetNegativeTracknTPCClusters());
1311 }
1312
1313
1314
d7d7e825 1315 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1316 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1317 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1318 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
ebcfaa7e 1319 fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1320 fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
5ce758b0 1321 if(fV0Reader->GetPositiveTracknTPCFClusters()!=0 && (Double_t)fV0Reader->GetPositiveTracknTPCClusters()!=0 ){
1322 Double_t PclsToF= (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
1323 fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), PclsToF);
1324 fHistograms->FillHistogram("ESD_P_TPCchi2", fV0Reader->GetPositiveTrackTPCchi2()/(Double_t)fV0Reader->GetPositiveTracknTPCClusters());
1325 }
1326
1327
d7d7e825 1328 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1329 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1330 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1331 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1332 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1333 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1334 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1335 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1336 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1337 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1338
1339 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1340 fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
9640a3d1 1341
1342 fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1343 fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1344 fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1345 fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1346
1347 fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1348 fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1349 fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1350 fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1351
9c1cb6f7 1352 Double_t negPID=0;
1353 Double_t posPID=0;
1354 fV0Reader->GetPIDProbability(negPID,posPID);
1355 fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1356 fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1357
1358 Double_t negPIDmupi=0;
1359 Double_t posPIDmupi=0;
1360 fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1361 fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1362 fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1363
48682642 1364 Double_t armenterosQtAlfa[2];
1365 fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(),
1366 fV0Reader-> GetPositiveKFParticle(),
1367 fV0Reader->GetMotherCandidateKFCombination(),
1368 armenterosQtAlfa);
1369
1370 fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1371
037dc2db 1372
d7d7e825 1373 // begin mapping
1374 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
9640a3d1 1375 Int_t zBin = fHistograms->GetZBin(fV0Reader->GetZ());
d7d7e825 1376 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
9c1cb6f7 1377 Double_t rFMD=30;
1378
e60f3265 1379
9640a3d1 1380
ebcfaa7e 1381 // Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
d7d7e825 1382
1383 TString nameESDMappingPhiR="";
e158cbc3 1384 nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 1385 //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1386
1387 TString nameESDMappingPhi="";
e158cbc3 1388 nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
ebcfaa7e 1389 //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1390
1391 TString nameESDMappingR="";
e158cbc3 1392 nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
ebcfaa7e 1393 //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1394
1395 TString nameESDMappingPhiInR="";
e158cbc3 1396 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 1397 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1398 fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1399
1400 TString nameESDMappingZInR="";
1401 nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1402 fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1403
1404 TString nameESDMappingPhiInZ="";
1405 nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1406 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1407 fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1408
9c1cb6f7 1409 if(fV0Reader->GetXYRadius()<rFMD){
1410 TString nameESDMappingFMDPhiInZ="";
1411 nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1412 fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1413 }
1414
1415
9640a3d1 1416 TString nameESDMappingRInZ="";
1417 nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1418 fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1419
1420 if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1421 TString nameESDMappingMidPtPhiInR="";
1422 nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1423 fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1424
1425 TString nameESDMappingMidPtZInR="";
1426 nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1427 fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1428
1429 TString nameESDMappingMidPtPhiInZ="";
1430 nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1431 fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
9c1cb6f7 1432 if(fV0Reader->GetXYRadius()<rFMD){
1433 TString nameESDMappingMidPtFMDPhiInZ="";
1434 nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1435 fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1436 }
1437
9640a3d1 1438
1439 TString nameESDMappingMidPtRInZ="";
1440 nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1441 fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1442
1443 }
1444
1445
d7d7e825 1446 // end mapping
1447
6c84d371 1448 new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
037dc2db 1449 fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
6c84d371 1450 // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
d7d7e825 1451 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1452 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
a0b94e5c 1453
d7d7e825 1454 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1455 if(fDoMCTruth){
1456
a68437fb 1457 if(fV0Reader->HasSameMCMother() == kFALSE){
1458 continue;
1459 }
d7d7e825 1460 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
a68437fb 1461 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1462
1463 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1464 continue;
1465 }
1466
1467 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1468 continue;
1469 }
bd6d9fa3 1470 if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
1471 (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays
037dc2db 1472 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1473 fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1474 }
1475 }
a68437fb 1476
1477 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1478 continue;
1479 }
a0b94e5c 1480
d7d7e825 1481 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a68437fb 1482
26923b22 1483 if(fDoCF){
1484 Double_t containerInput[3];
1485 containerInput[0] = fV0Reader->GetMotherCandidatePt();
1486 containerInput[1] = fV0Reader->GetMotherCandidateEta();
1487 containerInput[2] = fV0Reader->GetMotherCandidateMass();
1488 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1489 }
1490
d7d7e825 1491 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1492 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1493 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1494 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1495 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1496 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1497 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1498 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1499 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1500 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1501 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1502 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1503 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1504 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
a0b94e5c 1505
d7d7e825 1506 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1507 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
a0b94e5c 1508
d7d7e825 1509
1510 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1511 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1512 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1513 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1514
1515 fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1516 fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1517 fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1518 fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
96ade8ef 1519 if (fV0Reader->GetMotherCandidateP() != 0) {
1520 fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1521 fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1522 } else { cout << "Error::fV0Reader->GetNegativeTrackP() == 0 !!!" << endl; }
1523
9640a3d1 1524 fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1525 fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
70ef88b5 1526 fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1527
9640a3d1 1528
a0b94e5c 1529
d7d7e825 1530 //store MCTruth properties
1531 fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1532 fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1533 fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
48682642 1534
d7d7e825 1535 //resolution
1536 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
1537 Double_t esdpt = fV0Reader->GetMotherCandidatePt();
9c1cb6f7 1538 Double_t resdPt = 0.;
4a6157dc 1539 if(mcpt > 0){
9c1cb6f7 1540 resdPt = ((esdpt - mcpt)/mcpt)*100.;
d7d7e825 1541 }
4a6157dc 1542 else if(mcpt < 0){
1543 cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
1544 }
d7d7e825 1545
ca6d4600 1546 fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
d7d7e825 1547 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1548 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
ca6d4600 1549 fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
d7d7e825 1550
9c1cb6f7 1551 Double_t resdZ = 0.;
d7d7e825 1552 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
9c1cb6f7 1553 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
d7d7e825 1554 }
9c1cb6f7 1555 Double_t resdZAbs = 0.;
037dc2db 1556 resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1557
1558 fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
d7d7e825 1559 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1560 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1561 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
ca6d4600 1562
1563 // new for dPt_Pt-histograms for Electron and Positron
1564 Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
9c1cb6f7 1565 Double_t resEdPt = 0.;
ca6d4600 1566 if (mcEpt > 0){
9c1cb6f7 1567 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
ca6d4600 1568 }
1569 UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus();
1570 // AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
1571 UInt_t kTRDoutN = (statusN & AliESDtrack::kTRDout);
1572
1573 Int_t ITSclsE= fV0Reader->GetNegativeTracknITSClusters();
1574 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1575 switch(ITSclsE){
1576 case 0: // 0 ITS clusters
1577 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
1578 break;
1579 case 1: // 1 ITS cluster
1580 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
1581 break;
1582 case 2: // 2 ITS clusters
1583 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
1584 break;
1585 case 3: // 3 ITS clusters
1586 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
1587 break;
1588 case 4: // 4 ITS clusters
1589 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
1590 break;
1591 case 5: // 5 ITS clusters
1592 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
1593 break;
1594 case 6: // 6 ITS clusters
1595 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
1596 break;
1597 }
1598 //Filling histograms with respect to Electron resolution
1599 fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
1600 fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
1601 if(kTRDoutN){
1602 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1603 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1604 fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1605 fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1606 fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
1607 }
1608
1609 Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
1610 Double_t resPdPt = 0;
1611 if (mcPpt > 0){
9c1cb6f7 1612 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
ca6d4600 1613 }
1614
1615 UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus();
1616// AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
1617 UInt_t kTRDoutP = (statusP & AliESDtrack::kTRDout);
1618
1619 Int_t ITSclsP = fV0Reader->GetPositiveTracknITSClusters();
1620 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1621 switch(ITSclsP){
1622 case 0: // 0 ITS clusters
1623 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
1624 break;
1625 case 1: // 1 ITS cluster
1626 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
1627 break;
1628 case 2: // 2 ITS clusters
1629 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
1630 break;
1631 case 3: // 3 ITS clusters
1632 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
1633 break;
1634 case 4: // 4 ITS clusters
1635 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
1636 break;
1637 case 5: // 5 ITS clusters
1638 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
1639 break;
1640 case 6: // 6 ITS clusters
1641 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
1642 break;
1643 }
1644 //Filling histograms with respect to Positron resolution
1645 fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
1646 fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
1647 if(kTRDoutP){
1648 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1649 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1650 fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1651 fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1652 fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
1653 }
1654
1655
9c1cb6f7 1656 Double_t resdR = 0.;
d7d7e825 1657 if(fV0Reader->GetNegativeMCParticle()->R() != 0){
9c1cb6f7 1658 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
d7d7e825 1659 }
9c1cb6f7 1660 Double_t resdRAbs = 0.;
037dc2db 1661 resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1662
1663 fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
d7d7e825 1664 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1665 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1666 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
ca6d4600 1667 fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1668
9c1cb6f7 1669 Double_t resdPhiAbs=0.;
1670 resdPhiAbs=0.;
037dc2db 1671 resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1672 fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1673
d7d7e825 1674 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
d7d7e825 1675 }//if(fDoMCTruth)
1676 }//while(fV0Reader->NextV0)
1677 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1678 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
b5832f95 1679 fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
cb90a330 1680
1681 fV0Reader->ResetV0IndexNumber();
d7d7e825 1682}
1683
1684void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1685 // Fill AOD with reconstructed Gamma
a0b94e5c 1686
6c84d371 1687 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1688 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
d7d7e825 1689 //Create AOD particle object from AliKFParticle
d7d7e825 1690 //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1691 //but this means that I have to work a little bit more in my side.
1692 //AODPWG4Particle objects are simpler and lighter, I think
332f1f44 1693 /*
1694 AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
d7d7e825 1695 AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
332f1f44 1696 //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1697 gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
a0b94e5c 1698 gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
332f1f44 1699 gamma.SetPdg(AliPID::kEleCon); //photon id
a0b94e5c 1700 gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
332f1f44 1701 gamma.SetChi2(gammakf->Chi2());
a0b94e5c 1702 Int_t i = fAODBranch->GetEntriesFast();
1703 new((*fAODBranch)[i]) AliAODPWG4Particle(gamma);
d7d7e825 1704 */
332f1f44 1705
6c84d371 1706 AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
d7d7e825 1707 AliGammaConversionAODObject aodObject;
1708 aodObject.SetPx(gammakf->GetPx());
1709 aodObject.SetPy(gammakf->GetPy());
1710 aodObject.SetPz(gammakf->GetPz());
1711 aodObject.SetLabel1(fElectronv1[gammaIndex]);
1712 aodObject.SetLabel2(fElectronv2[gammaIndex]);
332f1f44 1713 aodObject.SetChi2(gammakf->Chi2());
04bf4381 1714 aodObject.SetE(gammakf->E());
1715 Int_t i = fAODGamma->GetEntriesFast();
1716 new((*fAODGamma)[i]) AliGammaConversionAODObject(aodObject);
a0b94e5c 1717 }
1718
d7d7e825 1719}
1720
6272370b 1721void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
1722 // omega meson analysis pi0+gamma decay
1723 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
48682642 1724 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
6272370b 1725 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
48682642 1726
6272370b 1727 AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1728 if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
1729 continue;
1730 }
1731
48682642 1732 AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
6272370b 1733 Double_t massOmegaCandidate = 0.;
1734 Double_t widthOmegaCandidate = 0.;
1735
48682642 1736 omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
6272370b 1737
e9aea39f 1738 if ( massOmegaCandidate > 733 && massOmegaCandidate < 833 ) {
c59360eb 1739 AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
1740 }
1741
48682642 1742 fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
6272370b 1743 fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
48682642 1744
1745 //delete omegaCandidate;
1746
1747 }// end of omega reconstruction in pi0+gamma channel
1748
1749 if(fDoJet == kTRUE){
1750 AliKFParticle* negPiKF=NULL;
1751 AliKFParticle* posPiKF=NULL;
1752
1753 // look at the pi+pi+pi0 channel
1754 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1755 AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1756 if (posTrack->GetSign()<0) continue;
c0d9051e 1757 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion))>2.) continue;
48682642 1758 if (posPiKF) delete posPiKF; posPiKF=NULL;
1759 posPiKF = new AliKFParticle( *(posTrack) ,211);
1760
1761 for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
1762 AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
1763 if( negTrack->GetSign()>0) continue;
c0d9051e 1764 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion))>2.) continue;
48682642 1765 if (negPiKF) delete negPiKF; negPiKF=NULL;
1766 negPiKF = new AliKFParticle( *(negTrack) ,-211);
1767 AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
1768 Double_t massOmegaCandidatePipPinPi0 = 0.;
1769 Double_t widthOmegaCandidatePipPinPi0 = 0.;
1770
1771 omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
ef2e2748 1772
1773 if ( massOmegaCandidatePipPinPi0 > 733 && massOmegaCandidatePipPinPi0 < 833 ) {
1774 AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
1775 }
48682642 1776
1777 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
1778 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
1779
1780 // delete omegaCandidatePipPinPi0;
1781 }
1782 }
1783 } // checking ig gammajet because in that case the chargedparticle list is created
1784
1785
6272370b 1786
6272370b 1787 }
48682642 1788
1789 if(fCalculateBackground){
5ce758b0 1790
1791 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
1792
1793 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
1794 Int_t mbin = 0;
1795 if(fUseTrackMultiplicityForBG == kTRUE){
1796 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1797 }
1798 else{
1799 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
1800 }
1801
1802 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
1803
48682642 1804 // Background calculation for the omega
1805 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 1806 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1807
1808 if(fMoveParticleAccordingToVertex == kTRUE){
1809 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
1810 }
48682642 1811 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1812 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 1813
1814 if(fMoveParticleAccordingToVertex == kTRUE){
1815 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1816 }
1817
48682642 1818 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
1819 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
1820 AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
1821 Double_t massOmegaBckCandidate = 0.;
1822 Double_t widthOmegaBckCandidate = 0.;
1823
1824 omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
ef2e2748 1825
1826
48682642 1827 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
1828 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
1829
1830 delete omegaBckCandidate;
1831 }
1832 }
1833 }
1834 } // end of checking if background calculation is available
6272370b 1835}
d7d7e825 1836
04bf4381 1837
1838void AliAnalysisTaskGammaConversion::AddOmegaToAOD(AliKFParticle * omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
1839 //See header file for documentation
1840 AliGammaConversionAODObject omega;
1841 omega.SetPx(omegakf->Px());
1842 omega.SetPy(omegakf->Py());
1843 omega.SetPz(omegakf->Pz());
1844 omega.SetChi2(omegakf->GetChi2());
1845 omega.SetE(omegakf->E());
1846 omega.SetIMass(mass);
1847 omega.SetLabel1(omegaDaughter);
1848 //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
1849 omega.SetLabel2(gammaDaughter);
1850 new((*fAODOmega)[fAODOmega->GetEntriesFast()]) AliGammaConversionAODObject(omega);
1851}
1852
1853
1854
d7d7e825 1855void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1856 // see header file for documentation
1857
6c84d371 1858 // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1859 // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
037dc2db 1860
5ce758b0 1861 fESDEvent = fV0Reader->GetESDEvent();
1862
037dc2db 1863 if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
1864 cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
1865 }
1866
6c84d371 1867 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1868 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
d7d7e825 1869
6c84d371 1870 // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1871 // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
a0b94e5c 1872
6c84d371 1873 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1874 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
a0b94e5c 1875
d7d7e825 1876 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1877 continue;
1878 }
1879 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1880 continue;
1881 }
a0b94e5c 1882
d7d7e825 1883 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1884
1885 Double_t massTwoGammaCandidate = 0.;
1886 Double_t widthTwoGammaCandidate = 0.;
1887 Double_t chi2TwoGammaCandidate =10000.;
1888 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
d707e3cf 1889 // if(twoGammaCandidate->GetNDF()>0){
1890 // chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1891 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2();
d7d7e825 1892
70ef88b5 1893 fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
2e2da371 1894 if((chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
d7d7e825 1895
1896 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1897 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1898
1899 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
1900 Double_t rapidity;
96ade8ef 1901 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() <= 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() <= 0){
1902 cout << "Error: |Pz| > E !!!! " << endl;
d7d7e825 1903 rapidity=0;
1904 }
1905 else{
1906 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1907 }
48682642 1908 Double_t alfa=0.0;
1909 if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
1910 alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
1911 /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
1912 }
d7d7e825 1913
1e7846f4 1914 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
1915 delete twoGammaCandidate;
1916 continue; // minimum opening angle to avoid using ghosttracks
1917 }
1918
d7d7e825 1919 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1920 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1921 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1922 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1923 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
1924 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1925 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
48682642 1926 fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
d7d7e825 1927 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
1928 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1929 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1930 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
67381a40 1931 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
9c1cb6f7 1932 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1933 }
d7d7e825 1934 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
9640a3d1 1935
9c1cb6f7 1936 if(fCalculateBackground){
5ce758b0 1937 /* Kenneth, just for testing*/
1938 AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
1939
1940 Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
1941 Int_t mbin=0;
1942 Int_t multKAA=0;
1943 if(fUseTrackMultiplicityForBG == kTRUE){
1944 multKAA=fV0Reader->CountESDTracks();
1945 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1946 }
1947 else{// means we use #v0s for multiplicity
1948 multKAA=fV0Reader->GetNGoodV0s();
1949 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
1950 }
1951 // cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
1952 // cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
1953 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
1954 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1955 /* end Kenneth, just for testing*/
1956 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1957 }
1958 /* if(fCalculateBackground){
9c1cb6f7 1959 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
1960 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1961 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
5ce758b0 1962 }*/
ebcfaa7e 1963 // if(fDoNeutralMesonV0MCCheck){
1964 if(fDoMCTruth){
037dc2db 1965 //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
1966 Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
1967 if(indexKF1<fV0Reader->GetNumberOfV0s()){
1968 fV0Reader->GetV0(indexKF1);//updates to the correct v0
1969 Double_t eta1 = fV0Reader->GetMotherCandidateEta();
1970 Bool_t isRealPi0=kFALSE;
10e3319b 1971 Bool_t isRealEta=kFALSE;
037dc2db 1972 Int_t gamma1MotherLabel=-1;
1973 if(fV0Reader->HasSameMCMother() == kTRUE){
1974 //cout<<"This v0 is a real v0!!!!"<<endl;
1975 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1976 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1977 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1978 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1979 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1980 gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1981 }
1982 }
1983 }
1984 }
1985 Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
1986 if(indexKF1 == indexKF2){
1987 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
1988 }
1989 if(indexKF2<fV0Reader->GetNumberOfV0s()){
1990 fV0Reader->GetV0(indexKF2);
1991 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
1992 Int_t gamma2MotherLabel=-1;
1993 if(fV0Reader->HasSameMCMother() == kTRUE){
1994 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1995 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1996 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1997 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1998 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1999 gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2000 }
2001 }
2002 }
2003 }
2004 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2005 if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2006 isRealPi0=kTRUE;
2007 }
10e3319b 2008 if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2009 isRealEta=kTRUE;
2010 }
2011
037dc2db 2012 }
2013
2014 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
ebcfaa7e 2015 // fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2016 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
10e3319b 2017 if(isRealPi0 || isRealEta){
037dc2db 2018 fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2019 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2020 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
ebcfaa7e 2021 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2022 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2023 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2024 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2025 }
037dc2db 2026 }
bd6d9fa3 2027 if(!isRealPi0 && !isRealEta){
2028 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2029 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2030 }else{
2031 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2032 }
2033 }
2034
037dc2db 2035 }
2036 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
ebcfaa7e 2037 // fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2038 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
10e3319b 2039 if(isRealPi0 || isRealEta){
037dc2db 2040 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2041 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2042 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
ebcfaa7e 2043 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2044 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2045 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2046 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2047 }
037dc2db 2048 }
bd6d9fa3 2049 if(!isRealPi0 && !isRealEta){
2050 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2051 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2052 }else{
2053 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2054 }
2055 }
037dc2db 2056 }
2057 else{
ebcfaa7e 2058 // fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2059 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
10e3319b 2060 if(isRealPi0 || isRealEta){
037dc2db 2061 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2062 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2063 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
ebcfaa7e 2064 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2065 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2066 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2067 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2068 }
bd6d9fa3 2069 if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
2070 fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
2071 }
2072 }
2073 if(!isRealPi0 && !isRealEta){
2074 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2075 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2076 }else{
2077 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2078 }
037dc2db 2079 }
2080 }
2081 }
2082 }
2083 }
9640a3d1 2084 if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2085 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2086 fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2087 }
ebcfaa7e 2088
2089 if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2090 fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2091 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2092 }
2093 else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2094 fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2095 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2096 }
2097 else{
2098 fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2099 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2100 }
6272370b 2101 Double_t lowMassPi0=0.1;
2102 Double_t highMassPi0=0.15;
2103 if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
2104 new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()]) AliKFParticle(*twoGammaCandidate);
2105 fGammav1.push_back(firstGammaIndex);
2106 fGammav2.push_back(secondGammaIndex);
c59360eb 2107 AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
6272370b 2108 }
2109
2e2da371 2110 }
d707e3cf 2111 //}
2112 delete twoGammaCandidate;
d7d7e825 2113 }
2114 }
2115}
2116
04bf4381 2117void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
2118 //See header file for documentation
2119 AliGammaConversionAODObject pion;
2120 pion.SetPx(pionkf->Px());
2121 pion.SetPy(pionkf->Py());
2122 pion.SetPz(pionkf->Pz());
2123 pion.SetChi2(pionkf->GetChi2());
2124 pion.SetE(pionkf->E());
2125 pion.SetIMass(mass);
2126 pion.SetLabel1(daughter1);
2127 //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
2128 pion.SetLabel2(daughter2);
2129 new((*fAODPi0)[fAODPi0->GetEntriesFast()]) AliGammaConversionAODObject(pion);
2130
2131}
2132
2133
2134
6272370b 2135void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
48682642 2136/*
6272370b 2137 // see header file for documentation
2138 // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2139
2140
2141
2142 Double_t vtx[3];
2143 vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2144 vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2145 vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2146
2147
2148 // Loop over all CaloClusters and consider only the PHOS ones:
2149 AliESDCaloCluster *clu;
2150 TLorentzVector pPHOS;
2151 TLorentzVector gammaPHOS;
2152 TLorentzVector gammaGammaConv;
2153 TLorentzVector pi0GammaConvPHOS;
2154 TLorentzVector gammaGammaConvBck;
2155 TLorentzVector pi0GammaConvPHOSBck;
2156
2157
2158 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2159 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2160 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2161 clu ->GetMomentum(pPHOS ,vtx);
2162 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2163 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2164 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2165 gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2166 pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2167 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2168 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2169
2170 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2171 TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2172 Double_t opanConvPHOS= v3D0.Angle(v3D1);
2173 if ( opanConvPHOS < 0.35){
2174 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2175 }else{
2176 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2177 }
2178
2179 }
2180
2181 // Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2182 }
2183 //==== End of the PHOS cluster selection ============
2184 TLorentzVector pEMCAL;
2185 TLorentzVector gammaEMCAL;
2186 TLorentzVector pi0GammaConvEMCAL;
2187 TLorentzVector pi0GammaConvEMCALBck;
2188
2189 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2190 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2191 if ( !clu->IsEMCAL() || clu->E()<0.1 ) continue;
2192 if (clu->GetNCells() <= 1) continue;
2193 if ( clu->GetTOF()*1e9 < 550 || clu->GetTOF()*1e9 > 750) continue;
2194
2195 clu ->GetMomentum(pEMCAL ,vtx);
2196 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2197 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2198 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2199 twoGammaDecayCandidateDaughter0->Py(),
2200 twoGammaDecayCandidateDaughter0->Pz(),0.);
2201 gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2202 pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2203 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2204 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2205 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2206 twoGammaDecayCandidateDaughter0->Py(),
2207 twoGammaDecayCandidateDaughter0->Pz());
2208 TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2209
2210
2211 Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2212 if ( opanConvEMCAL < 0.35){
2213 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2214 }else{
2215 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2216 }
2217
2218 }
48682642 2219 if(fCalculateBackground){
2220 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2221 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2222 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2223 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2224 gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2225 previousGoodV0.Py(),
2226 previousGoodV0.Pz(),0.);
2227 pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2228 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2229 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2230 pi0GammaConvEMCALBck.Pt());
2231 }
6272370b 2232 }
48682642 2233
2234 // Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2235 } // end of checking if background photons are available
2236 }
6272370b 2237 //==== End of the PHOS cluster selection ============
48682642 2238*/
6272370b 2239}
2240
5ce758b0 2241void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle *particle,AliGammaConversionBGHandler::GammaConversionVertex *vertex){
2242 //see header file for documentation
2243
2244 Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2245 Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2246 Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
2247
2248 // cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
2249 particle->X() = particle->GetX() - dx;
2250 particle->Y() = particle->GetY() - dy;
2251 particle->Z() = particle->GetZ() - dz;
2252}
2253
111d75df 2254void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
2255
2256 Double_t c = cos(angle);
2257 Double_t s = sin(angle);
2258
2259 Double_t A[7][ 7];
2260 for( Int_t i=0; i<7; i++ ){
2261 for( Int_t j=0; j<7; j++){
2262 A[i][j] = 0;
2263 }
2264 }
2265 for( int i=0; i<7; i++ ){
2266 A[i][i] = 1;
2267 }
2268 A[0][0] = c; A[0][1] = s;
2269 A[1][0] = -s; A[1][1] = c;
2270 A[3][3] = c; A[3][4] = s;
2271 A[4][3] = -s; A[4][4] = c;
2272
2273 Double_t AC[7][7];
2274 Double_t Ap[7];
2275
2276 for( Int_t i=0; i<7; i++ ){
2277 Ap[i] = 0;
2278 for( Int_t k=0; k<7; k++){
2279 Ap[i]+=A[i][k] * kfParticle->GetParameter(k);
2280 }
2281 }
2282
2283 for( Int_t i=0; i<7; i++){
2284 kfParticle->Parameter(i) = Ap[i];
2285 }
2286
2287 for( Int_t i=0; i<7; i++ ){
2288 for( Int_t j=0; j<7; j++ ){
2289 AC[i][j] = 0;
2290 for( Int_t k=0; k<7; k++ ){
2291 AC[i][j]+= A[i][k] * kfParticle->GetCovariance(k,j);
2292 }
2293 }
2294 }
2295
2296 for( Int_t i=0; i<7; i++ ){
2297 for( Int_t j=0; j<=i; j++ ){
2298 Double_t xx = 0;
2299 for( Int_t k=0; k<7; k++){
2300 xx+= AC[i][k]*A[j][k];
2301 }
2302 kfParticle->Covariance(i,j) = xx;
2303 }
2304 }
2305}
2306
2307
d7d7e825 2308void AliAnalysisTaskGammaConversion::CalculateBackground(){
2309 // see header file for documentation
5e55d806 2310
2311
2312 TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2313
5ce758b0 2314 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2315
2316 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2317 Int_t mbin = 0;
2318 if(fUseTrackMultiplicityForBG == kTRUE){
2319 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2320 }
2321 else{
2322 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2323 }
2324
111d75df 2325 if(fDoRotation == kTRUE){
2326 TRandom *random = new TRandom();
2327 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2328 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2329 for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
2330 for(Int_t nRandom=0;nRandom<nRandomEventsForBG;nRandom++){
2331
2332 AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
2333
2334 Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
2335
2336 Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
5ce758b0 2337
111d75df 2338 RotateKFParticle(&currentEventGoodV02,rotationValue);
5ce758b0 2339
111d75df 2340 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
5e55d806 2341
5ce758b0 2342 Double_t massBG =0.;
2343 Double_t widthBG = 0.;
2344 Double_t chi2BG =10000.;
2345 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2346
2347 if(fCheckBGProbability == kTRUE){
2348 if(massBG>0.1 && massBG<1.4){
2349 if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2350 delete backgroundCandidate;
2351 continue;
2352 }
2353 }
5ce758b0 2354 }
111d75df 2355
2356 // if(backgroundCandidate->GetNDF()>0){
2357 chi2BG = backgroundCandidate->GetChi2();
2358 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2359
2360 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2361 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2362
2363 Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2364
2365 Double_t rapidity;
2366 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2367 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2368
2369
2370 Double_t alfa=0.0;
2371 if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2372 alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2373 /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2374 }
2375
2376
2377 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2378 delete backgroundCandidate;
2379 continue; // minimum opening angle to avoid using ghosttracks
2380 }
2381
2382 // original
2383 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2384 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2385 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2386 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2387 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2388 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2389 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2390 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2391 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2392 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2393 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2394 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2395
2396 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2397 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2398 }
2399 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2400 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2401 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2402 }
2403
2404 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2405 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2406 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2407 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2408 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2409 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2410 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2411 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2412 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2413 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2414 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2415 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2416
2417 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2418 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2419 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2420 }
2e2da371 2421 }
111d75df 2422 //}
5ce758b0 2423 delete backgroundCandidate;
2424 }
2425 }
2426 }
111d75df 2427 delete random;
5ce758b0 2428 }
111d75df 2429 else{ // means no rotation
2430 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2431
2432 if(fUseTrackMultiplicityForBG){
2433 // cout<<"Using charged track multiplicity for background calculation"<<endl;
2434 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 2435
111d75df 2436 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2437
5ce758b0 2438 if(fMoveParticleAccordingToVertex == kTRUE){
2439 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2440 }
2441
2442 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2443 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2444 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2445 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
111d75df 2446 AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 2447
111d75df 2448 //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2449 //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2450
2451 //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2452 if(fMoveParticleAccordingToVertex == kTRUE){
2453 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2454 }
111d75df 2455 //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2456
2457 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
111d75df 2458
5ce758b0 2459 Double_t massBG =0.;
2460 Double_t widthBG = 0.;
2461 Double_t chi2BG =10000.;
2462 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2463 // if(backgroundCandidate->GetNDF()>0){
2464 // chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2465 chi2BG = backgroundCandidate->GetChi2();
2466 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2467
2468 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2469 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2470
2471 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2472
2473 Double_t rapidity;
2474
2475 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2476 cout << "Error: |Pz| > E !!!! " << endl;
2477 rapidity=0;
2478 } else {
2479 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2480 }
2481
2482 Double_t alfa=0.0;
2483 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2484 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2485 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2486 }
2487
2488
2489 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2490 delete backgroundCandidate;
2491 continue; // minimum opening angle to avoid using ghosttracks
2492 }
2493
2494 // original
2495 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2496 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2497 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2498 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2499 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2500 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2501 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2502 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2503 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2504 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2505 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2506 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2507
2508 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2509 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2510 }
2511 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2512 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2513 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2514 }
2515
2516
2517 // test
2518 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2519 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2520 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2521 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2522 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2523 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2524 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2525 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2526 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2527 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2528 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2529 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2530
2531 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2532 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2533 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2534 }
2535 // }
2536 }
2537 delete backgroundCandidate;
2538 }
2539 }
2540 }
2541 }
2542 else{ // means using #V0s for multiplicity
2543
2544 // cout<<"Using the v0 multiplicity to calculate background"<<endl;
2545
2546 fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2547 fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2548
2549 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2550 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2551 if(previousEventV0s){
2552
2553 if(fMoveParticleAccordingToVertex == kTRUE){
2554 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2555 }
2556
2557 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2558 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2559 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2560 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2561
2562 if(fMoveParticleAccordingToVertex == kTRUE){
2563 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
5ce758b0 2564 }
111d75df 2565
2566 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2567 Double_t massBG =0.;
2568 Double_t widthBG = 0.;
2569 Double_t chi2BG =10000.;
2570 backgroundCandidate->GetMass(massBG,widthBG);
2571 /* if(backgroundCandidate->GetNDF()>0){
2572 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2573 {//remember to remove
2574 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2575 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2576
2577 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2578 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2579 }
2580 */
2e2da371 2581 chi2BG = backgroundCandidate->GetChi2();
2582 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
5ce758b0 2583 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2584 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2585
2586 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2587
2588 Double_t rapidity;
2589 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2590 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2591
2592
2593 Double_t alfa=0.0;
2594 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2595 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2596 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2597 }
2598
2599
2600 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2601 delete backgroundCandidate;
2602 continue; // minimum opening angle to avoid using ghosttracks
2603 }
2604
2605 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2606 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2607 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2608 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2609 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2610 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2611 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2612 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2613 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2614 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2615 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2616 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2617
2618 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2619 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2620 }
2621 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2622 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2623 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2624 }
2625
2626 if(massBG>0.5 && massBG<0.6){
2627 fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
2628 }
2629 if(massBG>0.3 && massBG<0.4){
2630 fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
2631 }
037dc2db 2632
5ce758b0 2633 // test
2634 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2635 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2636 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2637 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2638 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2639 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2640 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2641 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2642 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2643 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2644 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2645 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2646
2647 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2648 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2649 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2650 }
2e2da371 2651 // }
5ce758b0 2652 }
111d75df 2653 delete backgroundCandidate;
2654 }
5ce758b0 2655 }
2656 }
d7d7e825 2657 }
111d75df 2658 } // end else (means use #v0s as multiplicity)
2659 } // end no rotation
d7d7e825 2660}
2661
2662
d7d7e825 2663void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
2664 //ProcessGammasForGammaJetAnalysis
a0b94e5c 2665
d7d7e825 2666 Double_t distIsoMin;
a0b94e5c 2667
d7d7e825 2668 CreateListOfChargedParticles();
a0b94e5c 2669
2670
6c84d371 2671 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
2672 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2673 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 2674 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 2675
01b7fdcc 2676 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
d7d7e825 2677 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
a0b94e5c 2678
d7d7e825 2679 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
01b7fdcc 2680 CalculateJetCone(gammaIndex);
d7d7e825 2681 }
2682 }
2683 }
2684}
2685
6272370b 2686//____________________________________________________________________
2687Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
2688{
2689//
2690// check whether particle has good DCAr(Pt) impact
2691// parameter. Only for TPC+ITS tracks (7*sigma cut)
2692// Origin: Andrea Dainese
2693//
2694
2695Float_t d0z0[2],covd0z0[3];
2696track->GetImpactParameters(d0z0,covd0z0);
2697Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
2698Float_t d0max = 7.*sigma;
2699if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
2700
2701return kFALSE;
2702}
2703
2704
d7d7e825 2705void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
2706 // CreateListOfChargedParticles
a0b94e5c 2707
d7d7e825 2708 fESDEvent = fV0Reader->GetESDEvent();
037dc2db 2709 Int_t numberOfESDTracks=0;
d7d7e825 2710 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2711 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 2712
d7d7e825 2713 if(!curTrack){
2714 continue;
2715 }
d707e3cf 2716 // Not needed if Standard function used.
2717// if(!IsGoodImpPar(curTrack)){
2718// continue;
2719// }
a0b94e5c 2720
d7d7e825 2721 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
6c84d371 2722 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
2723 // fChargedParticles.push_back(curTrack);
d7d7e825 2724 fChargedParticlesId.push_back(iTracks);
037dc2db 2725 numberOfESDTracks++;
d7d7e825 2726 }
2727 }
037dc2db 2728 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
e40fd7e2 2729
2730 if (fV0Reader->GetNumberOfContributorsVtx()>=1){
2731 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
2732 }
d7d7e825 2733}
9c1cb6f7 2734void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
2735
2736 Double_t massE=0.00051099892;
2737 TLorentzVector curElecPos;
2738 TLorentzVector curElecNeg;
2739 TLorentzVector curGamma;
2740
2741 TLorentzVector curGammaAt;
2742 TLorentzVector curElecPosAt;
2743 TLorentzVector curElecNegAt;
2744 AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
2745 AliKFVertex primVtxImprovedGamma = primVtxGamma;
2746
2747 const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
2748
2749 Double_t xPrimaryVertex=vtxT3D->GetXv();
2750 Double_t yPrimaryVertex=vtxT3D->GetYv();
2751 Double_t zPrimaryVertex=vtxT3D->GetZv();
70ef88b5 2752 // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
9c1cb6f7 2753
2754 Float_t nsigmaTPCtrackPos;
2755 Float_t nsigmaTPCtrackNeg;
2756 Float_t nsigmaTPCtrackPosToPion;
2757 Float_t nsigmaTPCtrackNegToPion;
2758 AliKFParticle* negKF=NULL;
2759 AliKFParticle* posKF=NULL;
2760
2761 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2762 AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
2763 if(!posTrack){
2764 continue;
2765 }
2766 if (posKF) delete posKF; posKF=NULL;
2767 if(posTrack->GetSign()<0) continue;
2768 if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2769 if(posTrack->GetKinkIndex(0)>0 ) continue;
2770 if(posTrack->GetNcls(1)<50)continue;
2771 Double_t momPos[3];
2772 // posTrack->GetConstrainedPxPyPz(momPos);
2773 posTrack->GetPxPyPz(momPos);
2774 AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
2775 curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
2776 if(TMath::Abs(curElecPos.Eta())<0.9) continue;
2777 posKF = new AliKFParticle( *(posTrack),-11);
2778
2779 nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
2780 nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
2781
2782 if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
2783 continue;
2784 }
2785
2786 if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
2787 continue;
2788 }
2789
2790
2791
2792 for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
2793 AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
2794 if(!negTrack){
2795 continue;
2796 }
2797 if (negKF) delete negKF; negKF=NULL;
2798 if(negTrack->GetSign()>0) continue;
2799 if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2800 if(negTrack->GetKinkIndex(0)>0 ) continue;
2801 if(negTrack->GetNcls(1)<50)continue;
2802 Double_t momNeg[3];
2803 // negTrack->GetConstrainedPxPyPz(momNeg);
2804 negTrack->GetPxPyPz(momNeg);
2805
2806 nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
2807 nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
2808 if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
2809 continue;
2810 }
2811 if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
2812 continue;
2813 }
2814 AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
2815 curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
2816 if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
2817 negKF = new AliKFParticle( *(negTrack) ,11);
2818
2819 Double_t b=fESDEvent->GetMagneticField();
2820 Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
2821 AliExternalTrackParam nt(*ntrk), pt(*ptrk);
2822 nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
2823
2824
2825 //--- Like in ITSV0Finder
2826 AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
2827 Double_t xxP,yyP,alphaP;
2828 Double_t rP[3];
2829
2830 // if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
2831 if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
2832 xxP=rP[0];
2833 yyP=rP[1];
2834 alphaP = TMath::ATan2(yyP,xxP);
2835
2836
2837 ptAt0.Propagate(alphaP,0,b);
2838 Float_t ptfacP = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
2839
2840 // Double_t distP = ptAt0.GetY();
2841 Double_t normP = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
2842 Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
2843 Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
2844 Double_t normdistP = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
2845
2846
2847 Double_t xxN,yyN,alphaN;
2848 Double_t rN[3];
2849 // if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
2850 if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
2851 xxN=rN[0];
2852 yyN=rN[1];
2853
2854 alphaN = TMath::ATan2(yyN,xxN);
2855
2856 ntAt0.Propagate(alphaN,0,b);
2857
2858 Float_t ptfacN = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
2859 // Double_t distN = ntAt0.GetY();
2860 Double_t normN = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
2861 Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
2862 Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
2863 Double_t normdistN = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
2864
2865 //-----------------------------
2866
2867 Double_t momNegAt[3];
2868 nt.GetPxPyPz(momNegAt);
2869 curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
2870
2871 Double_t momPosAt[3];
2872 pt.GetPxPyPz(momPosAt);
2873 curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
2874 if(dca>1){
2875 continue;
2876 }
2877
2878 // Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2879 // Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2880 AliESDv0 vertex(nt,jTracks,pt,iTracks);
2881
2882
2883 Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
2884
9c1cb6f7 2885
70ef88b5 2886
9c1cb6f7 2887 // cout<< "v0Rr::"<< v0Rr<<endl;
2888 // if (pvertex.GetRr()<0.5){
2889 // continue;
2890 //}
2891 // cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
2892 if(cpa<0.9)continue;
2893 // if (vertex.GetChi2V0() > 30) continue;
2894 // cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
2895 if ((xn+xp) < 0.4) continue;
2896 if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
2897 if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
2898
2899 //cout<<"pass"<<endl;
2900
2901 AliKFParticle v0GammaC;
2902 v0GammaC+=(*negKF);
2903 v0GammaC+=(*posKF);
2904 v0GammaC.SetMassConstraint(0,0.001);
2905 primVtxImprovedGamma+=v0GammaC;
2906 v0GammaC.SetProductionVertex(primVtxImprovedGamma);
2907
2908
2909 curGamma=curElecNeg+curElecPos;
2910 curGammaAt=curElecNegAt+curElecPosAt;
2911
2912 // invariant mass versus pt of K0short
2913
2914 Double_t chi2V0GammaC=100000.;
2915 if( v0GammaC.GetNDF() != 0) {
2916 chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
2917 }else{
2918 cout<< "ERROR::v0K0C.GetNDF()" << endl;
2919 }
2920
2921 if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
2922 if(fHistograms != NULL){
2923 fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
2924 fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
2925 fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
2926 fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
2927 fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
2928 fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
9c1cb6f7 2929 fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
2930 fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
2931
2932 new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()]) AliKFParticle(v0GammaC);
2933 fElectronRecalculatedv1.push_back(iTracks);
2934 fElectronRecalculatedv2.push_back(jTracks);
2935 }
2936 }
2937 }
2938 }
2939
2940 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
2941 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
2942 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);