]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/GammaConv/AliAnalysisTaskPi0Reconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changes from gsi svn
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskPi0Reconstruction.cxx
CommitLineData
2eedd4ed 1#include <exception>
2#include "TChain.h"
3#include "TTree.h"
4#include "TVector3.h"
5#include "TCanvas.h"
6#include "TRandom3.h"
7#include "TH2.h"
8#include "TH1.h"
9#include "TH3.h"
10
11#include "AliLog.h"
12#include "AliAnalysisTask.h"
13#include "AliAnalysisManager.h"
14#include "AliAnalysisTaskPi0Reconstruction.h"
15
16#include "AliESDEvent.h"
17#include "AliESDInputHandler.h"
18
19#include "AliAODCentrality.h"
20#include "AliTPCPIDResponse.h"
21#include "AliAODFmdCluster.h"
22#include "AliKFParticle.h"
23#include "AliTracker.h"
24#include "AliV0Reader.h"
25#include "AliAODv0.h"
26#include "AliAODPid.h"
27#include "AliAODEvent.h"
28#include "AliAODTrack.h"
29#include "AliESDtrack.h"
30#include "AliAODHandler.h"
31#include "AliConversionAODBGHandlerRP.h"
32#include "AliV0Reader.h"
33#include "AliV0ReaderV1.h"
34#include "AliCDBEntry.h"
35#include "TObjArray.h"
36#include "AliEventplane.h"
37#include "AliCentrality.h"
38#include "AliVCluster.h"
39#include "AliESDCaloCluster.h"
40#include <iostream>
41
42#include <exception>
43
44// Author Daniel Lohner (Daniel.Lohner@cern.ch)
45
46using namespace std;
47
48ClassImp(AliAnalysisTaskPi0Reconstruction)
49
50
51//________________________________________________________________________
52AliAnalysisTaskPi0Reconstruction::AliAnalysisTaskPi0Reconstruction(const char *name) : AliV0ReaderV1(name),
53 fDeltaAODBranchName("GammaConv_gamma"),
54 fBGHandler(NULL),
55 Pi0MassRange(NULL),
56 kEventMixing(kTRUE),
57 fIsHeavyIon(kTRUE),
58 fRandomizer(NULL),
59 kUseSatelliteAODs(kTRUE),
60 fPHOSGammas(NULL),
61 fEMCALGammas(NULL),
62 fPi0Candidates(NULL),
63 fBGPi0s(NULL),
64 fMCTruePi0s(NULL),
65 fRapidityCut(0.9),
66 fAlphaCut(1),
67 kUseOnlyTaggedPhotons(kFALSE),
68 fNRandomEventsForBGRotation(15)
69{
70
71 fRandomizer= new TRandom3();
72
73 // If V0 Reader is available use it
74 SetUseAODConversionPhoton(kTRUE);
75
76 // Default Values
77 Pi0MassRange=new Double_t[2];
78 Pi0MassRange[0]=0.;
79 Pi0MassRange[1]=0.3;
80
81 DefineInput(0, TChain::Class());
82 DefineOutput(1, TList::Class());
83}
84//________________________________________________________________________
85void AliAnalysisTaskPi0Reconstruction::SetBGHandler(AliConversionAODBGHandlerRP *bgHandler)
86{
87 fBGHandler=bgHandler;
88 fNCentralityBins=fBGHandler->GetNCentralityBins();
89}
90//________________________________________________________________________
91void AliAnalysisTaskPi0Reconstruction::SetDefaultBGHandler()
92{
93 AliInfo("Setting up default BGHandler");
94 fBGHandler=new AliConversionAODBGHandlerRP();
95 fNCentralityBins=fBGHandler->GetNCentralityBins();
96}
97
98//________________________________________________________________________
99AliAnalysisTaskPi0Reconstruction::~AliAnalysisTaskPi0Reconstruction(){
100
101 if(fBGHandler){delete fBGHandler;
102 fBGHandler=0x0;}
103 if(fRandomizer){delete fRandomizer;fRandomizer=0x0;}
104
105 if(Pi0MassRange)delete[] Pi0MassRange;
106
107}
108
109//________________________________________________________________________
110void AliAnalysisTaskPi0Reconstruction::UserCreateOutputObjects()
111{
112 // Create histograms
113 // Called once
114 // Create the output container
115 /* if(fOutputList != NULL){
116 delete fOutputList;
117 fOutputList = NULL;
118 }
119 if(fOutputList == NULL){
120 fOutputList = new TList();
121 fOutputList->SetOwner(kTRUE);
122 }*/
123 AliV0ReaderV1::UserCreateOutputObjects();
124
125 TList *fPi0List=new TList();
126 fPi0List->SetName("Pi0Reconstruction");
127 fPi0List->SetOwner(kTRUE);
128 fOutputList->Add(fPi0List);
129
130
131 hPi0Cuts=new TH1F("Pi0Cand_Cuts","Pi0Cand_Cuts" ,10,-0.5,9.5);
132 fPi0List->Add(hPi0Cuts);
133 hPi0BGCuts=new TH1F("Pi0Cand_BG_Cuts","Pi0Cand_BG_Cuts" ,10,-0.5,9.5);
134 fPi0List->Add(hPi0BGCuts);
135 hPi0Alpha=new TH1F("Pi0Cand_Alpha" ,"Pi0Cand Alpha" ,200,0,1);
136 fPi0List->Add(hPi0Alpha);
137 hPi0OpeningAngle=new TH1F("Pi0Cand_OpeningAngle" ,"Pi0Cand OpeningAngle" ,400,0,TMath::Pi());
138 fPi0List->Add(hPi0OpeningAngle);
139 hPi0Rapidity=new TH1F("Pi0Cand_Rapidity" ,"Pi0Cand_Rapidity" ,100,-1,1);
140 fPi0List->Add(hPi0Rapidity);
141 fPi0List->Add(hPi0BGCuts);
142 hPi0AlphaRejected=new TH1F("Pi0Cand_Alpha_Rejected" ,"Pi0Cand Alpha Rejected" ,200,0,1);
143 fPi0List->Add(hPi0AlphaRejected);
144 hPi0OpeningAngleRejected=new TH1F("Pi0Cand_OpeningAngle_Rejected" ,"Pi0Cand OpeningAngle Rejected" ,400,0,TMath::Pi());
145 fPi0List->Add(hPi0OpeningAngleRejected);
146 hPi0RapidityRejected=new TH1F("Pi0Cand_Rapidity_Rejected" ,"Pi0Cand_Rapidity_Rejected" ,100,-1,1);
147 fPi0List->Add(hPi0RapidityRejected);
148
149 hPool=new TH3F("BGPool","BGPool",fNCentralityBins,-0.5,fNCentralityBins-0.5,7,-0.5,6.5,fBGHandler->GetNRPBins(),-0.5,fBGHandler->GetNRPBins()-0.5);
150 fPi0List->Add(hPool);
151
152 // MC
153
154 if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
155
156 Int_t kGCnXBinsSpectra = Int_t((Pi0MassRange[1]-Pi0MassRange[0])*500); //500 for range 0 - 1
157 Double_t kGCfirstXBinSpectra = Pi0MassRange[0];
158 Double_t kGClastXBinSpectra = Pi0MassRange[1];
159
160 Int_t kGCnYBinsSpectra = 250;
161 Double_t kGCfirstYBinSpectra = 0.;
162 Double_t kGClastYBinSpectra = 25.;
163
164 hPi0TRUE=new TH1F*[fNCentralityBins];
165 hPi0RECOTRUE=new TH2F*[fNCentralityBins];
166
167 for(Int_t m=0;m<fNCentralityBins;m++){
168
169 hPi0TRUE[m]=new TH1F(Form("%d_TRUE_Pi0_Pt",m),Form("%d_TRUE_Pi0_Pt",m),kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
170 hPi0TRUE[m]->Sumw2();
171 fPi0List->Add(hPi0TRUE[m]);
172
173 hPi0RECOTRUE[m]=new TH2F(Form("%d_RECOTRUE_Pi0_Pt",m),Form("%d_RECOTRUE_Pi0_Pt",m),kGCnXBinsSpectra, kGCfirstXBinSpectra, kGClastXBinSpectra,kGCnYBinsSpectra, kGCfirstYBinSpectra, kGClastYBinSpectra);
174 hPi0RECOTRUE[m]->Sumw2();
175 fPi0List->Add(hPi0RECOTRUE[m]);
176
177 }
178 }
179
180 PostData(1, fOutputList);
181}
182
183//________________________________________________________________________
184void AliAnalysisTaskPi0Reconstruction::UserExec(Option_t *)
185{
186
187
188 AliV0ReaderV1::UserExec("");
189
190 // Check if BG Handler exists
191
192 if(!fBGHandler){SetDefaultBGHandler();}
193
194 // Process Event Info
195
196 if(fEventIsSelected){
197
198 GetConversionGammas();
199
200 CalculatePi0Candidates();
201
202 CalculateBackground();
203
204 ProcessMCMesons();
205 }
206
207 PostData(1, fOutputList);
208}
209
210//________________________________________________________________________
211
212void AliAnalysisTaskPi0Reconstruction::ProcessMCMesons(){
213
214 if(fMCStack){
215
216 if(fMCTruePi0s == NULL){
217 fMCTruePi0s = new TClonesArray("TParticle",100);
218 }
219 fMCTruePi0s->Delete();//Reset the TClonesArray
220
221 for (Int_t part=0; part<fMCStack->GetNprimary(); part++) {
222
223 TParticle *fMCMother = fMCStack->Particle(part);
224
225 if(fMCMother){
226
227 // Process only Pions
228
229 if(fMCMother->GetPdgCode()==111){
230
231 // Histogram for pi0->2y efficiency (take acceptance)
232 if(IsMCMesonInReconstructionAcceptance(fMCMother)){
233 hPi0TRUE[fCentralityBin]->Fill(fMCMother->Pt());
234 fMCTruePi0s->AddAt(fMCMother,fMCTruePi0s->GetEntriesFast());
235 }
236 }
237 }
238 }
239 }
240}
241
242//________________________________________________________________________
243Bool_t AliAnalysisTaskPi0Reconstruction::IsMCMesonInReconstructionAcceptance(TParticle *fMCMother){
244
245 // Select only pi0->2y
246 if(fMCMother->GetNDaughters()!=2){return kFALSE;}
247
248 // PseudoRapidity Cut
249 if(TMath::Abs(fMCMother->Eta())>fRapidityCut){return kFALSE;}
250
251 Int_t NGammaAcceptance=0;
252
253 for(Int_t i=0;i<2;i++){
254 TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
255
256 // Is Daughter a Photon?
257 if(MDaughter->GetPdgCode()==22){
258
259 // Gamma Acceptance Cut
260 if(IsMCConversionGammaInAcceptance(MDaughter)){
261 NGammaAcceptance++;}
262 }
263 }
264
265 if(NGammaAcceptance!=2){return kFALSE;}
266
267 return kTRUE;
268}
269
270//________________________________________________________________________
271
272Bool_t AliAnalysisTaskPi0Reconstruction::IsTruePi0(AliAODConversionMother *pi0){
273
274 if(fMCStack){
275
276 AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(pi0->GetLabel(0)));
277 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(pi0->GetLabel(1)));
278
279 TParticle *mcgam0=gamma0->GetMCParticle(fMCStack);
280 TParticle *mcgam1=gamma1->GetMCParticle(fMCStack);
281 TParticle *mcpi0=NULL;
282
283 if(mcgam0&&mcgam1){
284 if(mcgam0->GetMother(0)==mcgam1->GetMother(0)){
285 mcpi0=fMCStack->Particle(mcgam0->GetMother(0));
286
287 if(mcpi0){
288 if(mcpi0->GetPdgCode()==111){return kTRUE;}
289 }
290 }
291 }
292
293 }
294return kFALSE;
295}
296
297//________________________________________________________________________
298Bool_t AliAnalysisTaskPi0Reconstruction::IsMCPhotonReconstructed(TParticle *MDaughter){
299
300 for(Int_t ii=0;ii<fConversionGammas->GetEntriesFast();ii++){
301
302 AliAODConversionPhoton *gamma=(AliAODConversionPhoton *)fConversionGammas->At(ii);
303
304 TParticle *fMCParticle=gamma->GetMCParticle(fMCStack);
305
306 if(fMCParticle){
307
308 if(fMCParticle->GetDaughter(0)==MDaughter->GetDaughter(0)&&fMCParticle->GetDaughter(1)==MDaughter->GetDaughter(1)){
309 return kTRUE;
310 }
311 }
312 }
313 return kFALSE;
314}
315
316//________________________________________________________________________
317void AliAnalysisTaskPi0Reconstruction::CalculatePi0Candidates(){
318
319 if(fPi0Candidates == NULL){
320 fPi0Candidates = new TClonesArray("AliAODConversionMother",100);
321 }
322 fPi0Candidates->Delete();//Reset the TClonesArray
323
324 if(fConversionGammas->GetEntriesFast()>1){
325
326 for(Int_t firstGammaIndex=0;firstGammaIndex<fConversionGammas->GetEntriesFast()-1;firstGammaIndex++){
327
328 AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(firstGammaIndex));
329
330 // Combine Photons
331
332 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fConversionGammas->GetEntriesFast();secondGammaIndex++){
333
334 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fConversionGammas->At(secondGammaIndex));
335
336 //Check for same Electron ID
337 if(gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelNegative()
338 ||gamma0->GetTrackLabelNegative()==gamma1->GetTrackLabelPositive()||gamma0->GetTrackLabelPositive()==gamma1->GetTrackLabelNegative())continue;
339
340 AliAODConversionMother *pi0cand=new AliAODConversionMother(gamma0,gamma1);
341 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
342
343
344 if((CheckPi0(pi0cand))){
345 if (pi0cand->M() > Pi0MassRange[0] && pi0cand->M() < Pi0MassRange[1] ){
346
347 // Add Pi0 to Stack
348 new((*fPi0Candidates)[fPi0Candidates->GetEntriesFast()]) AliAODConversionMother(*pi0cand);
349
350 // Test MC truth
351 if(IsTruePi0(pi0cand)){
352 hPi0RECOTRUE[fCentralityBin]->Fill(pi0cand->M(),pi0cand->Pt());
353 }
354
355 }
356 }
357
358 delete pi0cand;
359 pi0cand=0x0;
360
361 }
362
363 }
364
365 }
366}
367
368
369//________________________________________________________________________
370void AliAnalysisTaskPi0Reconstruction::Terminate(Option_t *)
371{
372
373 // Draw result to the screen
374 fOutputList->Print();
375 // Called once at the end of the query
376}
377
378
379//________________________________________________________________________
380void AliAnalysisTaskPi0Reconstruction::RotateParticle(AliAODConversionPhoton *gamma,Double_t angle){
381 Int_t fNDegreesPMBackground=15;
382 Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
383 Double_t rotationValue = fRandomizer->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
384 gamma->RotateZ(rotationValue);
385
386}
387
388//________________________________________________________________________
389
390void AliAnalysisTaskPi0Reconstruction::CalculateBackground(){
391
392 if(fBGPi0s == NULL){
393 fBGPi0s = new TClonesArray("AliAODConversionMother",100);
394 }
395 fBGPi0s->Delete();//Reset the TClonesArray
396
397 AliAODConversionMother *BGcandidate=NULL;
398 TClonesArray *currentEventGammas=fConversionGammas;
399
400 //Rotation Method
401 if(!kEventMixing){
402
403 // Correct for the number of rotations
404 // BG is for rotation the same, except for factor NRotations
405 Double_t weight=1./Double_t(fNRandomEventsForBGRotation);
406
407 for(Int_t firstGammaIndex=0;firstGammaIndex<currentEventGammas->GetEntriesFast();firstGammaIndex++){
408
409 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(currentEventGammas->At(firstGammaIndex));
410
411 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<currentEventGammas->GetEntriesFast();secondGammaIndex++){
412 for(Int_t nRandom=0;nRandom<fNRandomEventsForBGRotation;nRandom++){
413
414 AliAODConversionPhoton *gamma2=dynamic_cast<AliAODConversionPhoton*>(currentEventGammas->At(secondGammaIndex));
415
416 RotateParticle(gamma2);
417
418 BGcandidate=new AliAODConversionMother(gamma1,gamma2);
419
420 if(CheckPi0(BGcandidate,kFALSE)){
421
422 if (BGcandidate->M() > Pi0MassRange[0] && BGcandidate->M() < Pi0MassRange[1] ){
423
424 new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(*BGcandidate);
425
426 dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
427 }
428 }
429 delete BGcandidate;
430 BGcandidate=NULL;
431
432 }
433
434 }
435 }
436 }
437
438 else{
439 // Do Event Mixing
440
441 Int_t psibin=fBGHandler->GetRPBinIndex(fEPAngle);
442
443 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler->GetNBGEvents(psibin,fBGHandler->GetZBinIndex(fVertexZ),fCentralityBin);nEventsInBG++){
444
445 AliGammaConversionPhotonVector *previousEventGammas = fBGHandler->GetBGGoodGammas(psibin,fBGHandler->GetZBinIndex(fVertexZ),fCentralityBin,nEventsInBG);
446
447 if(previousEventGammas){
448
449 // test weighted background
450 Double_t weight=1.0;///Double_t(fBGHandler->GetNBGEvents(psibin,zbin,fCentralityBin));
451 // Correct for the number of eventmixing:
452 // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1)) using sum formula sum(i)=N*(N-1)/2 -> N*(N-1)/2
453 // real combinations (since you cannot combine a photon with its own)
454 // but BG leads to N_{a}*N_{b} combinations
455 weight*=0.5*(Double_t(currentEventGammas->GetEntriesFast()-1))/Double_t(previousEventGammas->size());
456
457 for(Int_t iCurrent=0;iCurrent<currentEventGammas->GetEntriesFast();iCurrent++){
458 AliAODConversionPhoton *currentEventGamma = (AliAODConversionPhoton*)(currentEventGammas->At(iCurrent));
459 for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
460
461 AliAODConversionPhoton *previousEventGamma = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
462
463 BGcandidate=new AliAODConversionMother(previousEventGamma,currentEventGamma);
464
465 if(CheckPi0(BGcandidate,kFALSE)){
466
467 if (BGcandidate->M() > Pi0MassRange[0] && BGcandidate->M() < Pi0MassRange[1] ){
468
469 new((*fBGPi0s)[fBGPi0s->GetEntriesFast()]) AliAODConversionMother(*BGcandidate);
470 dynamic_cast<AliAODConversionMother*>(fBGPi0s->At(fBGPi0s->GetEntriesFast()-1))->SetWeight(weight);
471 }
472 }
473 delete BGcandidate;
474 BGcandidate=NULL;
475 }
476 }
477 }
478 }
479
480 // Add Current Event to BG POOL
481 fBGHandler->AddEvent(fConversionGammas,fEPAngle,fVertexZ,Int_t(fCentrality));
482 hPool->Fill(fCentralityBin,fBGHandler->GetZBinIndex(fVertexZ),fBGHandler->GetRPBinIndex(fEPAngle));
483
484 }
485}
486
487//________________________________________________________________________
488Bool_t AliAnalysisTaskPi0Reconstruction::GetPHOSGammas(){
489
490 // TClones Array for local working copy
491
492 /*
493 if(fPHOSGammas == NULL){
494 fPHOSGammas = new TClonesArray("AliAODConversionPhoton",100);
495 }
496 fPHOSGammas->Delete();//Reset the TClonesArray
497
498 //AliESDCaloCluster *clu=NULL;
499 AliVCluster *clu=NULL;
500 TLorentzVector pPHOS;
501 AliAODConversionPhoton *gammaPHOS=NULL;
502
503 Double_t vtx[3];
504 vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
505 vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
506 vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
507
508
509 for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++){
510 clu = fESDEvent->GetCaloCluster(i);
511
512 // Cuts from AliAnalysisTaskCaloConv
513 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
514 clu ->GetMomentum(pPHOS ,vtx);
515 if(pPHOS.Energy()<0.25){
516 continue ;
517 }
518 if(clu->GetNCells()<=2){
519 continue ;
520 }
521
522 Bool_t isNeutral = kTRUE ;
523 Bool_t isDispOK = kTRUE ;
524 Bool_t isTOFOK = kTRUE ;
525 Int_t iMod,iX,iZ ;
526
527 isNeutral = clu->GetEmcCpvDistance()>5. ; //To be improved
528 isDispOK = kFALSE ;
529 Double_t l0=clu->GetM02(),l1=clu->GetM20() ;
530 if(l1>= 0 && l0>= 0 && l1 < 0.1 && l0 < 0.1) isDispOK=kFALSE ;
531 if(l1>= 0 && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) isDispOK=kTRUE ;
532 if(l1>= 0 && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) isDispOK=kFALSE ;
533 if(l1>= 0 && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) isDispOK=kFALSE ;
534 if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) isDispOK=kTRUE ;
535 if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) isDispOK=kFALSE ;
536 if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) isDispOK=kFALSE ;
537 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) isDispOK=kTRUE ;
538 if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) isDispOK=kTRUE ;
539 if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) isDispOK=kTRUE ;
540
541 Float_t xyz[3] = {0,0,0};
542 clu->GetPosition(xyz); //Global position in ALICE system
543 TVector3 global(xyz) ;
544 Int_t relid[4] ;
545 if(!fPHOSgeom->GlobalPos2RelId(global,relid)){
546 printf("PHOS_beyond: x=%f, y=%f, z=%f \n",xyz[0],xyz[1],xyz[2]) ;
547 continue ;
548 }
549 iMod=relid[0] ;
550 iX=relid[2];
551 iZ=relid[3] ;
552 if(!IsGoodChannel("PHOS",iMod,iX,iZ))
553 continue ;
554
555 // Add selected cluster to Photonqueue
556 gammaPHOS=new AliAODConversionPhoton();
557 gammaPHOS->SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
558 new((*fPHOSGammas)[fPHOSGammas->GetEntriesFast()]) AliAODConversionPhoton(*gammaPHOS);
559 }
560
561 if(fPHOSGammas->GetEntries()){return kTRUE;}
562*/
563 return kFALSE;
564
565
566}
567
568//________________________________________________________________________
569Bool_t AliAnalysisTaskPi0Reconstruction::GetConversionGammas(){
570
571 // ESD mode
572
573 if(fESDEvent){
574
575 fConversionGammas=GetReconstructedGammas();
576
577 if(!fConversionGammas)AliError("No Gammas from V0 Reader");
578 }
579
580
581 // AOD mode
582
583 if(fAODEvent&&kUseSatelliteAODs){
584
585 if(fConversionGammas == NULL){
586 fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
587 }
588 fConversionGammas->Delete();//Reset the TClonesArray
589
590 //Get Gammas from satellite AOD gamma branch
591
592 AliAODConversionPhoton *gamma=0x0;
593
594 TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
595 if(!fInputGammas){FindDeltaAODBranchName();
596 fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
597 if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
598
599 // Apply Selection Cuts to Gammas and create local working copy
600 if(fInputGammas){
601 for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
602 gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
603 if(gamma){
604 if(IsGammaCandidate(gamma)){
605 new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
606 }
607 }
608 }
609
610 }
611
612 if(fConversionGammas->GetEntries()){return kTRUE;}
613
614 return kFALSE;
615}
616
617//________________________________________________________________________
618void AliAnalysisTaskPi0Reconstruction::FindDeltaAODBranchName(){
619 TList *list=fAODEvent->GetList();
620 for(Int_t ii=0;ii<list->GetEntries();ii++){
621 TString name((list->At(ii))->GetName());
622 if(name.BeginsWith("GammaConv")&&name.EndsWith("gamma")){
623 fDeltaAODBranchName=name;
624 AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
625 return;
626 }
627 }
628}
629
630
631///________________________________________________________________________
632Bool_t AliAnalysisTaskPi0Reconstruction::CheckPi0(AliAODConversionMother *pi0,Bool_t IsSignal)
633{
634
635 Bool_t passcuts=kTRUE;
636
637 TH1 *hist=0x0;
638
639 if(IsSignal){hist=hPi0Cuts;
640 }
641 else{hist=hPi0BGCuts;}
642
643 // Undefined Rapidity -> Floating Point exception
644 if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
645 hist->Fill(0);
646 passcuts=kFALSE;
647 }
648 else{
649 // PseudoRapidity Cut
650 if(TMath::Abs(pi0->PseudoRapidity())>fRapidityCut){
651 hist->Fill(1);
652 passcuts=kFALSE;
653 }
654 }
655
656 // Opening Angle Cut
657 if(pi0->GetOpeningAngle()<0.005){
658 hist->Fill(2);
659 passcuts=kFALSE;
660 }
661
662 // Alpha (Energy Asymmetry) Cut
663 if(pi0->GetAlpha()>fAlphaCut){
664 hist->Fill(3);
665 passcuts=kFALSE;
666 }
667
668 // Fill Histograms
669
670 if(passcuts){
671 hist->Fill(9);
672 if(IsSignal){
673 hPi0Rapidity->Fill(pi0->PseudoRapidity());
674 hPi0Alpha->Fill(pi0->GetAlpha());
675 hPi0OpeningAngle->Fill(pi0->GetOpeningAngle());
676
677 }
678 }
679 else{
680 hist->Fill(8);
681 if(IsSignal){
682 hPi0RapidityRejected->Fill(pi0->PseudoRapidity());
683 hPi0AlphaRejected->Fill(pi0->GetAlpha());
684 hPi0OpeningAngleRejected->Fill(pi0->GetOpeningAngle());
685 }
686 }
687
688 return passcuts;
689}