]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/PHOSTasks/PHOS_Correlation/AliPHOSCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates on V0 correlation
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_Correlation / AliPHOSCorrelations.cxx
CommitLineData
67ef08bd 1/**************************************************************************
2 * Copyright(c) 1998-2014, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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// Analysis task for identified PHOS cluster from pi0 and take korrelation betwen hadron-pi0 angel's.
17// Authors: Daniil Ponomarenko (Daniil.Ponomarenko@cern.ch)
18// Dmitry Blau
19// 07-Feb-2014
20
21#include "THashList.h"
22#include "TObjArray.h"
23#include "TClonesArray.h"
24#include "TH1F.h"
25#include "TH2F.h"
26#include "TH2I.h"
27#include "TH3F.h"
28#include "TH3D.h"
29#include "TMath.h"
30#include "TVector3.h"
31#include "TProfile.h"
32
33#include "AliAnalysisManager.h"
34#include "AliAnalysisTaskSE.h"
35#include "AliPHOSCorrelations.h"
36#include "AliPHOSGeometry.h"
37#include "AliESDEvent.h"
38#include "AliESDCaloCluster.h"
39#include "AliESDVertex.h"
40#include "AliESDtrackCuts.h"
41#include "AliESDtrack.h"
42#include "AliAODTrack.h"
43#include "AliVTrack.h"
44#include "AliPID.h"
45#include "AliTriggerAnalysis.h"
46#include "AliPIDResponse.h"
47#include "AliPHOSEsdCluster.h"
48#include "AliCDBManager.h"
49#include "AliPHOSCalibData.h"
50#include "AliCentrality.h"
51#include "AliAnalysisTaskSE.h"
52#include "AliEventplane.h"
53#include "AliOADBContainer.h"
54#include "AliAODEvent.h"
55#include "AliAODCaloCells.h"
56#include "AliAODCaloCluster.h"
57#include "AliCaloPhoton.h"
58#include "AliAODVertex.h"
59
60ClassImp(AliPHOSCorrelations)
61
62//_______________________________________________________________________________
63AliPHOSCorrelations::AliPHOSCorrelations()
64:AliAnalysisTaskSE(),
65 fOutputContainer(0x0),
66 fMinClusterEnergy(0.3),
67 fMinBCDistance(0),
68 fMinNCells(3),
69 fMinM02(0.2),
70 fTOFCutEnabled(1),
71 fTOFCut(100.e-9),
72 fNVtxZBins(1),
73 fCentEdges(10),
74 fCentNMixed(),
75 fNEMRPBins(9),
76 fAssocBins(),
77 fEventsMultiplicityOnly(false),
78 fCheckHibridGlobal(kOnlyHibridTracks),
79 fPeriod(kUndefinedPeriod),
80 fInternalTriggerSelection(kNoSelection),
81 fMaxAbsVertexZ(10.),
82 fManualV0EPCalc(false),
83 fCentCutoffDown(0.),
84 fCentCutoffUp(90),
85 fMassInvMin(0.125),
86 fMassInvMax(0.145),
87 fEvent(0x0),
88 fEventESD(0x0),
89 fEventAOD(0x0),
90 fESDtrackCuts(0x0),
91 fRunNumber(-999),
92 fInternalRunNumber(0),
93 fMultV0(0x0),
94 fV0Cpol(0.),fV0Apol(0.),
95 fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"),
96 fVertexVector(),
97 fVtxBin(0),
98 fCentralityEstimator("V0M"),
99 fCentrality(0.),
100 fCentBin(0),
101 fHaveTPCRP(0),
102 fRP(0.),
103 fEMRPBin(0),
104 fCaloPhotonsPHOS(0x0),
105 fTracksTPC(0x0),
106 fCaloPhotonsPHOSLists(0x0),
107 fTracksTPCLists(0x0)
108{
109 //Deafult constructor, no memory allocations here
110
111}
112
113//_______________________________________________________________________________
114AliPHOSCorrelations::AliPHOSCorrelations(const char *name, Period period)
115:AliAnalysisTaskSE(name),
116 fOutputContainer(0x0),
117 fMinClusterEnergy(0.3),
118 fMinBCDistance(0),
119 fMinNCells(3),
120 fMinM02(0.2),
121 fTOFCutEnabled(1),
122 fTOFCut(100.e-9),
123 fNVtxZBins(1),
124 fCentEdges(10),
125 fCentNMixed(),
126 fNEMRPBins(9),
127 fAssocBins(),
128 fEventsMultiplicityOnly(false),
129 fCheckHibridGlobal(kOnlyHibridTracks),
130 fPeriod(period),
131 fInternalTriggerSelection(kNoSelection),
132 fMaxAbsVertexZ(10.),
133 fManualV0EPCalc(false),
134 fCentCutoffDown(0.),
135 fCentCutoffUp(90),
136 fMassInvMin(0.125),
137 fMassInvMax(0.145),
138 fEvent(0x0),
139 fEventESD(0x0),
140 fEventAOD(0x0),
141 fESDtrackCuts(0x0),
142 fRunNumber(-999),
143 fInternalRunNumber(0),
144 fMultV0(0x0),
145 fV0Cpol(0.),fV0Apol(0.),
146 fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"),
147 fVertexVector(),
148 fVtxBin(0),
149 fCentralityEstimator("V0M"),
150 fCentrality(0.),
151 fCentBin(0),
152 fHaveTPCRP(0),
153 fRP(0.),
154 fEMRPBin(0),
155 fCaloPhotonsPHOS(0x0),
156 fTracksTPC(0x0),
157 fCaloPhotonsPHOSLists(0x0),
158 fTracksTPCLists(0x0)
159
160{
161 // Constructor
162 // Output slots #0 write into a TH1 container
163 DefineOutput(1,THashList::Class());
164
165 const Int_t nPtAssoc=10 ;
166 Double_t ptAssocBins[nPtAssoc]={0.,0.5,1.0,1.5,2.0,3.,5.,7.,10.,16} ;
167 fAssocBins.Set(nPtAssoc,ptAssocBins) ;
168
169 const int nbins = 9;
170 Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
171 TArrayD centEdges(nbins+1, edges);
172 Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
173 //Int_t nMixed[nbins] = {100,100,100,100,100,100,100,100,100};
174 TArrayI centNMixed(nbins, nMixed);
175 SetCentralityBinning(centEdges, centNMixed);
176
177 fVertex[0]=0; fVertex[1]=0; fVertex[2]=0;
178}
179//_______________________________________________________________________________
180AliPHOSCorrelations::~AliPHOSCorrelations()
181{
182
183 if(fCaloPhotonsPHOS){
184 delete fCaloPhotonsPHOS;
185 fCaloPhotonsPHOS=0x0;
186 }
187
188 if(fTracksTPC){
189 delete fTracksTPC;
190 fTracksTPC=0x0;
191 }
192
193 if(fCaloPhotonsPHOSLists){
194 fCaloPhotonsPHOSLists->SetOwner() ;
195 delete fCaloPhotonsPHOSLists;
196 fCaloPhotonsPHOSLists=0x0;
197 }
198
199 if(fTracksTPCLists){
200 fTracksTPCLists->SetOwner() ;
201 delete fTracksTPCLists;
202 fTracksTPCLists=0x0 ;
203 }
204
205 if( fESDtrackCuts){
206 delete fESDtrackCuts;
207 fESDtrackCuts=0x0 ;
208 }
209
210 if(fOutputContainer){
211 delete fOutputContainer;
212 fOutputContainer=0x0;
213 }
214
215}
216//_______________________________________________________________________________
217void AliPHOSCorrelations::UserCreateOutputObjects()
218{
219 // Create histograms
220 // Called once
221 const Int_t nRuns=200 ;
222
223 // Create histograms
224 if(fOutputContainer != NULL) { delete fOutputContainer; }
225 fOutputContainer = new THashList();
226 fOutputContainer->SetOwner(kTRUE);
227
228 //Event selection
229 fOutputContainer->Add(new TH1F("hTriggerPassedEvents","Event selection passed Cuts", 20, 0., 20.) );
230
231 fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", kTotalSelected+3, 0., double(kTotalSelected+3))) ;
232 if (!fEventsMultiplicityOnly)
233 {
234 fOutputContainer->Add(new TH2F("hSelEvents","Event selection", kTotalSelected+1, 0., double(kTotalSelected+1), nRuns,0.,float(nRuns))) ;
235 fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
236 fOutputContainer->Add(new TH2F("phiRPflat","RP distribution with TPC flat", 100, 0., 2.*TMath::Pi(),20,0.,100.)) ;
237 fOutputContainer->Add(new TH2F("massWindow","mean & sigma", 100,0.1,0.18,100,0.,0.5));
238 }
239
240 // Set hists, with track's and cluster's angle distributions.
241 if (!fEventsMultiplicityOnly)
242 {
243 SetHistEtaPhi();
244 SetHistCutDistribution();
245 SetHistPtAssoc();
246 }
247
248 // Setup photon lists
249 Int_t kapacity = fNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
250 fCaloPhotonsPHOSLists = new TObjArray(kapacity);
251 fCaloPhotonsPHOSLists->SetOwner();
252
253 fTracksTPCLists = new TObjArray(kapacity);
254 fTracksTPCLists->SetOwner();
255
256 PostData(1, fOutputContainer);
257}
258//_______________________________________________________________________________
259void AliPHOSCorrelations::SetHistEtaPhi()
260// Set hists, with track's and cluster's angle distributions.
261{
262// cout<<"\nSetting output SetHist_Eta_Phi()...";
263
264 Float_t pi = TMath::Pi();
265
266 //===
267 fOutputContainer->Add(new TH2F("clu_phieta","Cluster's Phi & Eta distribution", 300, double(-1.8), double(-0.6), 300, double(-0.2), double(0.2) ) );
268 TH2F * h = static_cast<TH2F*>(fOutputContainer->Last()) ;
269 h->GetXaxis()->SetTitle("#phi [rad]");
270 h->GetYaxis()->SetTitle("#eta [rad]");
271
272 //===
273 fOutputContainer->Add(new TH2F("clusingle_phieta","Cluster's single Phi & Eta distribution", 300, double(-1.8), double(-0.6), 300, double(-0.2), double(0.2) ) );
274 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
275 h->GetXaxis()->SetTitle("#phi [rad]");
276 h->GetYaxis()->SetTitle("#eta [rad]");
277
278 //===
279 fOutputContainer->Add(new TH2F("track_phieta","TPC track's Phi & Eta distribution", 100, double(-0.3), double(2*pi+0.3), 200, double(-0.9), double(0.9) ) );
280 h = static_cast<TH2F*>(fOutputContainer->FindObject("track_phieta")) ;
281 h->GetXaxis()->SetTitle("#phi [rad]");
282 h->GetYaxis()->SetTitle("#eta [rad]");
283
284// cout<<" OK!"<<endl;
285}
286//_______________________________________________________________________________
287void AliPHOSCorrelations::SetHistCutDistribution()
288// Set other histograms.
289{
290// cout<<"\nSetting output SetHist_CutDistribution...";
291
292 Int_t PtMult = 100;
293 Double_t PtMin = 0.;
294 Double_t PtMax = 20.;
295
296
297 // Real ++++++++++++++++++++++++++++++
298
299 fOutputContainer->Add(new TH2F("all_mpt"," Only standard cut's ", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
300 TH2F * h = static_cast<TH2F*>(fOutputContainer->Last()) ;
301 h->GetXaxis()->SetTitle("Mass [GeV]");
302 h->GetYaxis()->SetTitle("Pt [GEV]");
303
304 fOutputContainer->Add(new TH2F("cpv_mpt"," CPV cut ", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
305 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
306 h->GetXaxis()->SetTitle("Mass [GeV]");
307 h->GetYaxis()->SetTitle("Pt [GEV]");
308
309 fOutputContainer->Add(new TH2F("disp_mpt"," Disp cut ", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
310 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
311 h->GetXaxis()->SetTitle("Mass [GeV]");
312 h->GetYaxis()->SetTitle("Pt [GEV]");
313
314 fOutputContainer->Add(new TH2F("both_mpt"," Both cuts (CPV + Disp) ", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
315 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
316 h->GetXaxis()->SetTitle("Mass [GeV]");
317 h->GetYaxis()->SetTitle("Pt [GEV]");
318
319
320 // MIX +++++++++++++++++++++++++
321
322 fOutputContainer->Add(new TH2F("mix_all_mpt"," Only standard cut's (mix)", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
323 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
324 h->GetXaxis()->SetTitle("Mass [GeV]");
325 h->GetYaxis()->SetTitle("Pt [GEV]");
326
327 fOutputContainer->Add(new TH2F("mix_cpv_mpt"," CPV cut (mix)", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
328 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
329 h->GetXaxis()->SetTitle("Mass [GeV]");
330 h->GetYaxis()->SetTitle("Pt [GEV]");
331
332 fOutputContainer->Add(new TH2F("mix_disp_mpt"," Disp cut (mix)", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
333 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
334 h->GetXaxis()->SetTitle("Mass [GeV]");
335 h->GetYaxis()->SetTitle("Pt [GEV]");
336
337 fOutputContainer->Add(new TH2F("mix_both_mpt"," Both cuts (CPV + Disp) (mix)", 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
338 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
339 h->GetXaxis()->SetTitle("Mass [GeV]");
340 h->GetYaxis()->SetTitle("Pt [GEV]");
341
342
343 // Calibration Pi0peak {REAL}
344 for(Int_t mod=1; mod<4; mod++){
345 fOutputContainer->Add(new TH2F(Form("both%d_mpt",mod),Form("Both cuts (CPV + Disp) mod[%d]",mod), 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
346 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
347 h->GetXaxis()->SetTitle("Mass [GeV]");
348 h->GetYaxis()->SetTitle("Pt [GEV]");
349
350 // Calibration Pi0peak {MIX}
351 fOutputContainer->Add(new TH2F(Form("mix_both%d_mpt",mod),Form(" Both cuts (CPV + Disp) mod[%d]",mod), 100, fMassInvMin, fMassInvMax, PtMult, PtMin, PtMax ) );
352 h = static_cast<TH2F*>(fOutputContainer->FindObject("mix_both1_mpt")) ;
353 h->GetXaxis()->SetTitle("Mass [GeV]");
354 h->GetYaxis()->SetTitle("Pt [GEV]");
355
356 }
357
358// cout<<" OK!"<<endl;
359}
360//_______________________________________________________________________________
361void AliPHOSCorrelations::SetHistPtAssoc()
362{
363 Double_t pi = TMath::Pi();
364
365 Int_t PhiMult = 100;
366 Float_t PhiMin = -0.5*pi;
367 Float_t PhiMax = 1.5*pi;
368 Int_t EtaMult = 20;
369 Float_t EtaMin = -1.;
370 Float_t EtaMax = 1.;
371 Int_t PtTrigMult = 100;
372 Float_t PtTrigMin = 0.;
373 Float_t PtTrigMax = 20.;
374
375 TString spid[4]={"all","cpv","disp","both"} ;
376
377 for (int i = 0; i<fAssocBins.GetSize()-1; i++){
378 for(Int_t ipid=0; ipid<4; ipid++){
379 fOutputContainer->Add(new TH3F(Form("%s_ptphieta_ptAssoc_%3.1f",spid[ipid].Data(),fAssocBins.At(i+1)),
380 Form("%s_ptphieta_ptAssoc_%3.1f",spid[ipid].Data(),fAssocBins.At(i+1)),
381 PtTrigMult, PtTrigMin, PtTrigMax, PhiMult, PhiMin, PhiMax, EtaMult, EtaMin, EtaMax ) );
382 TH3F * h = static_cast<TH3F*>(fOutputContainer->Last()) ;
383 h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
384 h->GetYaxis()->SetTitle("#phi [rad]");
385 h->GetZaxis()->SetTitle("#eta [rad]");
386
387 fOutputContainer->Add(new TH3F(Form("mix_%s_ptphieta_ptAssoc_%3.1f",spid[ipid].Data(),fAssocBins.At(i+1)),
388 Form("Mixed %s_ptphieta_ptAssoc_%3.1f",spid[ipid].Data(),fAssocBins.At(i+1)),
389 PtTrigMult, PtTrigMin, PtTrigMax, PhiMult, PhiMin, PhiMax, EtaMult, EtaMin, EtaMax ) );
390 h = static_cast<TH3F*>(fOutputContainer->Last()) ;
391 h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
392 h->GetYaxis()->SetTitle("#phi [rad]");
393 h->GetZaxis()->SetTitle("#eta [rad]");
394
395 }
396 }
397}
398//_______________________________________________________________________________
399void AliPHOSCorrelations::UserExec(Option_t *)
400{
401// Main loop, called for each event analyze ESD/AOD
402 // Step 0: Event Objects
403 LogProgress(0);
404 fEvent = InputEvent();
405 if( ! fEvent ) {
406 AliError("Event could not be retrieved");
407 PostData(1, fOutputContainer);
408 return ;
409 }
410
411 fEventESD = dynamic_cast<AliESDEvent*>(fEvent);
412 fEventAOD = dynamic_cast<AliAODEvent*>(fEvent);
413
414 {
415 FillHistogram("hTriggerPassedEvents", 0);
416
417 Bool_t isMB = (fEvent->GetTriggerMask() & (ULong64_t(1)<<1));
418 Bool_t isCentral = (fEvent->GetTriggerMask() & (ULong64_t(1)<<4));
419 Bool_t isSemiCentral = (fEvent->GetTriggerMask() & (ULong64_t(1)<<7));
420
421 if (isMB) FillHistogram("hTriggerPassedEvents", 2.);
422 if (isCentral) FillHistogram("hTriggerPassedEvents", 3.);
423 if (isSemiCentral) FillHistogram("hTriggerPassedEvents", 4.);
424 }
425
426 // For first event from data only:
427 if( fRunNumber<0)
428 {
429// cout<<"Mean: "<< fMassInvMin + (fMassInvMax - fMassInvMin) /2. << " Sigma: "<< (fMassInvMax - fMassInvMin) /2.<<endl;
430 FillHistogram("massWindow", fMassInvMin + (fMassInvMax - fMassInvMin) /2., (fMassInvMax - fMassInvMin) /2.);
431 }
432
433 // Step 1(done once):
434 if( fRunNumber != fEvent->GetRunNumber() )
435 {
436 fRunNumber = fEvent->GetRunNumber();
437 fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
438 SetESDTrackCuts();
439 }
440 LogProgress(1);
441
442 if( RejectTriggerMaskSelection() )
443 {
444 PostData(1, fOutputContainer);
445 return; // Reject!
446 }
447 LogProgress(2);
448
449 // Step 2: Vertex
450 // fVertex, fVertexVector, fVtxBin
451 SetVertex();
452 if( RejectEventVertex() )
453 {
454 PostData(1, fOutputContainer);
455 return; // Reject!
456 }
457 LogProgress(3);
458
459 // Step 3: Centrality
460 // fCentrality, fCentBin
461 SetCentrality();
462 if( RejectEventCentrality() )
463 {
464 PostData(1, fOutputContainer);
465 return; // Reject!
466 }
467 FillHistogram("hCentrality",fCentrality,fInternalRunNumber-0.5) ;
468 LogProgress(4);
469
470 // Step 4: Reaction Plane
471 // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
472 EvalReactionPlane();
473 fEMRPBin = GetRPBin();
474 LogProgress(5);
475
476 // Step 5: Event Photons (PHOS Clusters) selectionMakeFlat
477 SelectPhotonClusters();
478 if( ! fCaloPhotonsPHOS->GetEntriesFast() ) LogSelection(kHasPHOSClusters, fInternalRunNumber);
479 LogProgress(6);
480
481 // Step 6: Event Associated particles (TPC Tracks) selection
482 SelectAccosiatedTracks();
483
484 if( ! fTracksTPC->GetEntriesFast() )
485 LogSelection(kHasTPCTracks, fInternalRunNumber);
486 LogSelection(kTotalSelected, fInternalRunNumber);
487 LogProgress(7);
488
489 // Step 7: Consider pi0 (photon/cluster) pairs.
490 ConsiderPi0s();
491
492 // Step 8; Mixing
493 ConsiderPi0sMix();
494
495 ConsiderTracksMix();
496 //this->ConsiderPi0sTracksMix(); // Read how make mix events!
497 LogProgress(8);
498
499 // Step 9: Make TPC's mask
500 FillTrackEtaPhi();
501 LogProgress(9);
502
503 // Step 10: Update lists
504 UpdatePhotonLists();
505 UpdateTrackLists();
506
507 LogProgress(10);
508
509 // Post output data.
510 PostData(1, fOutputContainer);
511}
512//_______________________________________________________________________________
513void AliPHOSCorrelations::SetESDTrackCuts()
514{
515 if( fEventESD ) {
516 // Create ESD track cut
517 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
518 //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
519 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
520 }
521}
522//_______________________________________________________________________________
523Int_t AliPHOSCorrelations::ConvertToInternalRunNumber(Int_t run){
524 if(fPeriod== kLHC11h){
525 switch(run)
526 {
527 case 170593 : return 179 ;
528 case 170572 : return 178 ;
529 case 170556 : return 177 ;
530 case 170552 : return 176 ;
531 case 170546 : return 175 ;
532 case 170390 : return 174 ;
533 case 170389 : return 173 ;
534 case 170388 : return 172 ;
535 case 170387 : return 171 ;
536 case 170315 : return 170 ;
537 case 170313 : return 169 ;
538 case 170312 : return 168 ;
539 case 170311 : return 167 ;
540 case 170309 : return 166 ;
541 case 170308 : return 165 ;
542 case 170306 : return 164 ;
543 case 170270 : return 163 ;
544 case 170269 : return 162 ;
545 case 170268 : return 161 ;
546 case 170267 : return 160 ;
547 case 170264 : return 159 ;
548 case 170230 : return 158 ;
549 case 170228 : return 157 ;
550 case 170208 : return 156 ;
551 case 170207 : return 155 ;
552 case 170205 : return 154 ;
553 case 170204 : return 153 ;
554 case 170203 : return 152 ;
555 case 170195 : return 151 ;
556 case 170193 : return 150 ;
557 case 170163 : return 149 ;
558 case 170162 : return 148 ;
559 case 170159 : return 147 ;
560 case 170155 : return 146 ;
561 case 170152 : return 145 ;
562 case 170091 : return 144 ;
563 case 170089 : return 143 ;
564 case 170088 : return 142 ;
565 case 170085 : return 141 ;
566 case 170084 : return 140 ;
567 case 170083 : return 139 ;
568 case 170081 : return 138 ;
569 case 170040 : return 137 ;
570 case 170038 : return 136 ;
571 case 170036 : return 135 ;
572 case 170027 : return 134 ;
573 case 169981 : return 133 ;
574 case 169975 : return 132 ;
575 case 169969 : return 131 ;
576 case 169965 : return 130 ;
577 case 169961 : return 129 ;
578 case 169956 : return 128 ;
579 case 169926 : return 127 ;
580 case 169924 : return 126 ;
581 case 169923 : return 125 ;
582 case 169922 : return 124 ;
583 case 169919 : return 123 ;
584 case 169918 : return 122 ;
585 case 169914 : return 121 ;
586 case 169859 : return 120 ;
587 case 169858 : return 119 ;
588 case 169855 : return 118 ;
589 case 169846 : return 117 ;
590 case 169838 : return 116 ;
591 case 169837 : return 115 ;
592 case 169835 : return 114 ;
593 case 169683 : return 113 ;
594 case 169628 : return 112 ;
595 case 169591 : return 111 ;
596 case 169590 : return 110 ;
597 case 169588 : return 109 ;
598 case 169587 : return 108 ;
599 case 169586 : return 107 ;
600 case 169584 : return 106 ;
601 case 169557 : return 105 ;
602 case 169555 : return 104 ;
603 case 169554 : return 103 ;
604 case 169553 : return 102 ;
605 case 169550 : return 101 ;
606 case 169515 : return 100 ;
607 case 169512 : return 99 ;
608 case 169506 : return 98 ;
609 case 169504 : return 97 ;
610 case 169498 : return 96 ;
611 case 169475 : return 95 ;
612 case 169420 : return 94 ;
613 case 169419 : return 93 ;
614 case 169418 : return 92 ;
615 case 169417 : return 91 ;
616 case 169415 : return 90 ;
617 case 169411 : return 89 ;
618 case 169238 : return 88 ;
619 case 169236 : return 87 ;
620 case 169167 : return 86 ;
621 case 169160 : return 85 ;
622 case 169156 : return 84 ;
623 case 169148 : return 83 ;
624 case 169145 : return 82 ;
625 case 169144 : return 81 ;
626 case 169143 : return 80 ;
627 case 169138 : return 79 ;
628 case 169099 : return 78 ;
629 case 169094 : return 77 ;
630 case 169091 : return 76 ;
631 case 169045 : return 75 ;
632 case 169044 : return 74 ;
633 case 169040 : return 73 ;
634 case 169035 : return 72 ;
635 case 168992 : return 71 ;
636 case 168988 : return 70 ;
637 case 168984 : return 69 ;
638 case 168826 : return 68 ;
639 case 168777 : return 67 ;
640 case 168514 : return 66 ;
641 case 168512 : return 65 ;
642 case 168511 : return 64 ;
643 case 168467 : return 63 ;
644 case 168464 : return 62 ;
645 case 168461 : return 61 ;
646 case 168460 : return 60 ;
647 case 168458 : return 59 ;
648 case 168362 : return 58 ;
649 case 168361 : return 57 ;
650 case 168356 : return 56 ;
651 case 168342 : return 55 ;
652 case 168341 : return 54 ;
653 case 168325 : return 53 ;
654 case 168322 : return 52 ;
655 case 168318 : return 51 ;
656 case 168311 : return 50 ;
657 case 168310 : return 49 ;
658 case 168213 : return 48 ;
659 case 168212 : return 47 ;
660 case 168208 : return 46 ;
661 case 168207 : return 45 ;
662 case 168206 : return 44 ;
663 case 168205 : return 43 ;
664 case 168204 : return 42 ;
665 case 168203 : return 41 ;
666 case 168181 : return 40 ;
667 case 168177 : return 39 ;
668 case 168175 : return 38 ;
669 case 168173 : return 37 ;
670 case 168172 : return 36 ;
671 case 168171 : return 35 ;
672 case 168115 : return 34 ;
673 case 168108 : return 33 ;
674 case 168107 : return 32 ;
675 case 168105 : return 31 ;
676 case 168104 : return 30 ;
677 case 168103 : return 29 ;
678 case 168076 : return 28 ;
679 case 168069 : return 27 ;
680 case 168068 : return 26 ;
681 case 168066 : return 25 ;
682 case 167988 : return 24 ;
683 case 167987 : return 23 ;
684 case 167986 : return 22 ;
685 case 167985 : return 21 ;
686 case 167921 : return 20 ;
687 case 167920 : return 19 ;
688 case 167915 : return 18 ;
689 case 167909 : return 17 ;
690 case 167903 : return 16 ;
691 case 167902 : return 15 ;
692 case 167818 : return 14 ;
693 case 167814 : return 13 ;
694 case 167813 : return 12 ;
695 case 167808 : return 11 ;
696 case 167807 : return 10 ;
697 case 167806 : return 9 ;
698 case 167713 : return 8 ;
699 case 167712 : return 7 ;
700 case 167711 : return 6 ;
701 case 167706 : return 5 ;
702 case 167693 : return 4 ;
703 case 166532 : return 3 ;
704 case 166530 : return 2 ;
705 case 166529 : return 1 ;
706
707 default : return 199;
708 }
709 }
710 if(fPeriod== kLHC10h){
711 switch(run){
712 case 139517 : return 137;
713 case 139514 : return 136;
714 case 139513 : return 135;
715 case 139511 : return 134;
716 case 139510 : return 133;
717 case 139507 : return 132;
718 case 139505 : return 131;
719 case 139504 : return 130;
720 case 139503 : return 129;
721 case 139470 : return 128;
722 case 139467 : return 127;
723 case 139466 : return 126;
724 case 139465 : return 125;
725 case 139440 : return 124;
726 case 139439 : return 123;
727 case 139438 : return 122;
728 case 139437 : return 121;
729 case 139360 : return 120;
730 case 139329 : return 119;
731 case 139328 : return 118;
732 case 139314 : return 117;
733 case 139311 : return 116;
734 case 139310 : return 115;
735 case 139309 : return 114;
736 case 139308 : return 113;
737 case 139173 : return 112;
738 case 139172 : return 111;
739 case 139110 : return 110;
740 case 139107 : return 109;
741 case 139105 : return 108;
742 case 139104 : return 107;
743 case 139042 : return 106;
744 case 139038 : return 105;
745 case 139037 : return 104;
746 case 139036 : return 103;
747 case 139029 : return 102;
748 case 139028 : return 101;
749 case 138983 : return 100;
750 case 138982 : return 99;
751 case 138980 : return 98;
752 case 138979 : return 97;
753 case 138978 : return 96;
754 case 138977 : return 95;
755 case 138976 : return 94;
756 case 138973 : return 93;
757 case 138972 : return 92;
758 case 138965 : return 91;
759 case 138924 : return 90;
760 case 138872 : return 89;
761 case 138871 : return 88;
762 case 138870 : return 87;
763 case 138837 : return 86;
764 case 138830 : return 85;
765 case 138828 : return 84;
766 case 138826 : return 83;
767 case 138796 : return 82;
768 case 138795 : return 81;
769 case 138742 : return 80;
770 case 138732 : return 79;
771 case 138730 : return 78;
772 case 138666 : return 77;
773 case 138662 : return 76;
774 case 138653 : return 75;
775 case 138652 : return 74;
776 case 138638 : return 73;
777 case 138624 : return 72;
778 case 138621 : return 71;
779 case 138583 : return 70;
780 case 138582 : return 69;
781 case 138579 : return 68;
782 case 138578 : return 67;
783 case 138534 : return 66;
784 case 138469 : return 65;
785 case 138442 : return 64;
786 case 138439 : return 63;
787 case 138438 : return 62;
788 case 138396 : return 61;
789 case 138364 : return 60;
790 case 138359 : return 59;
791 case 138275 : return 58;
792 case 138225 : return 57;
793 case 138201 : return 56;
794 case 138200 : return 55;
795 case 138197 : return 54;
796 case 138192 : return 53;
797 case 138190 : return 52;
798 case 138154 : return 51;
799 case 138153 : return 50;
800 case 138151 : return 49;
801 case 138150 : return 48;
802 case 138126 : return 47;
803 case 138125 : return 46;
804 case 137848 : return 45;
805 case 137847 : return 44;
806 case 137844 : return 43;
807 case 137843 : return 42;
808 case 137752 : return 41;
809 case 137751 : return 40;
810 case 137748 : return 39;
811 case 137724 : return 38;
812 case 137722 : return 37;
813 case 137718 : return 36;
814 case 137704 : return 35;
815 case 137693 : return 34;
816 case 137692 : return 33;
817 case 137691 : return 32;
818 case 137689 : return 31;
819 case 137686 : return 30;
820 case 137685 : return 29;
821 case 137639 : return 28;
822 case 137638 : return 27;
823 case 137608 : return 26;
824 case 137595 : return 25;
825 case 137549 : return 24;
826 case 137546 : return 23;
827 case 137544 : return 22;
828 case 137541 : return 21;
829 case 137539 : return 20;
830 case 137531 : return 19;
831 case 137530 : return 18;
832 case 137443 : return 17;
833 case 137441 : return 16;
834 case 137440 : return 15;
835 case 137439 : return 14;
836 case 137434 : return 13;
837 case 137432 : return 12;
838 case 137431 : return 11;
839 case 137430 : return 10;
840 case 137366 : return 9;
841 case 137243 : return 8;
842 case 137236 : return 7;
843 case 137235 : return 6;
844 case 137232 : return 5;
845 case 137231 : return 4;
846 case 137165 : return 3;
847 case 137162 : return 2;
848 case 137161 : return 1;
849 default : return 199;
850 }
851 }
852 if( kLHC13 == fPeriod )
853 {
854 switch(run)
855 {
856 case 195344 : return 1;
857 case 195346 : return 2;
858 case 195351 : return 3;
859 case 195389 : return 4;
860 case 195390 : return 5;
861 case 195391 : return 6;
862 case 195478 : return 7;
863 case 195479 : return 8;
864 case 195480 : return 9;
865 case 195481 : return 10;
866 case 195482 : return 11;
867 case 195483 : return 12;
868 case 195529 : return 13;
869 case 195531 : return 14;
870 case 195532 : return 15;
871 case 195566 : return 16;
872 case 195567 : return 17;
873 case 195568 : return 18;
874 case 195592 : return 19;
875 case 195593 : return 20;
876 case 195596 : return 21;
877 case 195633 : return 22;
878 case 195635 : return 23;
879 case 195644 : return 24;
880 case 195673 : return 25;
881 case 195675 : return 26;
882 case 195676 : return 27;
883 case 195677 : return 28;
884 case 195681 : return 29;
885 case 195682 : return 30;
886 case 195720 : return 31;
887 case 195721 : return 32;
888 case 195722 : return 33;
889 case 195724 : return 34;
890 case 195725 : return 34;
891 case 195726 : return 35;
892 case 195727 : return 36;
893 case 195760 : return 37;
894 case 195761 : return 38;
895 case 195765 : return 39;
896 case 195767 : return 40;
897 case 195783 : return 41;
898 case 195787 : return 42;
899 case 195826 : return 43;
900 case 195827 : return 44;
901 case 195829 : return 45;
902 case 195830 : return 46;
903 case 195831 : return 47;
904 case 195867 : return 48;
905 case 195869 : return 49;
906 case 195871 : return 50;
907 case 195872 : return 51;
908 case 195873 : return 52;
909 case 195935 : return 53;
910 case 195949 : return 54;
911 case 195950 : return 55;
912 case 195954 : return 56;
913 case 195955 : return 57;
914 case 195958 : return 58;
915 case 195989 : return 59;
916 case 195994 : return 60;
917 case 195998 : return 61;
918 case 196000 : return 62;
919 case 196006 : return 63;
920 case 196085 : return 64;
921 case 196089 : return 65;
922 case 196090 : return 66;
923 case 196091 : return 67;
924 case 196099 : return 68;
925 case 196105 : return 69;
926 case 196107 : return 70;
927 case 196185 : return 71;
928 case 196187 : return 72;
929 case 196194 : return 73;
930 case 196197 : return 74;
931 case 196199 : return 75;
932 case 196200 : return 76;
933 case 196201 : return 77;
934 case 196203 : return 78;
935 case 196208 : return 79;
936 case 196214 : return 80;
937 case 196308 : return 81;
938 case 196309 : return 82;
939 case 196310 : return 83;
940 case 196311 : return 84;
941 case 196433 : return 85;
942 case 196474 : return 86;
943 case 196475 : return 87;
944 case 196477 : return 88;
945 case 196528 : return 89;
946 case 196533 : return 90;
947 case 196535 : return 91;
948 case 196563 : return 92;
949 case 196564 : return 93;
950 case 196566 : return 94;
951 case 196568 : return 95;
952 case 196601 : return 96;
953 case 196605 : return 97;
954 case 196608 : return 98;
955 case 196646 : return 99;
956 case 196648 : return 100;
957 case 196701 : return 101;
958 case 196702 : return 102;
959 case 196703 : return 103;
960 case 196706 : return 104;
961 case 196714 : return 105;
962 case 196720 : return 106;
963 case 196721 : return 107;
964 case 196722 : return 108;
965 case 196772 : return 109;
966 case 196773 : return 110;
967 case 196774 : return 111;
968 case 196869 : return 112;
969 case 196870 : return 113;
970 case 196874 : return 114;
971 case 196876 : return 115;
972 case 196965 : return 116;
973 case 196967 : return 117;
974 case 196972 : return 118;
975 case 196973 : return 119;
976 case 196974 : return 120;
977 case 197003 : return 121;
978 case 197011 : return 122;
979 case 197012 : return 123;
980 case 197015 : return 124;
981 case 197027 : return 125;
982 case 197031 : return 126;
983 case 197089 : return 127;
984 case 197090 : return 128;
985 case 197091 : return 129;
986 case 197092 : return 130;
987 case 197094 : return 131;
988 case 197098 : return 132;
989 case 197099 : return 133;
990 case 197138 : return 134;
991 case 197139 : return 135;
992 case 197142 : return 136;
993 case 197143 : return 137;
994 case 197144 : return 138;
995 case 197145 : return 139;
996 case 197146 : return 140;
997 case 197147 : return 140;
998 case 197148 : return 141;
999 case 197149 : return 142;
1000 case 197150 : return 143;
1001 case 197152 : return 144;
1002 case 197153 : return 145;
1003 case 197184 : return 146;
1004 case 197189 : return 147;
1005 case 197247 : return 148;
1006 case 197248 : return 149;
1007 case 197254 : return 150;
1008 case 197255 : return 151;
1009 case 197256 : return 152;
1010 case 197258 : return 153;
1011 case 197260 : return 154;
1012 case 197296 : return 155;
1013 case 197297 : return 156;
1014 case 197298 : return 157;
1015 case 197299 : return 158;
1016 case 197300 : return 159;
1017 case 197302 : return 160;
1018 case 197341 : return 161;
1019 case 197342 : return 162;
1020 case 197348 : return 163;
1021 case 197349 : return 164;
1022 case 197351 : return 165;
1023 case 197386 : return 166;
1024 case 197387 : return 167;
1025 case 197388 : return 168;
1026 default : return 199;
1027 }
1028 }
1029 if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) )
1030 {
1031 AliWarning("Period not defined");
1032 }
1033 return 1;
1034 }
1035
1036//_______________________________________________________________________________
1037Bool_t AliPHOSCorrelations::RejectTriggerMaskSelection()
1038{
1039 const Bool_t REJECT = true;
1040 const Bool_t ACCEPT = false;
1041
1042 // No need to check trigger mask if no selection is done
1043 if( kNoSelection == fInternalTriggerSelection )
1044 return ACCEPT;
1045
1046 Bool_t reject = REJECT;
1047
1048 Bool_t isMB = (fEvent->GetTriggerMask() & (ULong64_t(1)<<1));
1049 Bool_t isCentral = (fEvent->GetTriggerMask() & (ULong64_t(1)<<4));
1050 Bool_t isSemiCentral = (fEvent->GetTriggerMask() & (ULong64_t(1)<<7));
1051
1052
1053 if( kCentralInclusive == fInternalTriggerSelection
1054 && isCentral ) reject = ACCEPT; // accept event.
1055 else if( kCentralExclusive == fInternalTriggerSelection
1056 && isCentral && !isSemiCentral && !isMB ) reject = ACCEPT; // accept event.
1057
1058 else if( kSemiCentralInclusive == fInternalTriggerSelection
1059 && isSemiCentral ) reject = ACCEPT; // accept event
1060 else if( kSemiCentralExclusive == fInternalTriggerSelection
1061 && isSemiCentral && !isCentral && !isMB ) reject = ACCEPT; // accept event.
1062
1063 else if( kMBInclusive == fInternalTriggerSelection
1064 && isMB ) reject = ACCEPT; // accept event.
1065 else if( kMBExclusive == fInternalTriggerSelection
1066 && isMB && !isCentral && !isSemiCentral ) reject = ACCEPT; // accept event.
1067
1068 if( REJECT == reject )
1069 return REJECT;
1070 else {
1071 LogSelection(kInternalTriggerMaskSelection, fInternalRunNumber);
1072 return ACCEPT;
1073 }
1074}
1075//_______________________________________________________________________________
1076void AliPHOSCorrelations::SetVertex()
1077{
1078 const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
1079 if( primaryVertex )
1080 {
1081 fVertex[0] = primaryVertex->GetX();
1082 fVertex[1] = primaryVertex->GetY();
1083 fVertex[2] = primaryVertex->GetZ();
1084 }
1085 else
1086 {
1087// AliError("Event has 0x0 Primary Vertex, defaulting to origo");
1088 fVertex[0] = 0;
1089 fVertex[1] = 0;
1090 fVertex[2] = 0;
1091 }
1092 fVertexVector = TVector3(fVertex);
1093
1094 fVtxBin=0 ;// No support for vtx binning implemented.
1095}
1096//_______________________________________________________________________________
1097Bool_t AliPHOSCorrelations::RejectEventVertex()
1098{
1099 if( ! fEvent->GetPrimaryVertex() )
1100 return true; // reject
1101 LogSelection(kHasVertex, fInternalRunNumber);
1102
1103 if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ )
1104 return true; // reject
1105 LogSelection(kHasAbsVertex, fInternalRunNumber);
1106
1107 return false; // accept event.
1108}
1109//_______________________________________________________________________________
1110void AliPHOSCorrelations::SetCentrality()
1111{
1112 AliCentrality *centrality = fEvent->GetCentrality();
1113 if( centrality )
1114 fCentrality=centrality->GetCentralityPercentile(fCentralityEstimator);
1115 else
1116 {
1117 AliError("Event has 0x0 centrality");
1118 fCentrality = -1.;
1119 }
1120
1121 //cout<<"fCentrality: "<<fCentrality<<endl;
1122 //FillHistogram("hCentrality",fCentrality,fInternalRunNumber-0.5) ;
1123 fCentBin = GetCentralityBin(fCentrality);
1124}
1125//_______________________________________________________________________________
1126Bool_t AliPHOSCorrelations::RejectEventCentrality()
1127{
1128 if (fCentrality<fCentCutoffDown)
1129 return true; //reject
1130 if(fCentrality>fCentCutoffUp)
1131 return true;
1132
1133 return false; // accept event.
1134}
1135//_______________________________________________________________________________
1136void AliPHOSCorrelations::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed){
1137// Define centrality bins by their edges
1138 for(int i=0; i<edges.GetSize()-1; ++i)
1139 if(edges.At(i) > edges.At(i+1)) AliFatal("edges are not sorted");
1140 if( edges.GetSize() != nMixed.GetSize()+1) AliFatal("edges and nMixed don't have appropriate relative sizes");
1141
1142 fCentEdges = edges;
1143 fCentNMixed = nMixed;
1144}
1145//_______________________________________________________________________________
1146Int_t AliPHOSCorrelations::GetCentralityBin(Float_t centralityV0M){
1147 int lastBinUpperIndex = fCentEdges.GetSize() -1;
1148 if( centralityV0M > fCentEdges[lastBinUpperIndex] ) {
1149 if( fDebug >= 1 )
1150 AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1151 return lastBinUpperIndex-1;
1152 }
1153 if( centralityV0M < fCentEdges[0] ) {
1154 if( fDebug >= 1 )
1155 AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1156 return 0;
1157 }
1158
1159 fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1160 return fCentBin;
1161}
1162//_______________________________________________________________________________
1163void AliPHOSCorrelations::SetCentralityBorders (double down , double up ){
1164 if (down < 0. || up > 100 || up<=down)
1165 AliError( Form("Warning. Bad value of centrality borders. Setting as default: fCentCutoffDown=%2.f, fCentCutoffUp=%2.f",fCentCutoffDown,fCentCutoffUp) );
1166 else{
1167 fCentCutoffDown = down;
1168 fCentCutoffUp = up;
1169 AliInfo( Form("Centrality border was set as fCentCutoffDown=%2.f, fCentCutoffUp=%2.f",fCentCutoffDown,fCentCutoffUp) );
1170 }
1171}
1172
1173//_______________________________________________________________________________
1174void AliPHOSCorrelations::EvalReactionPlane()
1175{
1176 // assigns: fHaveTPCRP and fRP
1177 // also does a few histogram fills
1178
1179 AliEventplane *eventPlane = fEvent->GetEventplane();
1180 if( ! eventPlane ) { AliError("Event has no event plane"); return; }
1181
1182 Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
1183
1184 if(reactionPlaneQ>=999 || reactionPlaneQ < 0.)
1185 {
1186 //reaction plain was not defined
1187 fHaveTPCRP = kFALSE;
1188 }
1189 else
1190 {
1191 fHaveTPCRP = kTRUE;
1192 }
1193
1194 if(fHaveTPCRP)
1195 fRP = reactionPlaneQ;
1196 else
1197 fRP = 0.;
1198
1199 FillHistogram("phiRPflat",fRP,fCentrality) ;
1200}
1201//_______________________________________________________________________________
1202Int_t AliPHOSCorrelations::GetRPBin()
1203{
1204 Double_t averageRP;
1205 averageRP = fRP ; // If possible, it is better to have EP bin from TPC
1206 // to have similar events for miximng (including jets etc) (fRPV0A+fRPV0C+fRP) /3.;
1207
1208 fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1209
1210 if( fEMRPBin > (Int_t)fNEMRPBins-1 )
1211 fEMRPBin = fNEMRPBins-1 ;
1212 else
1213 if(fEMRPBin < 0) fEMRPBin=0;
1214
1215 return fEMRPBin;
1216}
1217//_______________________________________________________________________________
1218void AliPHOSCorrelations::SelectPhotonClusters()
1219{
1220 //Selects PHOS clusters
1221
1222 // clear (or create) array for holding events photons/clusters
1223 if(fCaloPhotonsPHOS)
1224 fCaloPhotonsPHOS->Clear();
1225 else{
1226 fCaloPhotonsPHOS = new TClonesArray("AliCaloPhoton",200);
1227 }
1228
1229 Int_t nclu = fEvent->GetNumberOfCaloClusters() ;
1230 Int_t inPHOS=0 ;
1231 for (Int_t i=0; i<nclu; i++) {
1232 AliVCluster *clu = fEvent->GetCaloCluster(i);
1233 if ( !clu->IsPHOS() ) continue ;
1234 if( clu->E()< fMinClusterEnergy) continue; // reject cluster
1235
1236
1237 Double_t distBC=clu->GetDistanceToBadChannel();
1238 if(distBC<fMinBCDistance)
1239 continue ;
1240
1241 if(clu->GetNCells() < fMinNCells) continue ;
1242 if(clu->GetM02() < fMinM02) continue ;
1243
1244 if(fTOFCutEnabled){
1245 Double_t tof = clu->GetTOF();
1246 if(TMath::Abs(tof) > fTOFCut )
1247 continue ;
1248 }
1249 TLorentzVector lorentzMomentum;
1250 Double_t ecore = clu->GetMCEnergyFraction();
1251
1252 clu->GetMomentum(lorentzMomentum, fVertex);
1253 lorentzMomentum*=ecore/lorentzMomentum.E() ;
1254
1255 if(inPHOS>=fCaloPhotonsPHOS->GetSize()){
1256 fCaloPhotonsPHOS->Expand(inPHOS+50) ;
1257 }
1258
1259 AliCaloPhoton * ph =new((*fCaloPhotonsPHOS)[inPHOS]) AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E());
1260 inPHOS++ ;
1261 ph->SetCluster(clu);
1262
1263 Float_t cellId=clu->GetCellAbsId(0) ;
1264 Int_t mod = (Int_t)TMath:: Ceil(cellId/(56*64) ) ;
1265 ph->SetModule(mod) ;
1266
1267 ph->SetNCells(clu->GetNCells());
1268 ph->SetDispBit(clu->GetDispersion()<2.5) ;
1269 ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
1270 }
1271}
1272//_______________________________________________________________________________
1273void AliPHOSCorrelations::SelectAccosiatedTracks()
1274{
1275 // clear (or create) array for holding events photons/clusters
1276 if(fTracksTPC)
1277 fTracksTPC->Clear();
1278 else
1279 {
1280 fTracksTPC = new TClonesArray("TLorentzVector",12000);
1281 }
1282 Int_t iTracks=0 ;
1283 for (Int_t i=0; i<fEvent->GetNumberOfTracks(); i++)
1284 {
1285
1286 AliVParticle *track = fEvent->GetTrack(i);
1287 if(fEventESD){
1288 if(!SelectESDTrack((AliESDtrack*)track)) continue ;
1289 }
1290 else{
1291 if(!SelectAODTrack((AliAODTrack*)track)) continue ;
1292 }
1293 Double_t px = track->Px();
1294 Double_t py = track->Py();
1295 Double_t pz = track->Pz() ;
1296 Double_t e = track->E() ;
1297
1298 if(iTracks>=fTracksTPC->GetSize())
1299 fTracksTPC->Expand(iTracks+50) ;
1300
1301 new((*fTracksTPC)[iTracks]) TLorentzVector(px, py, pz,e);
1302 iTracks++ ;
1303 }
1304}
1305//_______________________________________________________________________________
1306void AliPHOSCorrelations::ConsiderPi0s()
1307{
1308
1309 const Int_t nPHOS=fCaloPhotonsPHOS->GetEntriesFast() ;
1310 for(Int_t i1=0; i1 < nPHOS-1; i1++){
1311 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1312 for (Int_t i2=i1+1; i2<nPHOS; i2++){
1313 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
1314 TLorentzVector p12 = *ph1 + *ph2;
1315
1316 Double_t phiTrigger=p12.Phi() ;
1317 Double_t etaTrigger=p12.Eta() ;
1318
1319 Double_t m=p12.M() ;
1320 Double_t pt=p12.Pt() ;
1321 int mod1 = ph1->Module() ;
1322 int mod2 = ph2->Module() ;
1323
1324// if(!TestMass(m,pt)) continue;
1325
1326 FillHistogram("clu_phieta",phiTrigger,etaTrigger);
1327 FillHistogram("clusingle_phieta",ph1->Phi(), ph1->Eta());
1328 FillHistogram("clusingle_phieta",ph2->Phi(), ph2->Eta());
1329
1330 FillHistogram("all_mpt",m, pt);
1331
1332 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1333 FillHistogram("cpv_mpt",m, pt);
1334
1335 if ( ph1->IsDispOK() && ph2->IsDispOK() ){
1336 FillHistogram("disp_mpt",m, pt);
1337 if ( ph1->IsCPVOK() && ph2->IsCPVOK() ) {
1338 FillHistogram("both_mpt",m, pt);
1339 if(mod1 == mod2) FillHistogram(Form("both%d_mpt",mod1),m, pt);
1340 }
1341 }
1342
1343 if(!TestMass(m,pt)) continue;
1344 // Take track's angles and compare with cluster's angles.
1345 for(Int_t i3=0; i3<fTracksTPC->GetEntriesFast(); i3++){
1346 TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i3);
1347
1348 Double_t phiAssoc = track->Phi();
1349 Double_t etaAssoc = track->Eta();
1350 Double_t ptAssoc = track->Pt();
1351
1352 Double_t dPhi = phiAssoc - phiTrigger;
1353 while (dPhi > 1.5*TMath::Pi()) dPhi-=2*TMath::Pi();
1354 while (dPhi < -.5*TMath::Pi()) dPhi+=2*TMath::Pi();
1355
1356 Double_t dEta = etaAssoc - etaTrigger;
1357
1358 Double_t ptAssocBin=GetAssocBin(ptAssoc) ;
1359 FillHistogram(Form("all_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1360 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1361 FillHistogram(Form("cpv_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1362
1363 if ( ph1->IsDispOK() && ph2->IsDispOK() ){
1364 FillHistogram(Form("disp_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1365 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1366 FillHistogram(Form("both_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1367 }
1368 }
1369 }
1370 }
1371}
1372
1373//_______________________________________________________________________________
1374void AliPHOSCorrelations::ConsiderPi0sMix()
1375{
1376 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1377 for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1378 TClonesArray * mixPHOS = static_cast<TClonesArray*>(arrayList->At(evi));
1379 for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast(); i1++){
1380 AliCaloPhoton * ph1 = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1381 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast(); i2++){
1382 AliCaloPhoton * ph2 = (AliCaloPhoton*)mixPHOS->At(i2) ;
1383 TLorentzVector p12 = *ph1 + *ph2;
1384 Double_t m=p12.M() ;
1385 Double_t pt=p12.Pt() ;
1386 int mod1 = ph1->Module() ;
1387 int mod2 = ph2->Module() ;
1388
1389// if ( m<fMassInvMin || m>fMassInvMax ) continue;
1390 FillHistogram("mix_all_mpt", m, pt);
1391 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1392 FillHistogram("mix_cpv_mpt",m, pt);
1393 if ( ph1->IsDispOK() && ph2->IsDispOK() ){
1394 FillHistogram("mix_disp_mpt",m, pt);
1395 if ( ph1->IsCPVOK() && ph2->IsCPVOK() ){
1396 FillHistogram("mix_both_mpt",m, pt);
1397 if (mod1 == mod2) FillHistogram(Form("mix_both%d_mpt",mod1),m, pt);
1398 }
1399 }
1400 }
1401 }
1402 }
1403}
1404//_______________________________________________________________________________
1405void AliPHOSCorrelations::ConsiderTracksMix()
1406{
1407 TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
1408 for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
1409 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1410 for (Int_t i2=0; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++){
1411 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
1412 TLorentzVector p12 = *ph1 + *ph2;
1413 Double_t phiTrigger=p12.Phi() ;
1414 Double_t etaTrigger=p12.Eta() ;
1415
1416 Double_t m=p12.M() ;
1417 Double_t pt=p12.Pt() ;
1418
1419 if(!TestMass(m,pt)) continue;
1420 for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1421 TClonesArray * mixTracks = static_cast<TClonesArray*>(arrayList->At(evi));
1422 for(Int_t i3=0; i3<mixTracks->GetEntriesFast(); i3++){
1423 TLorentzVector * track = (TLorentzVector*)mixTracks->At(i3);
1424
1425 Double_t phiAssoc = track->Phi();
1426 Double_t etaAssoc = track->Eta();
1427 Double_t ptAssoc = track->Pt();
1428
1429 Double_t ptAssocBin=GetAssocBin(ptAssoc) ;
1430
1431 Double_t dPhi = phiAssoc - phiTrigger;
1432 while (dPhi > 1.5*TMath::Pi()) dPhi-=2*TMath::Pi();
1433 while (dPhi < -.5*TMath::Pi()) dPhi+=2*TMath::Pi();
1434
1435 Double_t dEta = etaAssoc - etaTrigger;
1436
1437 FillHistogram(Form("mix_all_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1438 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1439 FillHistogram(Form("mix_cpv_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1440
1441 if ( ph1->IsDispOK() && ph2->IsDispOK() ){
1442 FillHistogram(Form("mix_disp_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1443 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1444 FillHistogram(Form("mix_both_ptphieta_ptAssoc_%3.1f",ptAssocBin), pt, dPhi, dEta);
1445 }
1446 }
1447 }
1448 }
1449 }
1450}
1451//_______________________________________________________________________________
1452TList* AliPHOSCorrelations::GetCaloPhotonsPHOSList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin){
1453
1454 int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1455 if( fCaloPhotonsPHOSLists->At(offset) ) {
1456 TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1457 return list;
1458 }
1459 else{ // no list for this bin has been created, yet
1460 TList* list = new TList();
1461 fCaloPhotonsPHOSLists->AddAt(list, offset);
1462 return list;
1463 }
1464}
1465//_______________________________________________________________________________
1466TList* AliPHOSCorrelations::GetTracksTPCList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin){
1467
1468 int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1469 if( fTracksTPCLists->At(offset) ) { // list exists
1470 TList* list = dynamic_cast<TList*> (fTracksTPCLists->At(offset));
1471 return list;
1472 }
1473 else { // no list for this bin has been created, yet
1474 TList* list = new TList();
1475 fTracksTPCLists->AddAt(list, offset);
1476 return list;
1477 }
1478}
1479//_______________________________________________________________________________
1480Double_t AliPHOSCorrelations::GetAssocBin(Double_t pt){
1481 //Calculates bin
1482 for(Int_t i=1; i<fAssocBins.GetSize(); i++){
1483 if(pt>fAssocBins.At(i-1) && pt<fAssocBins.At(i))
1484 return fAssocBins.At(i) ;
1485 }
1486 return fAssocBins.At(fAssocBins.GetSize()-1) ;
1487}
1488//_______________________________________________________________________________
1489void AliPHOSCorrelations::FillTrackEtaPhi()
1490// Distribution TPC's tracks by angles.
1491{
1492 for (Int_t i1=0; i1<fTracksTPC->GetEntriesFast(); i1++){
1493 TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i1);
1494 FillHistogram( "track_phieta", track->Phi(), track->Eta() );
1495 }
1496}
1497
1498//_______________________________________________________________________________
1499void AliPHOSCorrelations::UpdatePhotonLists()
1500{
1501 //Now we either add current events to stack or remove
1502 //If no photons in current event - no need to add it to mixed
1503
1504 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1505 if( fDebug >= 2 )
1506 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1507 if(fCaloPhotonsPHOS->GetEntriesFast()>0)
1508 {
1509 arrayList->AddFirst(fCaloPhotonsPHOS) ;
1510 fCaloPhotonsPHOS=0x0;
1511 if(arrayList->GetEntries() > fCentNMixed[fCentBin])
1512 { // Remove redundant events
1513 TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
1514 arrayList->RemoveLast() ;
1515 delete tmp;
1516 }
1517 }
1518}
1519//_______________________________________________________________________________
1520void AliPHOSCorrelations::UpdateTrackLists()
1521{
1522 //Now we either add current events to stack or remove
1523 //If no photons in current event - no need to add it to mixed
1524
1525 TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
1526
1527 if( fDebug >= 2 )
1528 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1529 if(fTracksTPC->GetEntriesFast()>0)
1530 {
1531
1532 arrayList->AddFirst(fTracksTPC) ;
1533 fTracksTPC=0x0;
1534 if(arrayList->GetEntries() > fCentNMixed[fCentBin])
1535 { // Remove redundant events
1536 TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
1537 arrayList->RemoveLast() ;
1538 delete tmp;
1539 }
1540 }
1541}
1542//_______________________________________________________________________________
1543Bool_t AliPHOSCorrelations::SelectESDTrack(AliESDtrack * t) const
1544// Estimate if this track can be used for the RP calculation. If all right - return "TRUE"
1545{
1546 Float_t pt=t->Pt();
1547 if(pt<0.5 || pt>10.) return kFALSE ;
1548 if(fabs( t->Eta() )>0.8) return kFALSE;
1549 if(!fESDtrackCuts->AcceptTrack(t)) return kFALSE ;
1550 return kTRUE ;
1551}
1552//_______________________________________________________________________________
1553Bool_t AliPHOSCorrelations::SelectAODTrack(AliAODTrack * t) const
1554// Estimate if this track can be used for the RP calculation. If all right - return "TRUE"
1555{
1556 Float_t pt=t->Pt();
1557 if(pt<0.5 || pt>10.) return kFALSE ;
1558 if(fabs( t->Eta() )>0.8) return kFALSE;
1559 if(fCheckHibridGlobal == kOnlyHibridTracks)
1560 {
1561 if(!t->IsHybridGlobalConstrainedGlobal())
1562 return kFALSE ;
1563 }
1564
1565 if (fCheckHibridGlobal == kWithOutHibridTracks)
1566 {
1567 if(t->IsHybridGlobalConstrainedGlobal())
1568 return kFALSE ;
1569 }
1570
1571 return kTRUE ;
1572}
1573//_______________________________________________________________________________
1574void AliPHOSCorrelations::SetPeriod(Period period)
1575{
1576 fPeriod = period;
1577 if( kLHC11h == fPeriod ) {
1578 fCentralityEstimator = "V0M";
1579 cout<<"Set Centraity Estimator as \"V0M\" "<<endl;
1580 }
1581 else
1582 if( kLHC10h == fPeriod ) {
1583 fCentralityEstimator = "V0M";
1584 cout<<"Set Centraity Estimator as \"V0M\" "<<endl;
1585 }
1586 else
1587 if( kLHC13 == fPeriod ) {
1588 fCentralityEstimator = "V0A";
1589 cout<<"Set Centraity Estimator as \"V0A\" "<<endl;
1590 }
1591 else
1592 {
1593 AliWarning("Period not defined. fCentralityEstimator set as default \"V0M\".");
1594 fCentralityEstimator = "V0M";
1595 }
1596}
1597
1598//_______________________________________________________________________________
1599void AliPHOSCorrelations::LogProgress(int step)
1600// Fill "step by step" hist
1601{
1602 //FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1603 FillHistogram("hTotSelEvents", step+0.5);
1604}
1605//_______________________________________________________________________________
1606void AliPHOSCorrelations::LogSelection(int step, int internalRunNumber)
1607{
1608 // the +0.5 is not realy neccisarry, but oh well... -henrik
1609 FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1610 //FillHistogram("hTotSelEvents", step+0.5);
1611}
1612//_______________________________________________________________________________
1613Bool_t AliPHOSCorrelations::TestMass(Double_t m, Double_t /*pt*/){
1614 //Check if mair in pi0 peak window
1615 //To make pT-dependent
1616 return (m>fMassInvMin && m<fMassInvMax) ;
1617
1618}
1619//_______________________________________________________________________________
1620void AliPHOSCorrelations::FillHistogram(const char * key,Double_t x)const
1621{
1622 //FillHistogram
1623 TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
1624 if(hist)
1625 hist->Fill(x) ;
1626 else
1627 AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
1628}
1629//_______________________________________________________________________________
1630void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x, Double_t y) const
1631{
1632 //Fills 2D histograms with key
1633 TObject * obj = fOutputContainer->FindObject(key);
1634
1635 TH2 * th2 = dynamic_cast<TH2*> (obj);
1636 if(th2) {
1637 th2->Fill(x, y) ;
1638 return;
1639 }
1640
1641 AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
1642}
1643//_______________________________________________________________________________
1644void AliPHOSCorrelations::FillHistogram(const char * key,Double_t x, Double_t y, Double_t z) const
1645{
1646 //Fills 3D histograms with key
1647 TObject * obj = fOutputContainer->FindObject(key);
1648
1649 TH3 * th3 = dynamic_cast<TH3*> (obj);
1650 if(th3) {
1651 th3->Fill(x, y, z) ;
1652 return;
1653 }
1654
1655 AliError(Form("can not find histogram (of instance TH3) <%s> ",key)) ;
1656}
1657//_______________________________________________________________________________