]>
Commit | Line | Data |
---|---|---|
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 | ||
46 | using namespace std; | |
47 | ||
48 | ClassImp(AliAnalysisTaskPi0Reconstruction) | |
49 | ||
50 | ||
51 | //________________________________________________________________________ | |
52 | AliAnalysisTaskPi0Reconstruction::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 | //________________________________________________________________________ | |
85 | void AliAnalysisTaskPi0Reconstruction::SetBGHandler(AliConversionAODBGHandlerRP *bgHandler) | |
86 | { | |
87 | fBGHandler=bgHandler; | |
88 | fNCentralityBins=fBGHandler->GetNCentralityBins(); | |
89 | } | |
90 | //________________________________________________________________________ | |
91 | void AliAnalysisTaskPi0Reconstruction::SetDefaultBGHandler() | |
92 | { | |
93 | AliInfo("Setting up default BGHandler"); | |
94 | fBGHandler=new AliConversionAODBGHandlerRP(); | |
95 | fNCentralityBins=fBGHandler->GetNCentralityBins(); | |
96 | } | |
97 | ||
98 | //________________________________________________________________________ | |
99 | AliAnalysisTaskPi0Reconstruction::~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 | //________________________________________________________________________ | |
110 | void 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 | //________________________________________________________________________ | |
184 | void 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 | ||
212 | void 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 | //________________________________________________________________________ | |
243 | Bool_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 | ||
272 | Bool_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 | } | |
294 | return kFALSE; | |
295 | } | |
296 | ||
297 | //________________________________________________________________________ | |
298 | Bool_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 | //________________________________________________________________________ | |
317 | void 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 | //________________________________________________________________________ | |
370 | void 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 | //________________________________________________________________________ | |
380 | void 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 | ||
390 | void 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 | //________________________________________________________________________ | |
488 | Bool_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 | //________________________________________________________________________ | |
569 | Bool_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 | //________________________________________________________________________ | |
618 | void 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 | ///________________________________________________________________________ | |
632 | Bool_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 | } |