]> 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
Track selection bit corrected; vertex z-bin bug fixed
[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
42800acf 16// Analysis task for identifion PHOS cluster from Pi0 and extracting pi0-hadron correlation.
17// Author: Daniil Ponomarenko <Daniil.Ponomarenko@cern.ch>
18// 20-Sept-2014
67ef08bd 19
faecffdf 20#include <Riostream.h>
67ef08bd 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"
67ef08bd 51#include "AliEventplane.h"
52#include "AliOADBContainer.h"
53#include "AliAODEvent.h"
54#include "AliAODCaloCells.h"
55#include "AliAODCaloCluster.h"
56#include "AliCaloPhoton.h"
57#include "AliAODVertex.h"
7deabbe7 58#include "AliInputEventHandler.h"
67ef08bd 59
faecffdf 60using std::cout;
61using std::endl;
62
67ef08bd 63ClassImp(AliPHOSCorrelations)
64
65//_______________________________________________________________________________
66AliPHOSCorrelations::AliPHOSCorrelations()
67:AliAnalysisTaskSE(),
9aed5fcd 68 fPHOSGeo(0x0),
42800acf 69 fOutputContainer(0x0),
70 fEvent(0x0),
71 fEventESD(0x0),
72 fEventAOD(0x0),
73 fEventHandler(0),
74 fCaloPhotonsPHOS(0x0),
75 fTracksTPC(0x0),
76 fCaloPhotonsPHOSLists(0x0),
77 fTracksTPCLists(0x0),
78 fRunNumber(-999),
79 fInternalRunNumber(0),
80 fPeriod(kUndefinedPeriod),
9aed5fcd 81 fPHOSEvent(false),
42800acf 82 fMBEvent(false),
83 fNVtxZBins(10),
08a327bf 84 fVtxEdges(10),
42800acf 85 fCentEdges(10),
86 fCentNMixed(),
08a327bf 87 fNEMRPBins(6),
42800acf 88 fAssocBins(),
67ef08bd 89 fVertexVector(),
42800acf 90 fVtxBin(0),
91 fCentralityEstimator("V0M"),
92 fCentrality(0.),
93 fCentBin(0),
94 fHaveTPCRP(0),
95 fRP(0.),
96 fEMRPBin(0),
97 fMaxAbsVertexZ(10.),
98 fCentralityLowLimit(0.),
99 fCentralityHightLimit(90),
7deabbe7 100 fMinClusterEnergy(0.3),
101 fMinBCDistance(0),
102 fMinNCells(3),
103 fMinM02(0.2),
104 fTOFCutEnabled(1),
105 fTOFCut(100.e-9),
08a327bf 106 fUseMassWindowParametrisation(true),
107 fMassInvMeanMin(0.13),
108 fMassInvMeanMax(0.145),
109 fNSigmaWidth(0.),
110 fUseEfficiency(true),
111 fESDtrackCuts(0x0),
112 fSelectHybridTracks(0),
113 fTrackStatus(0),
114 fTrackFilterMask(786),
115 fSelectSPDHitTracks(0),
116 fSelectFractionTPCSharedClusters(kTRUE),
117 fCutTPCSharedClustersFraction(0.4)
42800acf 118{
119 //Deafult constructor, no memory allocations here
42800acf 120}
121
122//_______________________________________________________________________________
123AliPHOSCorrelations::AliPHOSCorrelations(const char *name)
124:AliAnalysisTaskSE(name),
125 fPHOSGeo(0x0),
126 fOutputContainer(0x0),
7deabbe7 127 fEvent(0x0),
128 fEventESD(0x0),
129 fEventAOD(0x0),
130 fEventHandler(0),
42800acf 131 fCaloPhotonsPHOS(0x0),
132 fTracksTPC(0x0),
133 fCaloPhotonsPHOSLists(0x0),
134 fTracksTPCLists(0x0),
7deabbe7 135 fRunNumber(-999),
136 fInternalRunNumber(0),
42800acf 137 fPeriod(kUndefinedPeriod),
138 fPHOSEvent(false),
139 fMBEvent(false),
140 fNVtxZBins(10),
08a327bf 141 fVtxEdges(10),
42800acf 142 fCentEdges(10),
143 fCentNMixed(),
08a327bf 144 fNEMRPBins(6),
42800acf 145 fAssocBins(),
146 fVertexVector(),
7deabbe7 147 fVtxBin(0),
148 fCentralityEstimator("V0M"),
149 fCentrality(0.),
150 fCentBin(0),
151 fHaveTPCRP(0),
152 fRP(0.),
153 fEMRPBin(0),
42800acf 154 fMaxAbsVertexZ(10.),
155 fCentralityLowLimit(0.),
156 fCentralityHightLimit(90),
42800acf 157 fMinClusterEnergy(0.3),
158 fMinBCDistance(0),
159 fMinNCells(3),
160 fMinM02(0.2),
161 fTOFCutEnabled(1),
162 fTOFCut(100.e-9),
08a327bf 163 fUseMassWindowParametrisation(true),
164 fMassInvMeanMin(0.13),
165 fMassInvMeanMax(0.145),
166 fNSigmaWidth(0.),
167 fUseEfficiency(true),
168 fESDtrackCuts(0x0),
169 fSelectHybridTracks(0),
170 fTrackStatus(0),
171 fTrackFilterMask(786),
172 fSelectSPDHitTracks(0),
173 fSelectFractionTPCSharedClusters(kTRUE),
174 fCutTPCSharedClustersFraction(0.4)
7deabbe7 175{
176 // Constructor
177 // Output slots #0 write into a TH1 container
178 DefineOutput(1,THashList::Class());
179
42800acf 180 fMassMean[0] = 1.00796e-05 ;
181 fMassMean[1] = 0.136096 ;
08a327bf 182 fMassSigma[0] = 0.00383029 ;
183 fMassSigma[1] = 0.0041709 ;
184 fMassSigma[2] = 0.00468736 ;
42800acf 185
186 const Int_t nPtAssoc = 10 ;
187 Double_t ptAssocBins[nPtAssoc] = {0.,0.5,1.0,1.5,2.0,3.,5.,7.,10.,16} ;
7deabbe7 188 fAssocBins.Set(nPtAssoc,ptAssocBins) ;
189
190 const int nbins = 9;
191 Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
42800acf 192 TArrayD centEdges( nbins+1, edges );
7deabbe7 193 Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
7deabbe7 194 TArrayI centNMixed(nbins, nMixed);
195 SetCentralityBinning(centEdges, centNMixed);
08a327bf 196 SetVertexBinning();
7deabbe7 197
42800acf 198 fVertex[0] = 0;
199 fVertex[1] = 0;
200 fVertex[2] = 0;
7deabbe7 201
202 SetGeometry();
203
204 ZeroingVariables();
205}
206
67ef08bd 207//_______________________________________________________________________________
208AliPHOSCorrelations::AliPHOSCorrelations(const char *name, Period period)
209:AliAnalysisTaskSE(name),
9aed5fcd 210 fPHOSGeo(0x0),
42800acf 211 fOutputContainer(0x0),
212 fEvent(0x0),
213 fEventESD(0x0),
214 fEventAOD(0x0),
215 fEventHandler(0),
216 fCaloPhotonsPHOS(0x0),
217 fTracksTPC(0x0),
218 fCaloPhotonsPHOSLists(0x0),
219 fTracksTPCLists(0x0),
220 fRunNumber(-999),
221 fInternalRunNumber(0),
222 fPeriod(period),
9aed5fcd 223 fPHOSEvent(false),
42800acf 224 fMBEvent(false),
225 fNVtxZBins(10),
08a327bf 226 fVtxEdges(10),
42800acf 227 fCentEdges(10),
228 fCentNMixed(),
08a327bf 229 fNEMRPBins(6),
42800acf 230 fAssocBins(),
67ef08bd 231 fVertexVector(),
42800acf 232 fVtxBin(0),
233 fCentralityEstimator("V0M"),
234 fCentrality(0.),
235 fCentBin(0),
236 fHaveTPCRP(0),
237 fRP(0.),
238 fEMRPBin(0),
239 fMaxAbsVertexZ(10.),
240 fCentralityLowLimit(0.),
241 fCentralityHightLimit(90),
42800acf 242 fMinClusterEnergy(0.3),
243 fMinBCDistance(0),
244 fMinNCells(3),
245 fMinM02(0.2),
246 fTOFCutEnabled(1),
247 fTOFCut(100.e-9),
08a327bf 248 fUseMassWindowParametrisation(true),
249 fMassInvMeanMin(0.13),
250 fMassInvMeanMax(0.145),
251 fNSigmaWidth(0.),
252 fUseEfficiency(true),
253 fESDtrackCuts(0x0),
254 fSelectHybridTracks(0),
255 fTrackStatus(0),
256 fTrackFilterMask(786),
257 fSelectSPDHitTracks(0),
258 fSelectFractionTPCSharedClusters(kTRUE),
259 fCutTPCSharedClustersFraction(0.4)
67ef08bd 260{
261 // Constructor
262 // Output slots #0 write into a TH1 container
263 DefineOutput(1,THashList::Class());
264
42800acf 265 fMassMean[0] = 1.00796e-05 ;
266 fMassMean[1] = 0.136096 ;
08a327bf 267 fMassSigma[0] = 0.00383029 ;
268 fMassSigma[1] = 0.0041709 ;
269 fMassSigma[2] = 0.00468736 ;
42800acf 270
67ef08bd 271 const Int_t nPtAssoc=10 ;
272 Double_t ptAssocBins[nPtAssoc]={0.,0.5,1.0,1.5,2.0,3.,5.,7.,10.,16} ;
273 fAssocBins.Set(nPtAssoc,ptAssocBins) ;
274
275 const int nbins = 9;
276 Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
277 TArrayD centEdges(nbins+1, edges);
278 Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
279 //Int_t nMixed[nbins] = {100,100,100,100,100,100,100,100,100};
280 TArrayI centNMixed(nbins, nMixed);
281 SetCentralityBinning(centEdges, centNMixed);
08a327bf 282 SetVertexBinning();
67ef08bd 283
42800acf 284 fVertex[0] = 0;
285 fVertex[1] = 0;
286 fVertex[2] = 0;
d2c19ce3 287
9aed5fcd 288 SetGeometry();
289
290 ZeroingVariables();
67ef08bd 291}
42800acf 292
67ef08bd 293//_______________________________________________________________________________
294AliPHOSCorrelations::~AliPHOSCorrelations()
295{
67ef08bd 296 if(fCaloPhotonsPHOS){
297 delete fCaloPhotonsPHOS;
298 fCaloPhotonsPHOS=0x0;
299 }
300
301 if(fTracksTPC){
302 delete fTracksTPC;
303 fTracksTPC=0x0;
304 }
305
306 if(fCaloPhotonsPHOSLists){
307 fCaloPhotonsPHOSLists->SetOwner() ;
308 delete fCaloPhotonsPHOSLists;
309 fCaloPhotonsPHOSLists=0x0;
310 }
311
312 if(fTracksTPCLists){
313 fTracksTPCLists->SetOwner() ;
314 delete fTracksTPCLists;
315 fTracksTPCLists=0x0 ;
316 }
317
318 if( fESDtrackCuts){
319 delete fESDtrackCuts;
320 fESDtrackCuts=0x0 ;
321 }
322
323 if(fOutputContainer){
324 delete fOutputContainer;
325 fOutputContainer=0x0;
d2c19ce3 326 }
67ef08bd 327}
42800acf 328
67ef08bd 329//_______________________________________________________________________________
330void AliPHOSCorrelations::UserCreateOutputObjects()
331{
332 // Create histograms
333 // Called once
08a327bf 334 const Int_t nRuns = 200 ;
335 const Int_t ptMult = 300 ;
336 const Double_t ptMin = 0. ;
337 const Double_t ptMax = 30. ;
67ef08bd 338
339 // Create histograms
340 if(fOutputContainer != NULL) { delete fOutputContainer; }
341 fOutputContainer = new THashList();
342 fOutputContainer->SetOwner(kTRUE);
343
7deabbe7 344 // Event selection
08a327bf 345 fOutputContainer->Add(new TH1F( "hTriggerPassedEvents", "Event selection passed Cuts", 60, 0., 60.)) ;
7deabbe7 346 // Analysis event's progress
08a327bf 347 fOutputContainer->Add(new TH1F( "hTotSelEvents", "Event selection", 15, 0., 15 )) ;
348 fOutputContainer->Add(new TH2F( "hSelEvents", "Event selection", kTotalSelected+1, 0., double(kTotalSelected+1), nRuns,0., float(nRuns) )) ;
7deabbe7 349 // Centrality, Reaction plane selection
08a327bf 350 fOutputContainer->Add(new TH2F( "hCentrality", "Event centrality of all events", 100, 0., 100., nRuns,0., float(nRuns) )) ;
351 fOutputContainer->Add(new TH2F( "hCentralityTriggerEvent", "Event centrality trigger events", 100, 0., 100., nRuns,0., float(nRuns) )) ;
352 fOutputContainer->Add(new TH2F( "hCentralityMBEvent", "Event centrality MB events", 100, 0., 100., nRuns,0., float(nRuns) )) ;
353 fOutputContainer->Add(new TH2F( "phiRPflat", "RP distribution with TPC flat", 100, 0., 2.*TMath::Pi(), 20, 0., 100. )) ;
354
355 fOutputContainer->Add(new TH1F( "hCentralityBining", "Event centrality bining of all events", GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
356 fOutputContainer->Add(new TH1F( "hCentralityBiningTrigger", "Event centrality bining of trigger events", GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
357 fOutputContainer->Add(new TH1F( "hCentralityBiningMB", "Event centrality bining of MB events", GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
358
359 fOutputContainer->Add(new TH1F( "phiRPflatBining", "Event RP bining of all events", GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
360 fOutputContainer->Add(new TH1F( "phiRPflatBiningTrigger", "Event RP bining of trigger events", GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
361 fOutputContainer->Add(new TH1F( "phiRPflatBiningMB", "Event RP bining of MB events", GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
362
363 fOutputContainer->Add(new TH1F( "hVertexZBining", "Event vertex bining of all events", GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
364 fOutputContainer->Add(new TH1F( "hVertexZBiningTrigger", "Event vertex bining of trigger events", GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
365 fOutputContainer->Add(new TH1F( "hVertexZBiningMB", "Event vertex bining of MB events", GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
366
7deabbe7 367 // Mass selection
08a327bf 368 fOutputContainer->Add(new TH2F( "massWindow", "mean & sigma", 100., 0.095, 0.185, 500, 0., 0.05 )) ;
369 fOutputContainer->Add(new TH1F( "massWindowPass", "Mass selection", 10, 0., 10. )) ;
7deabbe7 370 // Cluster multiplisity
08a327bf 371 fOutputContainer->Add(new TH2F( "hCluEvsClu", "ClusterMult vs E", 200, 0., 10., 100, 0., 100. )) ;
d2c19ce3 372
67ef08bd 373
374 // Set hists, with track's and cluster's angle distributions.
9aed5fcd 375 SetHistPtNumTrigger(ptMult, ptMin, ptMax);
08a327bf 376 SetHistEtaPhi(ptMult, ptMin, ptMax);
d2c19ce3 377 SetHistPHOSClusterMap();
42800acf 378 SetHistMass (ptMult, ptMin, ptMax);
9aed5fcd 379 SetHistPtAssoc(ptMult, ptMin, ptMax);
67ef08bd 380
381 // Setup photon lists
382 Int_t kapacity = fNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
383 fCaloPhotonsPHOSLists = new TObjArray(kapacity);
384 fCaloPhotonsPHOSLists->SetOwner();
385
386 fTracksTPCLists = new TObjArray(kapacity);
387 fTracksTPCLists->SetOwner();
388
389 PostData(1, fOutputContainer);
390}
42800acf 391
67ef08bd 392//_______________________________________________________________________________
08a327bf 393void AliPHOSCorrelations::SetHistPtNumTrigger(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
9aed5fcd 394{
395 TString spid[4]={"all","cpv","disp","both"} ;
396 for(Int_t ipid=0; ipid<4; ipid++)
397 {
42800acf 398 fOutputContainer->Add(new TH1F( Form("nTrigger_%s", spid[ipid].Data()),
399 Form("Num of trigger particle %s", spid[ipid].Data()),
400 ptMult+300, ptMin, ptMax ) );
4ae540f6 401 TH1F *h = static_cast<TH1F*>(fOutputContainer->Last()) ;
83aecbbd 402 h->Sumw2();
403 h->GetXaxis()->SetTitle("Pt [GEV]");
404 //h->GetYaxis()->SetTitle("#varepsilon"); // 1/efficiensy
405 }
406 for(Int_t ipid=0; ipid<4; ipid++)
407 {
42800acf 408 fOutputContainer->Add(new TH1F( Form("nTrigger_%s_MB", spid[ipid].Data()),
409 Form("Num of trigger particle %s", spid[ipid].Data()),
410 ptMult+300, ptMin, ptMax ) );
83aecbbd 411 TH1F *h = static_cast<TH1F*>(fOutputContainer->Last()) ;
412 h->Sumw2();
9aed5fcd 413 h->GetXaxis()->SetTitle("Pt [GEV]");
4ae540f6 414 //h->GetYaxis()->SetTitle("#varepsilon"); // 1/efficiensy
9aed5fcd 415 }
416}
42800acf 417
9aed5fcd 418//_______________________________________________________________________________
08a327bf 419void AliPHOSCorrelations::SetHistEtaPhi(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
67ef08bd 420{
d2c19ce3 421 // Set hists, with track's and cluster's angle distributions.
67ef08bd 422
423 Float_t pi = TMath::Pi();
424
425 //===
42800acf 426 fOutputContainer->Add(new TH2F( "clu_phieta","Cluster's #phi & #eta distribution",
427 300, double(-1.8), double(-0.6),
428 300, double(-0.2), double(0.2) ) );
67ef08bd 429 TH2F * h = static_cast<TH2F*>(fOutputContainer->Last()) ;
42800acf 430 h->GetXaxis()->SetTitle("#phi [rad]");
d2c19ce3 431 h->GetYaxis()->SetTitle("#eta");
67ef08bd 432
433 //===
42800acf 434 fOutputContainer->Add(new TH2F( "clusingle_phieta","Cluster's #phi & #eta distribution",
435 300, double(-1.8), double(-0.6),
436 300, double(-0.2), double(0.2) ) );
67ef08bd 437 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
42800acf 438 h->GetXaxis()->SetTitle("#phi [rad]");
d2c19ce3 439 h->GetYaxis()->SetTitle("#eta");
67ef08bd 440
441 //===
08a327bf 442 fOutputContainer->Add(new TH3F( "track_phieta","TPC track's #phi & #eta distribution",
42800acf 443 200, double(-pi-0.3), double(pi+0.3),
08a327bf 444 200, double(-0.9), double(0.9),
445 ptMult, ptMin, ptMax ) );
446 TH3F* h3 = static_cast<TH3F*>(fOutputContainer->FindObject("track_phieta")) ;
447 h3->GetXaxis()->SetTitle("#phi [rad]");
448 h3->GetYaxis()->SetTitle("#eta");
449 h3->GetZaxis()->SetTitle("Pt [GeV/c]");
67ef08bd 450}
42800acf 451
67ef08bd 452//_______________________________________________________________________________
08a327bf 453void AliPHOSCorrelations::SetHistMass(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
67ef08bd 454{
42800acf 455 // Set mass histograms.
456
83aecbbd 457 Double_t binMult = 400;
458 Double_t massMin = 0.0;
459 Double_t massMax = 0.4;
67ef08bd 460
9aed5fcd 461 TString spid[4]={"all","cpv","disp","both"} ;
67ef08bd 462
9aed5fcd 463 TH2F * h;
67ef08bd 464
9aed5fcd 465 for(Int_t ipid=0; ipid<4; ipid++)
466 {
467 // Real ++++++++++++++++++++++++++++++
67ef08bd 468
42800acf 469 fOutputContainer->Add(new TH2F(Form("%s_mpt", spid[ipid].Data() ), "Real",
470 binMult, massMin, massMax,
471 ptMult, ptMin, ptMax ) );
9aed5fcd 472 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
83aecbbd 473 h->Sumw2();
9aed5fcd 474 h->GetXaxis()->SetTitle("Mass [GeV]");
475 h->GetYaxis()->SetTitle("Pt [GEV]");
476
477 // MIX +++++++++++++++++++++++++
478
42800acf 479 fOutputContainer->Add(new TH2F(Form("mix_%s_mpt", spid[ipid].Data() ), "Mix",
480 binMult, massMin, massMax,
481 ptMult, ptMin, ptMax ) );
9aed5fcd 482 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
83aecbbd 483 h->Sumw2();
42800acf 484 h->GetXaxis()->SetTitle("Mass [GeV]");
9aed5fcd 485 h->GetYaxis()->SetTitle("Pt [GEV]");
486 }
487
488 // Calibration PHOS Module Pi0peak {REAL}
67ef08bd 489 for(Int_t mod=1; mod<4; mod++){
42800acf 490 fOutputContainer->Add(new TH2F(Form( "both%d_mpt",mod), Form("Both cuts (CPV + Disp) mod[%d]",mod),
491 binMult, massMin, massMax,
492 ptMult, ptMin, ptMax ) );
67ef08bd 493 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
83aecbbd 494 h->Sumw2();
42800acf 495 h->GetXaxis()->SetTitle("Mass [GeV]");
67ef08bd 496 h->GetYaxis()->SetTitle("Pt [GEV]");
497
9aed5fcd 498 // Calibration PHOS Module Pi0peak {MIX}
42800acf 499 fOutputContainer->Add(new TH2F(Form( "mix_both%d_mpt",mod), Form(" Both cuts (CPV + Disp) mod[%d]",mod),
500 binMult, massMin, massMax,
501 ptMult, ptMin, ptMax ) );
83aecbbd 502 h = static_cast<TH2F*>(fOutputContainer->Last()) ;
503 h->Sumw2();
42800acf 504 h->GetXaxis()->SetTitle("Mass [GeV]");
67ef08bd 505 h->GetYaxis()->SetTitle("Pt [GEV]");
42800acf 506
67ef08bd 507 }
67ef08bd 508}
42800acf 509
67ef08bd 510//_______________________________________________________________________________
08a327bf 511void AliPHOSCorrelations::SetHistPtAssoc(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
67ef08bd 512{
513 Double_t pi = TMath::Pi();
514
42800acf 515 Int_t PhiMult = 100 ;
516 Float_t PhiMin = -0.5*pi ;
517 Float_t PhiMax = 1.5*pi ;
518 Int_t EtaMult = 20 ;
519 Float_t EtaMin = -1. ;
520 Float_t EtaMax = 1. ;
67ef08bd 521
522 TString spid[4]={"all","cpv","disp","both"} ;
523
524 for (int i = 0; i<fAssocBins.GetSize()-1; i++){
525 for(Int_t ipid=0; ipid<4; ipid++){
42800acf 526 // Main histo for ConsiderPi0s().
527 fOutputContainer->Add(new TH3F(Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
528 Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
529 ptMult, ptMin, ptMax,
530 PhiMult, PhiMin, PhiMax,
531 EtaMult, EtaMin, EtaMax ) );
67ef08bd 532 TH3F * h = static_cast<TH3F*>(fOutputContainer->Last()) ;
83aecbbd 533 h->Sumw2();
534 h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
535 h->GetYaxis()->SetTitle("#phi [rad]");
536 h->GetZaxis()->SetTitle("#eta");
537
42800acf 538 // For ConsiderPi0s_MBSelection().
539 fOutputContainer->Add(new TH3F(Form("%s_ptphieta_ptAssoc_%3.1f_MB", spid[ipid].Data(), fAssocBins.At(i+1)),
540 Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
541 ptMult, ptMin, ptMax,
542 PhiMult, PhiMin, PhiMax,
543 EtaMult, EtaMin, EtaMax ) );
83aecbbd 544 h = static_cast<TH3F*>(fOutputContainer->Last()) ;
545 h->Sumw2();
546 h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
67ef08bd 547 h->GetYaxis()->SetTitle("#phi [rad]");
d2c19ce3 548 h->GetZaxis()->SetTitle("#eta");
67ef08bd 549
83aecbbd 550 // For Mixed events in ConsiderTracksMix()
42800acf 551 fOutputContainer->Add(new TH3F(Form("mix_%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
552 Form("Mixed %s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
553 ptMult, ptMin, ptMax,
554 PhiMult, PhiMin, PhiMax,
555 EtaMult, EtaMin, EtaMax ) );
67ef08bd 556 h = static_cast<TH3F*>(fOutputContainer->Last()) ;
83aecbbd 557 h->Sumw2();
558 h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
67ef08bd 559 h->GetYaxis()->SetTitle("#phi [rad]");
d2c19ce3 560 h->GetZaxis()->SetTitle("#eta");
67ef08bd 561 }
562 }
563}
d2c19ce3 564
42800acf 565//_______________________________________________________________________________
d2c19ce3 566void AliPHOSCorrelations::SetHistPHOSClusterMap()
567{
42800acf 568 // Cluster X/Z/E distribution.
7deabbe7 569 for(int i = 0; i<5; i++)
d2c19ce3 570 {
42800acf 571 fOutputContainer->Add(new TH3F( Form("QA_cluXZE_mod%i", i), Form("PHOS Clusters XZE distribution of module %i", i),
572 70, 0, 70,
573 60, 0, 60,
574 200, 0, 20 ) );
d2c19ce3 575 TH3F *h = static_cast<TH3F*>(fOutputContainer->Last()) ;
83aecbbd 576 h->GetXaxis()->SetTitle("X");
d2c19ce3 577 h->GetYaxis()->SetTitle("Z");
578 h->GetZaxis()->SetTitle("E");
579 }
580}
42800acf 581
67ef08bd 582//_______________________________________________________________________________
583void AliPHOSCorrelations::UserExec(Option_t *)
584{
d2c19ce3 585 // Main loop, called for each event analyze ESD/AOD
67ef08bd 586 // Step 0: Event Objects
7deabbe7 587 LogProgress(0);
42800acf 588
d2c19ce3 589 fEvent = InputEvent();
590 if( ! fEvent )
591 {
592 AliError("Event could not be retrieved");
593 PostData(1, fOutputContainer);
594 return ;
595 }
7deabbe7 596 LogProgress(1);
67ef08bd 597
67ef08bd 598 // Step 1(done once):
599 if( fRunNumber != fEvent->GetRunNumber() )
600 {
601 fRunNumber = fEvent->GetRunNumber();
602 fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
603 SetESDTrackCuts();
604 }
42800acf 605 LogProgress(2);
606
607 // Step 2: Preparation variables for new event
608 ZeroingVariables();
67ef08bd 609
42800acf 610 fEventESD = dynamic_cast<AliESDEvent*>(fEvent);
611 fEventAOD = dynamic_cast<AliAODEvent*>(fEvent);
612
613 // Get Event-Handler for the trigger information
7deabbe7 614 fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
615 if (!fEventHandler)
67ef08bd 616 {
7deabbe7 617 AliError("Could not get InputHandler");
67ef08bd 618 PostData(1, fOutputContainer);
619 return; // Reject!
620 }
42800acf 621 LogProgress(3);
67ef08bd 622
42800acf 623 // Step 3: Event trigger selection
624 // fPHOSEvent, fMBEvent
7deabbe7 625 if( RejectTriggerMaskSelection() )
626 {
627 PostData(1, fOutputContainer);
628 return; // Reject!
629 }
42800acf 630 LogProgress(4);
631
632 // Step 4: Vertex
67ef08bd 633 // fVertex, fVertexVector, fVtxBin
634 SetVertex();
635 if( RejectEventVertex() )
636 {
637 PostData(1, fOutputContainer);
638 return; // Reject!
639 }
42800acf 640 LogProgress(5);
67ef08bd 641
42800acf 642 // Step 5: Centrality
67ef08bd 643 // fCentrality, fCentBin
644 SetCentrality();
645 if( RejectEventCentrality() )
646 {
647 PostData(1, fOutputContainer);
648 return; // Reject!
649 }
42800acf 650 LogProgress(6);
651 if(fPHOSEvent) FillHistogram( "hCentralityTriggerEvent", fCentrality, fInternalRunNumber-0.5 ) ;
652 if(fMBEvent) FillHistogram( "hCentralityMBEvent", fCentrality, fInternalRunNumber-0.5 ) ;
653 FillHistogram( "hCentrality", fCentrality, fInternalRunNumber-0.5 ) ;
67ef08bd 654
42800acf 655 // Step 6: Reaction Plane
67ef08bd 656 // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
657 EvalReactionPlane();
658 fEMRPBin = GetRPBin();
67ef08bd 659
42800acf 660 // Step 7: Event Photons (PHOS Clusters) selection
67ef08bd 661 SelectPhotonClusters();
9aed5fcd 662 if( ! fCaloPhotonsPHOS->GetEntriesFast() )
663 LogSelection(kHasPHOSClusters, fInternalRunNumber);
67ef08bd 664
42800acf 665 // Step 8: Event Associated particles (TPC Tracks) selection
67ef08bd 666 SelectAccosiatedTracks();
67ef08bd 667 if( ! fTracksTPC->GetEntriesFast() )
7deabbe7 668 LogSelection(kHasTPCTracks, fInternalRunNumber);
67ef08bd 669 LogSelection(kTotalSelected, fInternalRunNumber);
4ae540f6 670
08a327bf 671 // Step 9: Fill TPC's track mask and control bining hists.
7deabbe7 672 FillTrackEtaPhi();
08a327bf 673 FillEventBiningProperties();
42800acf 674 LogProgress(7);
4ae540f6 675
42800acf 676 // Step 10: Extract one most energetic pi0 candidate in this event.
677 SelectTriggerPi0ME();
83aecbbd 678
42800acf 679 // Step 11: Start correlation analysis.
680 if (fPHOSEvent)
681 {
682 ConsiderPi0s(); // Consider the most energetic Pi0 in this event with all tracks of this event.
683 LogProgress(8);
7deabbe7 684 }
42800acf 685
686 if(fPeriod == kLHC13 && fMBEvent)
7deabbe7 687 {
42800acf 688 ConsiderPi0s_MBSelection();
689 LogProgress(9);
7deabbe7 690 }
4ae540f6 691
7deabbe7 692 // Filling mixing histograms:
42800acf 693 if (fMBEvent)
7deabbe7 694 {
42800acf 695 ConsiderPi0sMix(); // Make background for extracting pi0 mass.
08a327bf 696 ConsiderTracksMix(); // Compare only one most energetic pi0 candidate with all tracks from previous MB events.
4ae540f6 697
08a327bf 698 // Update pull using MB events only!
42800acf 699 UpdatePhotonLists(); // Updating pull of photons.
700 UpdateTrackLists(); // Updating pull of tracks.
701 LogProgress(10);
7deabbe7 702 }
67ef08bd 703
08a327bf 704 LogProgress(11);
67ef08bd 705 // Post output data.
706 PostData(1, fOutputContainer);
707}
42800acf 708
67ef08bd 709//_______________________________________________________________________________
710void AliPHOSCorrelations::SetESDTrackCuts()
711{
42800acf 712 if( fEventESD )
713 {
714 // Create ESD track cut
715 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
716 fESDtrackCuts->SetRequireTPCRefit(kTRUE) ;
717 }
67ef08bd 718}
42800acf 719
67ef08bd 720//_______________________________________________________________________________
08a327bf 721Int_t AliPHOSCorrelations::ConvertToInternalRunNumber(const Int_t& run) const
42800acf 722{
723 // Manual setup using data from logbook.
08a327bf 724 if(fPeriod== kLHC11h)
725 {
726 switch(run)
727 {
728 case 170593 : return 179 ;
729 case 170572 : return 178 ;
730 case 170556 : return 177 ;
731 case 170552 : return 176 ;
732 case 170546 : return 175 ;
733 case 170390 : return 174 ;
734 case 170389 : return 173 ;
735 case 170388 : return 172 ;
736 case 170387 : return 171 ;
737 case 170315 : return 170 ;
738 case 170313 : return 169 ;
739 case 170312 : return 168 ;
740 case 170311 : return 167 ;
741 case 170309 : return 166 ;
742 case 170308 : return 165 ;
743 case 170306 : return 164 ;
744 case 170270 : return 163 ;
745 case 170269 : return 162 ;
746 case 170268 : return 161 ;
747 case 170267 : return 160 ;
748 case 170264 : return 159 ;
749 case 170230 : return 158 ;
750 case 170228 : return 157 ;
751 case 170208 : return 156 ;
752 case 170207 : return 155 ;
753 case 170205 : return 154 ;
754 case 170204 : return 153 ;
755 case 170203 : return 152 ;
756 case 170195 : return 151 ;
757 case 170193 : return 150 ;
758 case 170163 : return 149 ;
759 case 170162 : return 148 ;
760 case 170159 : return 147 ;
761 case 170155 : return 146 ;
762 case 170152 : return 145 ;
763 case 170091 : return 144 ;
764 case 170089 : return 143 ;
765 case 170088 : return 142 ;
766 case 170085 : return 141 ;
767 case 170084 : return 140 ;
768 case 170083 : return 139 ;
769 case 170081 : return 138 ;
770 case 170040 : return 137 ;
771 case 170038 : return 136 ;
772 case 170036 : return 135 ;
773 case 170027 : return 134 ;
774 case 169981 : return 133 ;
775 case 169975 : return 132 ;
776 case 169969 : return 131 ;
777 case 169965 : return 130 ;
778 case 169961 : return 129 ;
779 case 169956 : return 128 ;
780 case 169926 : return 127 ;
781 case 169924 : return 126 ;
782 case 169923 : return 125 ;
783 case 169922 : return 124 ;
784 case 169919 : return 123 ;
785 case 169918 : return 122 ;
786 case 169914 : return 121 ;
787 case 169859 : return 120 ;
788 case 169858 : return 119 ;
789 case 169855 : return 118 ;
790 case 169846 : return 117 ;
791 case 169838 : return 116 ;
792 case 169837 : return 115 ;
793 case 169835 : return 114 ;
794 case 169683 : return 113 ;
795 case 169628 : return 112 ;
796 case 169591 : return 111 ;
797 case 169590 : return 110 ;
798 case 169588 : return 109 ;
799 case 169587 : return 108 ;
800 case 169586 : return 107 ;
801 case 169584 : return 106 ;
802 case 169557 : return 105 ;
803 case 169555 : return 104 ;
804 case 169554 : return 103 ;
805 case 169553 : return 102 ;
806 case 169550 : return 101 ;
807 case 169515 : return 100 ;
808 case 169512 : return 99 ;
809 case 169506 : return 98 ;
810 case 169504 : return 97 ;
811 case 169498 : return 96 ;
812 case 169475 : return 95 ;
813 case 169420 : return 94 ;
814 case 169419 : return 93 ;
815 case 169418 : return 92 ;
816 case 169417 : return 91 ;
817 case 169415 : return 90 ;
818 case 169411 : return 89 ;
819 case 169238 : return 88 ;
820 case 169236 : return 87 ;
821 case 169167 : return 86 ;
822 case 169160 : return 85 ;
823 case 169156 : return 84 ;
824 case 169148 : return 83 ;
825 case 169145 : return 82 ;
826 case 169144 : return 81 ;
827 case 169143 : return 80 ;
828 case 169138 : return 79 ;
829 case 169099 : return 78 ;
830 case 169094 : return 77 ;
831 case 169091 : return 76 ;
832 case 169045 : return 75 ;
833 case 169044 : return 74 ;
834 case 169040 : return 73 ;
835 case 169035 : return 72 ;
836 case 168992 : return 71 ;
837 case 168988 : return 70 ;
838 case 168984 : return 69 ;
839 case 168826 : return 68 ;
840 case 168777 : return 67 ;
841 case 168514 : return 66 ;
842 case 168512 : return 65 ;
843 case 168511 : return 64 ;
844 case 168467 : return 63 ;
845 case 168464 : return 62 ;
846 case 168461 : return 61 ;
847 case 168460 : return 60 ;
848 case 168458 : return 59 ;
849 case 168362 : return 58 ;
850 case 168361 : return 57 ;
851 case 168356 : return 56 ;
852 case 168342 : return 55 ;
853 case 168341 : return 54 ;
854 case 168325 : return 53 ;
855 case 168322 : return 52 ;
856 case 168318 : return 51 ;
857 case 168311 : return 50 ;
858 case 168310 : return 49 ;
859 case 168213 : return 48 ;
860 case 168212 : return 47 ;
861 case 168208 : return 46 ;
862 case 168207 : return 45 ;
863 case 168206 : return 44 ;
864 case 168205 : return 43 ;
865 case 168204 : return 42 ;
866 case 168203 : return 41 ;
867 case 168181 : return 40 ;
868 case 168177 : return 39 ;
869 case 168175 : return 38 ;
870 case 168173 : return 37 ;
871 case 168172 : return 36 ;
872 case 168171 : return 35 ;
873 case 168115 : return 34 ;
874 case 168108 : return 33 ;
875 case 168107 : return 32 ;
876 case 168105 : return 31 ;
877 case 168104 : return 30 ;
878 case 168103 : return 29 ;
879 case 168076 : return 28 ;
880 case 168069 : return 27 ;
881 case 168068 : return 26 ;
882 case 168066 : return 25 ;
883 case 167988 : return 24 ;
884 case 167987 : return 23 ;
885 case 167986 : return 22 ;
886 case 167985 : return 21 ;
887 case 167921 : return 20 ;
888 case 167920 : return 19 ;
889 case 167915 : return 18 ;
890 case 167909 : return 17 ;
891 case 167903 : return 16 ;
892 case 167902 : return 15 ;
893 case 167818 : return 14 ;
894 case 167814 : return 13 ;
895 case 167813 : return 12 ;
896 case 167808 : return 11 ;
897 case 167807 : return 10 ;
898 case 167806 : return 9 ;
899 case 167713 : return 8 ;
900 case 167712 : return 7 ;
901 case 167711 : return 6 ;
902 case 167706 : return 5 ;
903 case 167693 : return 4 ;
904 case 166532 : return 3 ;
905 case 166530 : return 2 ;
906 case 166529 : return 1 ;
907
908 default : return 199;
909 }
910 }
911
912 if(fPeriod== kLHC10h)
913 {
914 switch(run)
915 {
916 case 139517 : return 137;
917 case 139514 : return 136;
918 case 139513 : return 135;
919 case 139511 : return 134;
920 case 139510 : return 133;
921 case 139507 : return 132;
922 case 139505 : return 131;
923 case 139504 : return 130;
924 case 139503 : return 129;
925 case 139470 : return 128;
926 case 139467 : return 127;
927 case 139466 : return 126;
928 case 139465 : return 125;
929 case 139440 : return 124;
930 case 139439 : return 123;
931 case 139438 : return 122;
932 case 139437 : return 121;
933 case 139360 : return 120;
934 case 139329 : return 119;
935 case 139328 : return 118;
936 case 139314 : return 117;
937 case 139311 : return 116;
938 case 139310 : return 115;
939 case 139309 : return 114;
940 case 139308 : return 113;
941 case 139173 : return 112;
942 case 139172 : return 111;
943 case 139110 : return 110;
944 case 139107 : return 109;
945 case 139105 : return 108;
946 case 139104 : return 107;
947 case 139042 : return 106;
948 case 139038 : return 105;
949 case 139037 : return 104;
950 case 139036 : return 103;
951 case 139029 : return 102;
952 case 139028 : return 101;
953 case 138983 : return 100;
954 case 138982 : return 99;
955 case 138980 : return 98;
956 case 138979 : return 97;
957 case 138978 : return 96;
958 case 138977 : return 95;
959 case 138976 : return 94;
960 case 138973 : return 93;
961 case 138972 : return 92;
962 case 138965 : return 91;
963 case 138924 : return 90;
964 case 138872 : return 89;
965 case 138871 : return 88;
966 case 138870 : return 87;
967 case 138837 : return 86;
968 case 138830 : return 85;
969 case 138828 : return 84;
970 case 138826 : return 83;
971 case 138796 : return 82;
972 case 138795 : return 81;
973 case 138742 : return 80;
974 case 138732 : return 79;
975 case 138730 : return 78;
976 case 138666 : return 77;
977 case 138662 : return 76;
978 case 138653 : return 75;
979 case 138652 : return 74;
980 case 138638 : return 73;
981 case 138624 : return 72;
982 case 138621 : return 71;
983 case 138583 : return 70;
984 case 138582 : return 69;
985 case 138579 : return 68;
986 case 138578 : return 67;
987 case 138534 : return 66;
988 case 138469 : return 65;
989 case 138442 : return 64;
990 case 138439 : return 63;
991 case 138438 : return 62;
992 case 138396 : return 61;
993 case 138364 : return 60;
994 case 138359 : return 59;
995 case 138275 : return 58;
996 case 138225 : return 57;
997 case 138201 : return 56;
998 case 138200 : return 55;
999 case 138197 : return 54;
1000 case 138192 : return 53;
1001 case 138190 : return 52;
1002 case 138154 : return 51;
1003 case 138153 : return 50;
1004 case 138151 : return 49;
1005 case 138150 : return 48;
1006 case 138126 : return 47;
1007 case 138125 : return 46;
1008 case 137848 : return 45;
1009 case 137847 : return 44;
1010 case 137844 : return 43;
1011 case 137843 : return 42;
1012 case 137752 : return 41;
1013 case 137751 : return 40;
1014 case 137748 : return 39;
1015 case 137724 : return 38;
1016 case 137722 : return 37;
1017 case 137718 : return 36;
1018 case 137704 : return 35;
1019 case 137693 : return 34;
1020 case 137692 : return 33;
1021 case 137691 : return 32;
1022 case 137689 : return 31;
1023 case 137686 : return 30;
1024 case 137685 : return 29;
1025 case 137639 : return 28;
1026 case 137638 : return 27;
1027 case 137608 : return 26;
1028 case 137595 : return 25;
1029 case 137549 : return 24;
1030 case 137546 : return 23;
1031 case 137544 : return 22;
1032 case 137541 : return 21;
1033 case 137539 : return 20;
1034 case 137531 : return 19;
1035 case 137530 : return 18;
1036 case 137443 : return 17;
1037 case 137441 : return 16;
1038 case 137440 : return 15;
1039 case 137439 : return 14;
1040 case 137434 : return 13;
1041 case 137432 : return 12;
1042 case 137431 : return 11;
1043 case 137430 : return 10;
1044 case 137366 : return 9;
1045 case 137243 : return 8;
1046 case 137236 : return 7;
1047 case 137235 : return 6;
1048 case 137232 : return 5;
1049 case 137231 : return 4;
1050 case 137165 : return 3;
1051 case 137162 : return 2;
1052 case 137161 : return 1;
1053 default : return 199;
1054 }
1055 }
1056
1057 if( kLHC13 == fPeriod )
1058 {
1059 switch(run)
1060 {
1061 case 195344 : return 1;
1062 case 195346 : return 2;
1063 case 195351 : return 3;
1064 case 195389 : return 4;
1065 case 195390 : return 5;
1066 case 195391 : return 6;
1067 case 195478 : return 7;
1068 case 195479 : return 8;
1069 case 195480 : return 9;
1070 case 195481 : return 10;
1071 case 195482 : return 11;
1072 case 195483 : return 12;
1073 case 195529 : return 13;
1074 case 195531 : return 14;
1075 case 195532 : return 15;
1076 case 195566 : return 16;
1077 case 195567 : return 17;
1078 case 195568 : return 18;
1079 case 195592 : return 19;
1080 case 195593 : return 20;
1081 case 195596 : return 21;
1082 case 195633 : return 22;
1083 case 195635 : return 23;
1084 case 195644 : return 24;
1085 case 195673 : return 25;
1086 case 195675 : return 26;
1087 case 195676 : return 27;
1088 case 195677 : return 28;
1089 case 195681 : return 29;
1090 case 195682 : return 30;
1091 case 195720 : return 31;
1092 case 195721 : return 32;
1093 case 195722 : return 33;
1094 case 195724 : return 34;
1095 case 195725 : return 34;
1096 case 195726 : return 35;
1097 case 195727 : return 36;
1098 case 195760 : return 37;
1099 case 195761 : return 38;
1100 case 195765 : return 39;
1101 case 195767 : return 40;
1102 case 195783 : return 41;
1103 case 195787 : return 42;
1104 case 195826 : return 43;
1105 case 195827 : return 44;
1106 case 195829 : return 45;
1107 case 195830 : return 46;
1108 case 195831 : return 47;
1109 case 195867 : return 48;
1110 case 195869 : return 49;
1111 case 195871 : return 50;
1112 case 195872 : return 51;
1113 case 195873 : return 52;
1114 case 195935 : return 53;
1115 case 195949 : return 54;
1116 case 195950 : return 55;
1117 case 195954 : return 56;
1118 case 195955 : return 57;
1119 case 195958 : return 58;
1120 case 195989 : return 59;
1121 case 195994 : return 60;
1122 case 195998 : return 61;
1123 case 196000 : return 62;
1124 case 196006 : return 63;
1125 case 196085 : return 64;
1126 case 196089 : return 65;
1127 case 196090 : return 66;
1128 case 196091 : return 67;
1129 case 196099 : return 68;
1130 case 196105 : return 69;
1131 case 196107 : return 70;
1132 case 196185 : return 71;
1133 case 196187 : return 72;
1134 case 196194 : return 73;
1135 case 196197 : return 74;
1136 case 196199 : return 75;
1137 case 196200 : return 76;
1138 case 196201 : return 77;
1139 case 196203 : return 78;
1140 case 196208 : return 79;
1141 case 196214 : return 80;
1142 case 196308 : return 81;
1143 case 196309 : return 82;
1144 case 196310 : return 83;
1145 case 196311 : return 84;
1146 case 196433 : return 85;
1147 case 196474 : return 86;
1148 case 196475 : return 87;
1149 case 196477 : return 88;
1150 case 196528 : return 89;
1151 case 196533 : return 90;
1152 case 196535 : return 91;
1153 case 196563 : return 92;
1154 case 196564 : return 93;
1155 case 196566 : return 94;
1156 case 196568 : return 95;
1157 case 196601 : return 96;
1158 case 196605 : return 97;
1159 case 196608 : return 98;
1160 case 196646 : return 99;
1161 case 196648 : return 100;
1162 case 196701 : return 101;
1163 case 196702 : return 102;
1164 case 196703 : return 103;
1165 case 196706 : return 104;
1166 case 196714 : return 105;
1167 case 196720 : return 106;
1168 case 196721 : return 107;
1169 case 196722 : return 108;
1170 case 196772 : return 109;
1171 case 196773 : return 110;
1172 case 196774 : return 111;
1173 case 196869 : return 112;
1174 case 196870 : return 113;
1175 case 196874 : return 114;
1176 case 196876 : return 115;
1177 case 196965 : return 116;
1178 case 196967 : return 117;
1179 case 196972 : return 118;
1180 case 196973 : return 119;
1181 case 196974 : return 120;
1182 case 197003 : return 121;
1183 case 197011 : return 122;
1184 case 197012 : return 123;
1185 case 197015 : return 124;
1186 case 197027 : return 125;
1187 case 197031 : return 126;
1188 case 197089 : return 127;
1189 case 197090 : return 128;
1190 case 197091 : return 129;
1191 case 197092 : return 130;
1192 case 197094 : return 131;
1193 case 197098 : return 132;
1194 case 197099 : return 133;
1195 case 197138 : return 134;
1196 case 197139 : return 135;
1197 case 197142 : return 136;
1198 case 197143 : return 137;
1199 case 197144 : return 138;
1200 case 197145 : return 139;
1201 case 197146 : return 140;
1202 case 197147 : return 140;
1203 case 197148 : return 141;
1204 case 197149 : return 142;
1205 case 197150 : return 143;
1206 case 197152 : return 144;
1207 case 197153 : return 145;
1208 case 197184 : return 146;
1209 case 197189 : return 147;
1210 case 197247 : return 148;
1211 case 197248 : return 149;
1212 case 197254 : return 150;
1213 case 197255 : return 151;
1214 case 197256 : return 152;
1215 case 197258 : return 153;
1216 case 197260 : return 154;
1217 case 197296 : return 155;
1218 case 197297 : return 156;
1219 case 197298 : return 157;
1220 case 197299 : return 158;
1221 case 197300 : return 159;
1222 case 197302 : return 160;
1223 case 197341 : return 161;
1224 case 197342 : return 162;
1225 case 197348 : return 163;
1226 case 197349 : return 164;
1227 case 197351 : return 165;
1228 case 197386 : return 166;
1229 case 197387 : return 167;
1230 case 197388 : return 168;
1231 default : return 199;
1232 }
1233 }
1234
1235 if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) )
1236 {
1237 AliWarning("Period not defined");
1238 }
1239 return 1;
42800acf 1240}
67ef08bd 1241
1242//_______________________________________________________________________________
1243Bool_t AliPHOSCorrelations::RejectTriggerMaskSelection()
1244{
7deabbe7 1245 // Analyse trigger event and reject it if it not intresting.
67ef08bd 1246 const Bool_t REJECT = true;
1247 const Bool_t ACCEPT = false;
1248
7deabbe7 1249 if( fDebug >= 2 )
1250 AliInfo( Form("Event passed offline phos trigger test: %s ", fEvent->GetFiredTriggerClasses().Data() ) );
1251
08a327bf 1252 const Int_t physSelMask = fEventHandler->IsEventSelected();
1253
1254 Bool_t isAny = physSelMask & AliVEvent::kAny; // to accept any trigger
1255
1256 Bool_t isPHI1 = physSelMask & AliVEvent::kPHI1; // PHOS trigger, CINT1 suite
1257 Bool_t isPHI7 = physSelMask & AliVEvent::kPHI7; // PHOS trigger, CINT7 suite
1258 Bool_t isPHI8 = physSelMask & AliVEvent::kPHI8; // PHOS trigger, CINT8 suite
1259 Bool_t isCentral = physSelMask & AliVEvent::kCentral; // PbPb central collision trigger
1260 Bool_t isSemiCentral = physSelMask & AliVEvent::kSemiCentral; // PbPb semicentral collision trigger
1261 Bool_t isPHOSPb = physSelMask & AliVEvent::kPHOSPb; // PHOS trigger, CINT8 suite for PbPb
1262
1263 Bool_t isMB = physSelMask & AliVEvent::kMB; // Minimum bias trigger, i.e. interaction trigger, offline SPD or V0 selection
1264 Bool_t isINT7 = physSelMask & AliVEvent::kINT7; // V0AND trigger, offline V0 selection
1265 Bool_t isAnyINT = physSelMask & AliVEvent::kAnyINT; // to accept any interaction (aka minimum bias) trigger
1266
1267 // Other triggers
1268 Bool_t isMUON = physSelMask & AliVEvent::kMUON; // Muon trigger, offline SPD or V0 selection
1269 Bool_t isHighMult = physSelMask & AliVEvent::kHighMult; // High-multiplicity trigger (threshold defined online), offline SPD or V0 selection
1270 Bool_t isEMC1 = physSelMask & AliVEvent::kEMC1; // EMCAL trigger
1271 Bool_t isCINT5 = physSelMask & AliVEvent::kCINT5; // Minimum bias trigger without SPD. i.e. interaction trigger, offline V0 selection
1272 Bool_t isCMUS5 = physSelMask & AliVEvent::kCMUS5; // Muon trigger, offline V0 selection
1273 Bool_t isMUSPB = physSelMask & AliVEvent::kMUSPB; // idem for PbPb
1274 Bool_t isMUSH7 = physSelMask & AliVEvent::kMUSH7; // Muon trigger: high pt, single muon, offline V0 selection, CINT7 suite
1275 Bool_t isMUSHPB = physSelMask & AliVEvent::kMUSHPB; // idem for PbPb
1276 Bool_t isMUL7 = physSelMask & AliVEvent::kMUL7; // Muon trigger: like sign dimuon, offline V0 selection, CINT7 suite
1277 Bool_t isMuonLikePB = physSelMask & AliVEvent::kMuonLikePB; // idem for PbPb
1278 Bool_t isMUU7 = physSelMask & AliVEvent::kMUU7; // Muon trigger, unlike sign dimuon, offline V0 selection, CINT7 suite
1279 Bool_t isMuonUnlikePB = physSelMask & AliVEvent::kMuonUnlikePB; // idem for PbPb
1280 Bool_t isEMC7 = physSelMask & AliVEvent::kEMC7; // EMCAL trigger, CINT7 suite
1281 Bool_t isEMC8 = physSelMask & AliVEvent::kEMC8; // EMCAL trigger, CINT8 suite
1282 Bool_t isMUS7 = physSelMask & AliVEvent::kMUS7; // Muon trigger: low pt, single muon, offline V0 selection, CINT7 suite
1283
1284 Bool_t isEMCEJE = physSelMask & AliVEvent::kEMCEJE; // EMCAL jet patch trigger
1285 Bool_t isEMCEGA = physSelMask & AliVEvent::kEMCEGA; // EMCAL gamma trigger
1286
1287 Bool_t isDG5 = physSelMask & AliVEvent::kDG5; // Double gap diffractive
1288 Bool_t isZED = physSelMask & AliVEvent::kZED; // ZDC electromagnetic dissociation
1289 Bool_t isSPI7 = physSelMask & AliVEvent::kSPI7; // Power interaction trigger
1290 Bool_t isSPI = physSelMask & AliVEvent::kSPI; // Power interaction trigger
1291 Bool_t isINT8 = physSelMask & AliVEvent::kINT8; // CINT8 trigger: 0TVX (T0 vertex) triger
1292 Bool_t isMuonSingleLowPt8 = physSelMask & AliVEvent::kMuonSingleLowPt8; // Muon trigger : single muon, low pt, T0 selection, CINT8 suite
1293 Bool_t isMuonSingleHighPt8 = physSelMask & AliVEvent::kMuonSingleHighPt8; // Muon trigger : single muon, high pt, T0 selection, CINT8 suite
1294 Bool_t isMuonLikeLowPt8 = physSelMask & AliVEvent::kMuonLikeLowPt8; // Muon trigger : like sign muon, low pt, T0 selection, CINT8 suite
1295 Bool_t isMuonUnlikeLowPt8 = physSelMask & AliVEvent::kMuonUnlikeLowPt8; // Muon trigger : unlike sign muon, low pt, T0 selection, CINT8 suite
1296 Bool_t isMuonUnlikeLowPt0 = physSelMask & AliVEvent::kMuonUnlikeLowPt0; // Muon trigger : unlike sign muon, low pt, no additional L0 requirement
1297 Bool_t isUserDefined = physSelMask & AliVEvent::kUserDefined; // Set when custom trigger classes are set in AliPhysicsSelection, offline SPD or V0 selection
1298 Bool_t isTRD = physSelMask & AliVEvent::kTRD; // TRD trigger
1299 Bool_t isFastOnly = physSelMask & AliVEvent::kFastOnly; // The fast cluster fired. This bit is set in to addition another trigger bit, e.g. kMB
67ef08bd 1300
7deabbe7 1301 // All input events
08a327bf 1302 FillHistogram("hTriggerPassedEvents", 0 );
42800acf 1303 if ( isAny ) FillHistogram("hTriggerPassedEvents", 1.) ;
67ef08bd 1304
7deabbe7 1305 // PHOS events.
42800acf 1306 if ( isPHI1 ) FillHistogram("hTriggerPassedEvents", 2. );
1307 if ( isPHI7 ) FillHistogram("hTriggerPassedEvents", 3. );
1308 if ( isPHI8 ) FillHistogram("hTriggerPassedEvents", 4. );
1309 if ( isCentral ) FillHistogram("hTriggerPassedEvents", 5. );
1310 if ( isSemiCentral ) FillHistogram("hTriggerPassedEvents", 6. );
1311 if ( isPHOSPb ) FillHistogram("hTriggerPassedEvents", 7. );
67ef08bd 1312
7deabbe7 1313 // MB events.
42800acf 1314 if ( isMB ) FillHistogram("hTriggerPassedEvents", 8. );
1315 if ( isINT7 ) FillHistogram("hTriggerPassedEvents", 9. );
1316 if ( isAnyINT ) FillHistogram("hTriggerPassedEvents", 10. );
7deabbe7 1317
08a327bf 1318 Bool_t isTriggerEvent = false;
1319 Bool_t isMIXEvent = false;
7deabbe7 1320
08a327bf 1321 // TODO: Set false by default.
42800acf 1322 fPHOSEvent = false ;
1323 fMBEvent = false ;
83aecbbd 1324
08a327bf 1325 // Choosing triggers for different periods.
83aecbbd 1326 if(fPeriod == kLHC13)
7deabbe7 1327 {
08a327bf 1328 // Working: MB + TriggerPHOS events in real events; MB only in mixed events.
83aecbbd 1329 isTriggerEvent = isPHI7 || isINT7 ;
1330 isMIXEvent = isINT7 ;
7deabbe7 1331 }
1332
83aecbbd 1333 if(fPeriod == kLHC11h)
7deabbe7 1334 {
08a327bf 1335 // Working: The same trigger in real and mixed events.
1336 isTriggerEvent = isCentral || isSemiCentral ;
1337 isMIXEvent = isCentral || isSemiCentral ;
1338 }
83aecbbd 1339
08a327bf 1340 // Select events.
1341 if(isTriggerEvent || isMIXEvent)
1342 {
1343 if ( isTriggerEvent )
83aecbbd 1344 {
83aecbbd 1345 FillHistogram("hTriggerPassedEvents", 17.);
08a327bf 1346 fPHOSEvent = true;
1347 }
1348
1349 if ( isMIXEvent )
1350 {
1351 FillHistogram("hTriggerPassedEvents", 18.);
1352 fMBEvent = true;
83aecbbd 1353 }
08a327bf 1354
1355 return ACCEPT;
67ef08bd 1356 }
7deabbe7 1357
1358 // other events
08a327bf 1359 FillHistogram("hTriggerPassedEvents", 19.);
1360
1361 //if ( isAny ) FillHistogram("hTriggerPassedEvents", 1.+20) ; // Not necessary
1362 // PHOS events.
1363 if ( isPHI1 ) FillHistogram("hTriggerPassedEvents", 2.+20 );
1364 if ( isPHI7 ) FillHistogram("hTriggerPassedEvents", 3.+20 );
1365 if ( isPHI8 ) FillHistogram("hTriggerPassedEvents", 4.+20 );
1366 if ( isCentral ) FillHistogram("hTriggerPassedEvents", 5.+20 );
1367 if ( isSemiCentral ) FillHistogram("hTriggerPassedEvents", 6.+20 );
1368 if ( isPHOSPb ) FillHistogram("hTriggerPassedEvents", 7.+20 );
1369 // MB events.
1370 if ( isMB ) FillHistogram("hTriggerPassedEvents", 8.+20 );
1371 if ( isINT7 ) FillHistogram("hTriggerPassedEvents", 9.+20 );
1372 if ( isAnyINT ) FillHistogram("hTriggerPassedEvents", 10.+20 );
1373 // Other rejected events
1374 if ( isMUON ) FillHistogram("hTriggerPassedEvents", 11.+20 );
1375 if ( isHighMult ) FillHistogram("hTriggerPassedEvents", 12.+20 );
1376 if ( isEMC1 ) FillHistogram("hTriggerPassedEvents", 13.+20 );
1377 if ( isCINT5 ) FillHistogram("hTriggerPassedEvents", 14.+20 );
1378 if ( isCMUS5 ) FillHistogram("hTriggerPassedEvents", 15.+20 );
1379 if ( isMUSPB ) FillHistogram("hTriggerPassedEvents", 16.+20 );
1380 if ( isMUSH7 ) FillHistogram("hTriggerPassedEvents", 17.+20 );
1381 if ( isMUSHPB ) FillHistogram("hTriggerPassedEvents", 18.+20 );
1382 if ( isMUL7 ) FillHistogram("hTriggerPassedEvents", 19.+20 );
1383 if ( isMuonLikePB ) FillHistogram("hTriggerPassedEvents", 20.+20 );
1384 if ( isMUU7 ) FillHistogram("hTriggerPassedEvents", 21.+20 );
1385 if ( isMuonUnlikePB ) FillHistogram("hTriggerPassedEvents", 22.+20 );
1386 if ( isEMC7 ) FillHistogram("hTriggerPassedEvents", 23.+20 );
1387 if ( isEMC8 ) FillHistogram("hTriggerPassedEvents", 24.+20 );
1388 if ( isMUS7 ) FillHistogram("hTriggerPassedEvents", 25.+20 );
1389 if ( isEMCEJE ) FillHistogram("hTriggerPassedEvents", 26.+20 );
1390 if ( isEMCEGA ) FillHistogram("hTriggerPassedEvents", 27.+20 );
1391 if ( isDG5 ) FillHistogram("hTriggerPassedEvents", 28.+20 );
1392 if ( isZED ) FillHistogram("hTriggerPassedEvents", 29.+20 );
1393 if ( isSPI7 ) FillHistogram("hTriggerPassedEvents", 30.+20 );
1394 if ( isSPI ) FillHistogram("hTriggerPassedEvents", 31.+20 );
1395 if ( isINT8 ) FillHistogram("hTriggerPassedEvents", 32.+20 );
1396 if ( isMuonSingleLowPt8 ) FillHistogram("hTriggerPassedEvents", 33.+20 );
1397 if ( isMuonSingleHighPt8 ) FillHistogram("hTriggerPassedEvents", 34.+20 );
1398 if ( isMuonLikeLowPt8 ) FillHistogram("hTriggerPassedEvents", 35.+20 );
1399 if ( isMuonUnlikeLowPt8 ) FillHistogram("hTriggerPassedEvents", 36.+20 );
1400 if ( isMuonUnlikeLowPt0 ) FillHistogram("hTriggerPassedEvents", 37.+20 );
1401 if ( isUserDefined ) FillHistogram("hTriggerPassedEvents", 38.+20 );
1402 if ( isTRD ) FillHistogram("hTriggerPassedEvents", 39.+20 );
1403 if ( isFastOnly ) FillHistogram("hTriggerPassedEvents", 40.+20 );
1404
7deabbe7 1405 return REJECT;
67ef08bd 1406}
42800acf 1407
67ef08bd 1408//_______________________________________________________________________________
1409void AliPHOSCorrelations::SetVertex()
1410{
1411 const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
1412 if( primaryVertex )
1413 {
1414 fVertex[0] = primaryVertex->GetX();
1415 fVertex[1] = primaryVertex->GetY();
1416 fVertex[2] = primaryVertex->GetZ();
1417 }
1418 else
1419 {
9aed5fcd 1420 //AliError("Event has 0x0 Primary Vertex, defaulting to origo");
67ef08bd 1421 fVertex[0] = 0;
1422 fVertex[1] = 0;
1423 fVertex[2] = 0;
1424 }
1425 fVertexVector = TVector3(fVertex);
1426
08a327bf 1427 fVtxBin=GetVertexBin(fVertexVector) ;
67ef08bd 1428}
42800acf 1429
67ef08bd 1430//_______________________________________________________________________________
08a327bf 1431Bool_t AliPHOSCorrelations::RejectEventVertex() const
67ef08bd 1432{
42800acf 1433 if( ! fEvent->GetPrimaryVertex() ) return true; // reject
67ef08bd 1434 LogSelection(kHasVertex, fInternalRunNumber);
1435
42800acf 1436 if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ ) return true; // reject
67ef08bd 1437 LogSelection(kHasAbsVertex, fInternalRunNumber);
1438
1439 return false; // accept event.
1440}
42800acf 1441
08a327bf 1442//_______________________________________________________________________________
1443void AliPHOSCorrelations::SetVertexBinning()
1444{
1445 // Define vertex bins by their edges
1446 const int nbins = fNVtxZBins+1;
1447 const double binWidth = 2*fMaxAbsVertexZ/fNVtxZBins;
1448 Double_t edges[nbins+1];
1449 for (int i = 0; i < nbins; ++i)
1450 {
1451 edges[i] = -1.*fNVtxZBins + binWidth*(double)i;
1452 }
1453
1454 TArrayD vtxEdges(nbins, edges);
1455
1456 for(int i=0; i<vtxEdges.GetSize()-1; ++i)
1457 {
1458 if(vtxEdges.At(i) > vtxEdges.At(i+1))
1459 AliFatal("edges are not sorted");
1460
1461 fVtxEdges = vtxEdges;
1462 }
1463}
1464
1465//_______________________________________________________________________________
1466Int_t AliPHOSCorrelations::GetVertexBin(const TVector3& vertexVector)
1467{
1468 int lastBinUpperIndex = fVtxEdges.GetSize() -1;
1469 if( vertexVector.z() > fVtxEdges[lastBinUpperIndex] )
1470 {
1471 if( fDebug >= 1 )
1472 AliWarning( Form("vertex (%f) larger then upper edge of last vertex bin (%f)!", vertexVector.z(), fVtxEdges[lastBinUpperIndex]) );
1473 return lastBinUpperIndex-1;
1474 }
1475 if( vertexVector.z() < fVtxEdges[0] )
1476 {
1477 if( fDebug >= 1 )
1478 AliWarning( Form("vertex (%f) smaller then lower edge of first bin (%f)!", vertexVector.z(), fVtxEdges[0]) );
1479 return 0;
1480 }
1481
1482 fVtxBin = TMath::BinarySearch<Double_t> ( GetNumberOfVertexBins(), fVtxEdges.GetArray(), vertexVector.z() );
1483 return fVtxBin;
1484}
1485
67ef08bd 1486//_______________________________________________________________________________
1487void AliPHOSCorrelations::SetCentrality()
1488{
1489 AliCentrality *centrality = fEvent->GetCentrality();
1490 if( centrality )
1491 fCentrality=centrality->GetCentralityPercentile(fCentralityEstimator);
1492 else
1493 {
1494 AliError("Event has 0x0 centrality");
1495 fCentrality = -1.;
1496 }
1497
67ef08bd 1498 fCentBin = GetCentralityBin(fCentrality);
1499}
42800acf 1500
67ef08bd 1501//_______________________________________________________________________________
08a327bf 1502Bool_t AliPHOSCorrelations::RejectEventCentrality() const
67ef08bd 1503{
42800acf 1504 if (fCentrality<fCentralityLowLimit)
1505 return true; //reject
1506 if(fCentrality>fCentralityHightLimit)
67ef08bd 1507 return true; //reject
67ef08bd 1508
1509 return false; // accept event.
1510}
42800acf 1511
67ef08bd 1512//_______________________________________________________________________________
42800acf 1513void AliPHOSCorrelations::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed)
1514{
1515 // Define centrality bins by their edges
1516 for(int i=0; i<edges.GetSize()-1; ++i)
1517 {
1518 if(edges.At(i) > edges.At(i+1)) AliFatal("edges are not sorted");
1519 if( edges.GetSize() != nMixed.GetSize()+1) AliFatal("edges and nMixed don't have appropriate relative sizes");
1520
1521 fCentEdges = edges;
1522 fCentNMixed = nMixed;
1523 }
67ef08bd 1524}
42800acf 1525
67ef08bd 1526//_______________________________________________________________________________
08a327bf 1527Int_t AliPHOSCorrelations::GetCentralityBin(const Float_t& centralityV0M)
42800acf 1528{
1529 int lastBinUpperIndex = fCentEdges.GetSize() -1;
1530 if( centralityV0M > fCentEdges[lastBinUpperIndex] )
1531 {
1532 if( fDebug >= 1 )
1533 AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1534 return lastBinUpperIndex-1;
1535 }
1536 if( centralityV0M < fCentEdges[0] )
1537 {
1538 if( fDebug >= 1 )
1539 AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1540 return 0;
1541 }
1542
1543 fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1544 return fCentBin;
67ef08bd 1545}
42800acf 1546
67ef08bd 1547//_______________________________________________________________________________
08a327bf 1548void AliPHOSCorrelations::SetCentralityBorders(const Double_t& downLimit , const Double_t& upLimit )
42800acf 1549{
1550 if (downLimit < 0. || upLimit > 100 || upLimit<=downLimit)
1551 AliError( Form("Warning. Bad value of centrality borders. Setting as default: fCentralityLowLimit=%2.f, fCentralityHightLimit=%2.f", fCentralityLowLimit, fCentralityHightLimit) );
1552 else
1553 {
1554 fCentralityLowLimit = downLimit;
1555 fCentralityHightLimit = upLimit;
1556 AliInfo( Form("Centrality border was set as fCentralityLowLimit=%2.f, fCentralityHightLimit=%2.f", fCentralityLowLimit, fCentralityHightLimit ) );
1557 }
67ef08bd 1558}
1559
1560//_______________________________________________________________________________
1561void AliPHOSCorrelations::EvalReactionPlane()
1562{
1563 // assigns: fHaveTPCRP and fRP
42800acf 1564 // also does RP histogram fill
67ef08bd 1565
1566 AliEventplane *eventPlane = fEvent->GetEventplane();
1567 if( ! eventPlane ) { AliError("Event has no event plane"); return; }
1568
1569 Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
1570
1571 if(reactionPlaneQ>=999 || reactionPlaneQ < 0.)
1572 {
1573 //reaction plain was not defined
1574 fHaveTPCRP = kFALSE;
1575 }
1576 else
1577 {
42800acf 1578 //reaction plain defined
67ef08bd 1579 fHaveTPCRP = kTRUE;
1580 }
1581
1582 if(fHaveTPCRP)
1583 fRP = reactionPlaneQ;
1584 else
1585 fRP = 0.;
1586
1587 FillHistogram("phiRPflat",fRP,fCentrality) ;
1588}
42800acf 1589
67ef08bd 1590//_______________________________________________________________________________
1591Int_t AliPHOSCorrelations::GetRPBin()
1592{
1593 Double_t averageRP;
42800acf 1594 averageRP = fRP ; // If possible, it is better to have EP bin from TPC
1595 // to have similar events for miximng (including jets etc) (fRPV0A+fRPV0C+fRP) /3.;
67ef08bd 1596 fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
42800acf 1597 if(fEMRPBin > (Int_t)fNEMRPBins-1)
67ef08bd 1598 fEMRPBin = fNEMRPBins-1 ;
1599 else
42800acf 1600 if(fEMRPBin < 0) fEMRPBin=0;
67ef08bd 1601
1602 return fEMRPBin;
1603}
42800acf 1604
67ef08bd 1605//_______________________________________________________________________________
1606void AliPHOSCorrelations::SelectPhotonClusters()
1607{
d2c19ce3 1608 //Selects PHOS clusters
67ef08bd 1609
d2c19ce3 1610 // clear (or create) array for holding events photons/clusters
1611 if(fCaloPhotonsPHOS)
1612 fCaloPhotonsPHOS->Clear();
9aed5fcd 1613 else
1614 {
d2c19ce3 1615 fCaloPhotonsPHOS = new TClonesArray("AliCaloPhoton",200);
9aed5fcd 1616 fCaloPhotonsPHOS->SetOwner();
d2c19ce3 1617 }
67ef08bd 1618
42800acf 1619 Int_t inPHOS = 0 ;
67ef08bd 1620
42800acf 1621 for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
9aed5fcd 1622 {
1623 AliVCluster *clu = fEvent->GetCaloCluster(i);
1624 if (!clu->IsPHOS() || clu->E()< fMinClusterEnergy) continue; // reject cluster
67ef08bd 1625
9aed5fcd 1626 Float_t position[3];
1627 clu->GetPosition(position);
1628 TVector3 global(position) ;
1629 Int_t relId[4] ;
1630 fPHOSGeo->GlobalPos2RelId(global,relId) ;
1631 Int_t modPHOS = relId[0] ;
1632 Int_t cellXPHOS = relId[2];
1633 Int_t cellZPHOS = relId[3] ;
1634
d2c19ce3 1635 Double_t distBC=clu->GetDistanceToBadChannel();
42800acf 1636 if(distBC<fMinBCDistance) continue ; // reject cluster
1637 if(clu->GetNCells() < fMinNCells) continue ; // reject cluster
1638 if(clu->GetM02() < fMinM02) continue ; // reject cluster
d2c19ce3 1639
42800acf 1640 if(fTOFCutEnabled)
1641 {
d2c19ce3 1642 Double_t tof = clu->GetTOF();
1643 if(TMath::Abs(tof) > fTOFCut ) continue ;
1644 }
1645 TLorentzVector lorentzMomentum;
9aed5fcd 1646 Double_t ecore = clu->GetCoreEnergy();
1647 //Double_t ecore = clu->E();
d2c19ce3 1648
42800acf 1649 FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
9aed5fcd 1650
1651 Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0) ?
9aed5fcd 1652 clu->GetMomentum(lorentzMomentum, origo);
1653
42800acf 1654 if(inPHOS>=fCaloPhotonsPHOS->GetSize())
d2c19ce3 1655 fCaloPhotonsPHOS->Expand(inPHOS+50) ;
d2c19ce3 1656
42800acf 1657 AliCaloPhoton * ph =new((*fCaloPhotonsPHOS)[inPHOS]) AliCaloPhoton(lorentzMomentum.X(), lorentzMomentum.Py(), lorentzMomentum.Z(), lorentzMomentum.E());
d2c19ce3 1658 inPHOS++ ;
1659 ph->SetCluster(clu);
1660
42800acf 1661 // Manual PHOS module number calculation
9aed5fcd 1662 /*Float_t cellId=clu->GetCellAbsId(0) ;
1663 Int_t mod = (Int_t)TMath:: Ceil(cellId/(56*64) ) ; */
1664 ph->SetModule(modPHOS) ;
1665
1666 lorentzMomentum*=ecore/lorentzMomentum.E() ;
d2c19ce3 1667
9aed5fcd 1668 //ph->SetNCells(clu->GetNCells());
1669 ph->SetMomV2(&lorentzMomentum) ;
42800acf 1670 ph->SetDispBit(clu->GetDispersion() < 2.5) ;
1671 ph->SetCPVBit(clu->GetEmcCpvDistance() > 2.) ;
d2c19ce3 1672
d2c19ce3 1673 FillHistogram(Form("QA_cluXZE_mod%i", modPHOS), cellXPHOS, cellZPHOS, lorentzMomentum.E() ) ;
1674 }
67ef08bd 1675}
42800acf 1676
67ef08bd 1677//_______________________________________________________________________________
1678void AliPHOSCorrelations::SelectAccosiatedTracks()
1679{
d2c19ce3 1680 // clear (or create) array for holding events tracks
67ef08bd 1681 if(fTracksTPC)
1682 fTracksTPC->Clear();
1683 else
1684 {
1685 fTracksTPC = new TClonesArray("TLorentzVector",12000);
1686 }
42800acf 1687 Int_t iTracks = 0 ;
1688 for (Int_t i = 0; i < fEvent->GetNumberOfTracks(); i++)
67ef08bd 1689 {
1690
08a327bf 1691 AliVTrack * track = (AliVTrack*)fEvent->GetTrack(i);
1692
1693 //Select tracks under certain conditions, TPCrefit, ITSrefit ...
1694 ULong_t status = track->GetStatus();
1695 if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1696 {
1697 if( fDebug > 2 ) printf("\t Reject track, status != fTrackStatus\n" );
1698 continue ;
1699 }
1700
42800acf 1701 if(fEventESD)
1702 {
08a327bf 1703 if(!SelectESDTrack((AliESDtrack*)track)) continue ; // reject track
67ef08bd 1704 }
42800acf 1705 else
1706 {
1707 if(!SelectAODTrack((AliAODTrack*)track)) continue ; // reject track
67ef08bd 1708 }
42800acf 1709
67ef08bd 1710 Double_t px = track->Px();
1711 Double_t py = track->Py();
42800acf 1712 Double_t pz = track->Pz();
1713 Double_t e = track->E() ;
67ef08bd 1714
42800acf 1715 if(iTracks >= fTracksTPC->GetSize())
1716 fTracksTPC->Expand(iTracks+50) ;
67ef08bd 1717
42800acf 1718 new((*fTracksTPC)[iTracks]) TLorentzVector(px, py, pz, e);
67ef08bd 1719 iTracks++ ;
1720 }
1721}
42800acf 1722
7deabbe7 1723//_______________________________________________________________________________
1724void AliPHOSCorrelations::SelectTriggerPi0ME()
1725{
42800acf 1726 const Int_t nPHOS = fCaloPhotonsPHOS->GetEntriesFast() ;
1727 for(Int_t i1 = 0; i1 < nPHOS-1; i1++)
7deabbe7 1728 {
42800acf 1729 AliCaloPhoton * ph1 = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1730 for (Int_t i2 = i1+1; i2 < nPHOS; i2++)
7deabbe7 1731 {
1732 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
42800acf 1733 TLorentzVector p12 = *ph1 + *ph2;
7deabbe7 1734
42800acf 1735 Double_t phiTrigger = p12.Phi() ;
1736 Double_t etaTrigger = p12.Eta() ;
7deabbe7 1737
42800acf 1738 Double_t m = p12.M() ;
1739 Double_t pt = p12.Pt();
7deabbe7 1740 Double_t eff = 1./GetEfficiency(pt);
1741 int mod1 = ph1->Module() ;
1742 int mod2 = ph2->Module() ;
1743
42800acf 1744 FillHistogram("clu_phieta", phiTrigger, etaTrigger );
1745 FillHistogram("clusingle_phieta", ph1->Phi(), ph1->Eta() );
1746 FillHistogram("clusingle_phieta", ph2->Phi(), ph2->Eta() );
7deabbe7 1747
42800acf 1748 FillHistogram("all_mpt", m, pt, eff );
83aecbbd 1749
7deabbe7 1750 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1751 {
42800acf 1752 FillHistogram("cpv_mpt", m, pt, eff );
83aecbbd 1753 }
7deabbe7 1754
1755 if ( ph1->IsDispOK() && ph2->IsDispOK() )
1756 {
42800acf 1757 FillHistogram("disp_mpt", m, pt, eff );
7deabbe7 1758 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1759 {
42800acf 1760 FillHistogram("both_mpt", m, pt, eff );
7deabbe7 1761 if(mod1 == mod2) // for each module
1762 {
42800acf 1763 FillHistogram(Form("both%d_mpt", mod1), m, pt, eff );
7deabbe7 1764 }
1765 }
1766 }
1767
42800acf 1768 if(!TestMass(m,pt)) continue; //reject this pair
7deabbe7 1769
1770 Int_t modCase = GetModCase(mod1, mod2);
1771
1772 //Now we choosing most energetic pi0.
1773 TestPi0ME(kPidAll, p12, modCase);
1774 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1775 TestPi0ME(kPidCPV, p12, modCase);
1776 if ( ph1->IsDispOK() && ph2->IsDispOK() )
1777 {
1778 TestPi0ME(kPidDisp, p12, modCase);
1779 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1780 TestPi0ME(kPidBoth, p12, modCase);
1781 }
1782 }
1783 }
1784}
1785
67ef08bd 1786//_______________________________________________________________________________
1787void AliPHOSCorrelations::ConsiderPi0s()
1788{
42800acf 1789 TString spid[4] = {"all","cpv","disp","both"} ;
7deabbe7 1790 // Counting number of trigger particles.
1791 for (int ipid = 0; ipid < 4; ipid++)
9aed5fcd 1792 {
7deabbe7 1793 if (fMEExists[ipid])
42800acf 1794 FillHistogram( Form("nTrigger_%s", spid[ipid].Data()), GetMEPt(ipid), 1./GetEfficiency(GetMEPt(ipid)) );
7deabbe7 1795 }
9aed5fcd 1796
7deabbe7 1797 // Take track's angles and compare with trigger's angles.
42800acf 1798 for(Int_t i3 = 0; i3 < fTracksTPC->GetEntriesFast(); i3++)
1799 {
7deabbe7 1800 TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i3);
9aed5fcd 1801
7deabbe7 1802 Double_t phiAssoc = track->Phi();
1803 Double_t etaAssoc = track->Eta();
42800acf 1804 Double_t ptAssoc = track->Pt() ;
9aed5fcd 1805
42800acf 1806 Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
7deabbe7 1807 Double_t dPhi(0.), dEta(0.);
9aed5fcd 1808
9aed5fcd 1809 for (int ipid = 0; ipid < 4; ipid++)
1810 {
7deabbe7 1811 if (GetMEExists(ipid))
9aed5fcd 1812 {
7deabbe7 1813 dPhi = GetMEPhi(ipid) - phiAssoc;
42800acf 1814 while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1815 while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
7deabbe7 1816 dEta = GetMEEta(ipid) - etaAssoc;
42800acf 1817 FillHistogram( Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), ptAssocBin), GetMEPt(ipid), dPhi, dEta, 1./GetEfficiency(GetMEPt(ipid)) );
7deabbe7 1818 }
1819 }
1820 }
9aed5fcd 1821}
42800acf 1822
9aed5fcd 1823//_______________________________________________________________________________
42800acf 1824void AliPHOSCorrelations::ConsiderPi0s_MBSelection()
83aecbbd 1825{
42800acf 1826 TString spid[4] = {"all","cpv","disp","both"} ;
83aecbbd 1827 // Counting number of trigger particles.
1828 for (int ipid = 0; ipid < 4; ipid++)
1829 {
1830 if (GetMEExists(ipid))
1831 {
1832
42800acf 1833 FillHistogram( Form("nTrigger_%s_MB", spid[ipid].Data()), GetMEPt(ipid), 1./GetEfficiency(GetMEPt(ipid)) );
83aecbbd 1834 }
1835 }
1836
1837 // Take track's angles and compare with trigger's angles.
42800acf 1838 for(Int_t i3 = 0; i3 < fTracksTPC->GetEntriesFast(); i3++)
1839 {
83aecbbd 1840 TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i3);
1841
1842 Double_t phiAssoc = track->Phi();
1843 Double_t etaAssoc = track->Eta();
42800acf 1844 Double_t ptAssoc = track->Pt();
83aecbbd 1845
42800acf 1846 Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
83aecbbd 1847 Double_t dPhi(0.), dEta(0.);
1848
1849 for (int ipid = 0; ipid < 4; ipid++)
1850 {
1851 if (GetMEExists(ipid))
1852 {
1853 dPhi = GetMEPhi(ipid) - phiAssoc;
42800acf 1854 while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1855 while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
83aecbbd 1856 dEta = GetMEEta(ipid) - etaAssoc;
1857 FillHistogram(Form("%s_ptphieta_ptAssoc_%3.1f_MB", spid[ipid].Data(), ptAssocBin), GetMEPt(ipid), dPhi, dEta, 1./GetEfficiency(GetMEPt(ipid)) );
1858 }
1859 }
1860 }
1861}
42800acf 1862
83aecbbd 1863//_______________________________________________________________________________
42800acf 1864void AliPHOSCorrelations::ConsiderPi0sMix()
9aed5fcd 1865{
42800acf 1866 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1867 for(Int_t evi = 0; evi < arrayList->GetEntries(); evi++)
1868 {
1869 TClonesArray * mixPHOS = static_cast<TClonesArray*>(arrayList->At(evi));
1870 for (Int_t i1 = 0; i1 < fCaloPhotonsPHOS->GetEntriesFast(); i1++)
1871 {
1872 AliCaloPhoton * ph1 = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1873 for(Int_t i2 = 0; i2 < mixPHOS->GetEntriesFast(); i2++)
1874 {
1875 AliCaloPhoton * ph2 = (AliCaloPhoton*)mixPHOS->At(i2) ;
1876 TLorentzVector p12 = *ph1 + *ph2;
1877 Double_t m = p12.M() ;
1878 Double_t pt = p12.Pt() ;
1879 Double_t eff = 1./GetEfficiency(pt);
1880
1881 int mod1 = ph1->Module() ;
1882 int mod2 = ph2->Module() ;
1883
1884 FillHistogram("mix_all_mpt", m, pt, eff);
1885 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1886 {
1887 FillHistogram("mix_cpv_mpt",m, pt, eff);
1888 }
1889 if ( ph1->IsDispOK() && ph2->IsDispOK() )
1890 {
1891 FillHistogram("mix_disp_mpt",m, pt, eff);
1892 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1893 {
1894 FillHistogram("mix_both_mpt",m, pt, eff);
1895 if (mod1 == mod2) // for each module
1896 {
1897 FillHistogram(Form("mix_both%d_mpt",mod1),m, pt, eff);
1898 }
1899 }
1900 }
1901 }
1902 }
1903 }
1904}
1905
1906//_______________________________________________________________________________
1907void AliPHOSCorrelations::ConsiderTracksMix()
1908{
1909 TString spid[4] = {"all","cpv","disp","both"} ;
9aed5fcd 1910
1911 TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
1912
42800acf 1913 for(Int_t evi = 0; evi < arrayList->GetEntries();evi++)
1914 {
9aed5fcd 1915 TClonesArray * mixTracks = static_cast<TClonesArray*>(arrayList->At(evi));
42800acf 1916 for(Int_t i3 = 0; i3 < mixTracks->GetEntriesFast(); i3++)
1917 {
9aed5fcd 1918 TLorentzVector * track = (TLorentzVector*)mixTracks->At(i3);
1919
1920 Double_t phiAssoc = track->Phi();
1921 Double_t etaAssoc = track->Eta();
08a327bf 1922 Double_t ptAssoc = track->Pt();
9aed5fcd 1923
42800acf 1924 Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
9aed5fcd 1925
42800acf 1926 Double_t ptTrigger(0.);
7deabbe7 1927
9aed5fcd 1928 Double_t dPhi(0.), dEta(0.);
1929
1930 for (int ipid = 0; ipid < 4; ipid++)
1931 {
7deabbe7 1932 if (GetMEExists(ipid))
9aed5fcd 1933 {
7deabbe7 1934 dPhi = GetMEPhi(ipid) - phiAssoc;
42800acf 1935 while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1936 while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
7deabbe7 1937 dEta = GetMEEta(ipid) - etaAssoc;
1938 ptTrigger = GetMEPt(ipid);
9aed5fcd 1939
7deabbe7 1940 FillHistogram(Form("mix_%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), ptAssocBin), ptTrigger, dPhi, dEta, 1./GetEfficiency(ptTrigger));
42800acf 1941 }
9aed5fcd 1942 }
1943 }
1944 }
1945}
1946
67ef08bd 1947//_______________________________________________________________________________
08a327bf 1948TList* AliPHOSCorrelations::GetCaloPhotonsPHOSList(const UInt_t vtxBin, const UInt_t centBin, const UInt_t rpBin)
42800acf 1949{
1950 int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1951 if( fCaloPhotonsPHOSLists->At(offset) )
1952 {
1953 // list exists
1954 TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1955 return list;
1956 }
1957 else
1958 {
1959 // no list for this bin has been created, yet
1960 TList* list = new TList();
1961 fCaloPhotonsPHOSLists->AddAt(list, offset);
1962 return list;
1963 }
67ef08bd 1964}
42800acf 1965
67ef08bd 1966//_______________________________________________________________________________
08a327bf 1967TList* AliPHOSCorrelations::GetTracksTPCList(const UInt_t vtxBin, const UInt_t centBin, const UInt_t rpBin)
1968{
42800acf 1969 int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1970 if( fTracksTPCLists->At(offset) )
1971 {
1972 // list exists
1973 TList* list = dynamic_cast<TList*> (fTracksTPCLists->At(offset));
1974 return list;
1975 }
1976 else
1977 {
1978 // no list for this bin has been created, yet
1979 TList* list = new TList();
1980 fTracksTPCLists->AddAt(list, offset);
1981 return list;
1982 }
67ef08bd 1983}
42800acf 1984
67ef08bd 1985//_______________________________________________________________________________
08a327bf 1986Double_t AliPHOSCorrelations::GetAssocBin(const Double_t& pt) const
9aed5fcd 1987{
42800acf 1988 //Calculates bin of associated particle pt.
1989 for(Int_t i=1; i<fAssocBins.GetSize(); i++)
1990 {
1991 if(pt>fAssocBins.At(i-1) && pt<fAssocBins.At(i))
1992 return fAssocBins.At(i) ;
1993 }
1994
1995 return fAssocBins.At(fAssocBins.GetSize()-1) ;
67ef08bd 1996}
42800acf 1997
67ef08bd 1998//_______________________________________________________________________________
08a327bf 1999void AliPHOSCorrelations::FillTrackEtaPhi() const
67ef08bd 2000{
d2c19ce3 2001 // Distribution TPC's tracks by angles.
42800acf 2002 for (Int_t i1=0; i1<fTracksTPC->GetEntriesFast(); i1++)
2003 {
d2c19ce3 2004 TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i1);
08a327bf 2005 Double_t phiAssoc = track->Phi();
2006 Double_t etaAssoc = track->Eta();
2007 Double_t ptAssoc = track->Pt() ;
2008
2009 FillHistogram( "track_phieta", phiAssoc, etaAssoc, ptAssoc );
d2c19ce3 2010 }
67ef08bd 2011}
2012
2013//_______________________________________________________________________________
2014void AliPHOSCorrelations::UpdatePhotonLists()
2015{
42800acf 2016 //Now we either add current events to stack or remove
2017 //If no photons in current event - no need to add it to mixed
2018
2019 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
2020 if( fDebug >= 2 )
2021 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
2022 if(fCaloPhotonsPHOS->GetEntriesFast()>0)
2023 {
2024 arrayList->AddFirst(fCaloPhotonsPHOS) ;
2025 fCaloPhotonsPHOS=0x0;
2026 if(arrayList->GetEntries() > fCentNMixed[fCentBin])
2027 {
2028 // Remove redundant events
2029 TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
2030 arrayList->RemoveLast() ;
2031 delete tmp;
2032 }
2033 }
67ef08bd 2034}
42800acf 2035
67ef08bd 2036//_______________________________________________________________________________
2037void AliPHOSCorrelations::UpdateTrackLists()
2038{
42800acf 2039 //Now we either add current events to stack or remove
2040 //If no photons in current event - no need to add it to mixed
2041
2042 TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
2043
2044 if( fDebug >= 2 )
2045 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
2046 if(fTracksTPC->GetEntriesFast()>0)
2047 {
2048 arrayList->AddFirst(fTracksTPC) ;
2049 fTracksTPC=0x0;
2050 if(arrayList->GetEntries() > fCentNMixed[fCentBin])
2051 {
2052 // Remove redundant events
2053 TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
2054 arrayList->RemoveLast() ;
2055 delete tmp;
2056 }
2057 }
67ef08bd 2058}
42800acf 2059
67ef08bd 2060//_______________________________________________________________________________
08a327bf 2061Bool_t AliPHOSCorrelations::SelectESDTrack(const AliESDtrack * t) const
67ef08bd 2062{
08a327bf 2063 // Estimate if this track can be used for correlation analisys. If all right - return "TRUE".
2064 Float_t pt = t->Pt();
2065 if( pt<0.5 || pt>20. ) return kFALSE ;
2066 if( TMath::Abs( t->Eta() ) > 0.8 ) return kFALSE;
2067 if( !fESDtrackCuts->AcceptTrack(t) ) return kFALSE ;
2068
2069 return kTRUE ;
67ef08bd 2070}
42800acf 2071
67ef08bd 2072//_______________________________________________________________________________
08a327bf 2073Bool_t AliPHOSCorrelations::SelectAODTrack(const AliAODTrack * t) const
67ef08bd 2074{
08a327bf 2075 // Estimate if this track can be used for correlation analisys. If all right - return "TRUE".
2076 Float_t pt = t->Pt() ;
2077 if( pt<0.5 || pt>20. ) return kFALSE ;
2078 if( TMath::Abs(t->Eta()) > 0.8 ) return kFALSE ;
42800acf 2079
08a327bf 2080 if(fSelectHybridTracks)
67ef08bd 2081 {
08a327bf 2082 if(!t->IsHybridGlobalConstrainedGlobal())
2083 {
2084 if( fDebug > 2 ) printf("\t Reject track, IsHybridGlobalConstrainedGlobal() is FALSE\n ");
2085 return kFALSE ;
2086 }
67ef08bd 2087 }
08a327bf 2088
2089 if(fSelectSPDHitTracks) //Not much sense to use with TPC only or Hybrid tracks
2090 {
2091 if(!t->HasPointOnITSLayer(0) && !t->HasPointOnITSLayer(1))
2092 {
2093 if( fDebug > 2 ) printf("\t Reject track, HasPointOnITSLayer(0) is %i and HasPointOnITSLayer(1) is %i\n", t->HasPointOnITSLayer(0), t->HasPointOnITSLayer(1) );
2094 return kFALSE ; ;
2095 }
2096 }
2097
2098 if(fSelectFractionTPCSharedClusters)
2099 {
2100 Double_t frac = Double_t(t->GetTPCnclsS()) / Double_t(t->GetTPCncls());
2101 if (frac > fCutTPCSharedClustersFraction)
2102 {
2103 if( fDebug > 2 ) printf("\t Reject track, shared cluster fraction %f > %f\n",frac, fCutTPCSharedClustersFraction);
2104 return kFALSE ;
2105 }
2106 }
67ef08bd 2107
2108 return kTRUE ;
2109}
42800acf 2110
67ef08bd 2111//_______________________________________________________________________________
2112void AliPHOSCorrelations::LogProgress(int step)
67ef08bd 2113{
42800acf 2114 // Fill "step by step" hist
2115 FillHistogram("hTotSelEvents", step+0.5);
67ef08bd 2116}
42800acf 2117
67ef08bd 2118//_______________________________________________________________________________
08a327bf 2119void AliPHOSCorrelations::LogSelection(const int& step, const int& internalRunNumber) const
67ef08bd 2120{
42800acf 2121 // the +0.5 is not realy neccisarry, but oh well... -henrik
2122 FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
08a327bf 2123}
42800acf 2124
67ef08bd 2125//_______________________________________________________________________________
08a327bf 2126Bool_t AliPHOSCorrelations::TestMass(const Double_t& m, const Double_t& pt) const
d2c19ce3 2127{
08a327bf 2128 // If pi0 candidate outside of mass peak, then reject.
2129 if (fUseMassWindowParametrisation && fNSigmaWidth)
9aed5fcd 2130 {
08a327bf 2131 // Parametrization
2132 FillHistogram("massWindow", MassMeanFunction(pt), MassSigmaFunction(pt)*fNSigmaWidth);
2133 if ( MassMeanFunction(pt)-MassSigmaFunction(pt)*fNSigmaWidth<m && m<MassMeanFunction(pt)+MassSigmaFunction(pt)*fNSigmaWidth )
7deabbe7 2134 {
2135 FillHistogram("massWindowPass", 1);
2136 return true;
2137 }
2138 else
2139 {
2140 FillHistogram("massWindowPass", 2);
2141 return false;
2142 }
08a327bf 2143 }
42800acf 2144 else
08a327bf 2145 {
2146 if(!fNSigmaWidth)
2147 {
2148 AliInfo( Form("fNSigmaWidth equel 0. Class will use default mass window: from %f to %f GeV.", fMassInvMeanMin, fMassInvMeanMax ) );
2149 }
2150
2151 FillHistogram("massWindow", fMassInvMeanMin, fMassInvMeanMax);
2152 if(fMassInvMeanMin>m && m<fMassInvMeanMax)
7deabbe7 2153 {
2154 FillHistogram("massWindowPass", 3);
2155 return true;
2156 }
2157 else
2158 {
2159 FillHistogram("massWindowPass", 4);
2160 return false;
2161 }
08a327bf 2162 }
67ef08bd 2163}
42800acf 2164
67ef08bd 2165//_______________________________________________________________________________
08a327bf 2166Double_t AliPHOSCorrelations::MassMeanFunction(const Double_t &pt) const
9aed5fcd 2167{
2168 // Parametrization mean of mass window
42800acf 2169 return ( fMassMean[0]*pt + fMassMean[1] );
9aed5fcd 2170}
42800acf 2171
9aed5fcd 2172//_______________________________________________________________________________
08a327bf 2173Double_t AliPHOSCorrelations::MassSigmaFunction(const Double_t &pt) const
67ef08bd 2174{
9aed5fcd 2175 // Parametrization sigma of mass window
08a327bf 2176 return ( TMath::Sqrt(fMassSigma[0]*fMassSigma[0]/pt + fMassSigma[1]*fMassSigma[1]/pt/pt + fMassSigma[2]*fMassSigma[2]));
9aed5fcd 2177}
42800acf 2178
9aed5fcd 2179//_____________________________________________________________________________
08a327bf 2180void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x) const
42800acf 2181{
2182 //FillHistogram
2183 TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
2184 if(hist)
2185 hist->Fill(x) ;
2186 else
2187 AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
67ef08bd 2188}
42800acf 2189
9aed5fcd 2190//_____________________________________________________________________________
08a327bf 2191void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x, Double_t y)const
42800acf 2192{
2193 //FillHistogram
2194 TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
2195 if(th1)
2196 th1->Fill(x, y) ;
2197 else
2198 AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
9aed5fcd 2199}
2200
2201//_____________________________________________________________________________
08a327bf 2202void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x, Double_t y, Double_t z) const
42800acf 2203{
2204 //Fills 1D histograms with key
2205 TObject * obj = fOutputContainer->FindObject(key);
2206
2207 TH2 * th2 = dynamic_cast<TH2*> (obj);
2208 if(th2)
2209 {
2210 th2->Fill(x, y, z) ;
2211 return;
2212 }
2213
2214 TH3 * th3 = dynamic_cast<TH3*> (obj);
2215 if(th3)
2216 {
2217 th3->Fill(x, y, z) ;
2218 return;
2219 }
2220
2221 AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
67ef08bd 2222}
42800acf 2223
9aed5fcd 2224//_____________________________________________________________________________
08a327bf 2225void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x,Double_t y, Double_t z, Double_t w) const
42800acf 2226{
2227 //Fills 1D histograms with key
2228 TObject * obj = fOutputContainer->FindObject(key);
2229
2230 TH3 * th3 = dynamic_cast<TH3*> (obj);
2231 if(th3)
2232 {
2233 th3->Fill(x, y, z, w) ;
2234 return;
2235 }
2236
2237 AliError(Form("can not find histogram (of instance TH3) <%s> ",key)) ;
67ef08bd 2238}
42800acf 2239
0acaf361 2240//_____________________________________________________________________________
2241void AliPHOSCorrelations::SetGeometry()
2242{
42800acf 2243 // Initialize the PHOS geometry
2244 //Init geometry
2245 if(!fPHOSGeo)
2246 {
2247 AliOADBContainer geomContainer("phosGeo");
2248 geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
2249 TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
2250 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
2251 for(Int_t mod=0; mod<5; mod++)
2252 {
2253 if(!matrixes->At(mod))
2254 {
2255 if( fDebug )
2256 AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2257 continue;
2258 }
2259 else
2260 {
2261 fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
2262 if( fDebug >1 )
2263 AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2264 }
2265 }
2266 }
0acaf361 2267}
42800acf 2268
9aed5fcd 2269//_____________________________________________________________________________
08a327bf 2270Double_t AliPHOSCorrelations::GetEfficiency(const Double_t& pt) const
42800acf 2271{
08a327bf 2272 Double_t x = pt;
9aed5fcd 2273 //Efficiency for Both2core only!
83aecbbd 2274 if (!fUseEfficiency)
2275 return 1.;
9aed5fcd 2276
2277 Double_t e =1.;
2278 // From 0 to 5 - 11h for different centrality.
2279 /*0: 0-5%
2280 1: 5-10%
2281 2: 10-20%
2282 3: 20-40%
2283 4: 40-60%
2284 5: 60-80%
2285 6: 0-20%
2286 7: 0-10%*/
42800acf 2287 Double_t par0[9] = {-798863, 339.714, 6407.1, -457.778, 1283.65, -117.075, -19.3764, 0, 0 };
2288 Double_t par1[9] = {-799344, -1852.1, 3326.29, -384.229, 504.046, 562.608, 130.518, 0, 0 };
2289 Double_t par2[9] = {-858904, -1923.28, 5350.74, -568.946, 945.497, 419.647, 101.911, 0, 0 };
2290 Double_t par3[9] = {-795652, -1495.97, 2926.46, -357.804, 478.961, 551.127, 128.86, 0, 0 };
2291 Double_t par4[9] = {-891951, 279626, -123110, -5464.75, 27470.8, 283264, 15355.1, 192762, 44828.6 };
2292 Double_t par5[9] = {-1.1094e+06, -986.915, 2127.71, -268.908, 375.594, 380.791, 89.4053, 0, 0 };
9aed5fcd 2293 // Double_t par6[7] = {4.86106e+09, 4.47013e+08, -1.48079e+09, 1.47233e+08, -2.62356e+08, -1.00639e+08, -2.45629e+07, 0, 0};
2294 // Double_t par7[7] = {-1.36243e+06, -26011.1, 135838, -12161.3, 24956.8, 4985.4, 1285.57, 0, 0};
2295
42800acf 2296 // 8 bin for pPb13 and 0-100%
2297 Double_t par8[9] = {6.87095e+06, 8.36553e+06, -3.29572e+06, 2.18688e+06, -739490, 521666, 106661, 0, 0};
9aed5fcd 2298
9aed5fcd 2299 Double_t* pFitPoint;
2300
08a327bf 2301 // TODO: fix functions value below 1 GeV/c.
2302 if(x < 1.) x = 1.;
2303
9aed5fcd 2304 if(fPeriod == kLHC11h)
2305 {
42800acf 2306 if (fCentrality <= 5) pFitPoint = &par0[0];
2307 if (fCentrality > 5 && fCentrality <= 10) pFitPoint = &par1[0];
2308 if (fCentrality > 10 && fCentrality <= 20) pFitPoint = &par2[0];
2309 if (fCentrality > 20 && fCentrality <= 40) pFitPoint = &par3[0];
2310 if (fCentrality > 40 && fCentrality <= 60) pFitPoint = &par4[0];
2311 if (fCentrality > 60) pFitPoint = &par5[0];
9aed5fcd 2312
2313 Double_t pFit[9];
7deabbe7 2314 for (int i = 0; i < 9; ++i)
9aed5fcd 2315 {
2316 pFit[i] = *(pFitPoint+i);
2317 }
2318
42800acf 2319 if (fCentrality > 40 && fCentrality <= 60)
9aed5fcd 2320 e = TMath::Exp(-(((((1.+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))+(pFit[7]*(x*(x*(x*x)))))/((((pFit[3]*x)+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x))))+(pFit[8]*(x*(x*(x*x))))))) ;
2321 else
2322 e = TMath::Exp(-((((1.+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))/(((pFit[3]*x)+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x)))))) ;
2323 }
2324 else
2325 if( fPeriod == kLHC13 )
2326 {
2327 pFitPoint = &par8[0];
2328 Double_t pFit[9];
1207640c 2329 for( int i = 0; i < 9; i++ )
9aed5fcd 2330 {
2331 pFit[i] = *(pFitPoint+i);
2332 }
2333
2334 e = TMath::Exp(-((((pFit[0]+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))/(((1.+(pFit[3]*x))+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x)))))) ;
2335 }
2336 else
2337 {
2338 // No case
83aecbbd 2339 AliWarning(Form("No efficiensy choise. Return 1"));
9aed5fcd 2340 e = 1.;
2341 }
2342
2343 return e;
4ae540f6 2344 // return 1.; // For test.
9aed5fcd 2345}
9aed5fcd 2346
42800acf 2347//_____________________________________________________________________________
08a327bf 2348Int_t AliPHOSCorrelations::GetModCase(const Int_t &mod1, const Int_t &mod2) const
42800acf 2349{
9aed5fcd 2350 // Return modules pair namber.
2351 if(mod1 == mod2)
2352 {
2353 if(mod1 == 1) return 1;
2354 if(mod1 == 2) return 2;
2355 if(mod1 == 3) return 3;
2356 }
2357 else
2358 {
2359 if(mod1 == 1 || mod2 == 1)
2360 if(mod1 == 2 || mod2 == 2)
2361 return 12;
2362
2363 if(mod1 == 1 || mod2 == 1)
2364 if(mod1 == 3 || mod2 == 3)
2365 return 13;
2366 if(mod1 == 2 || mod2 == 2)
2367 if(mod1 == 3 || mod2 == 3)
2368 return 23;
2369 }
2370
2371 AliError(Form("No choise for mod1 = %i, mod2 = %i", mod1, mod2));
2372 return 1;
2373}
42800acf 2374
9aed5fcd 2375//_____________________________________________________________________________
08a327bf 2376void AliPHOSCorrelations::TestPi0ME(const Int_t& ipid, const TLorentzVector& p12, const Int_t& modCase)
9aed5fcd 2377{
42800acf 2378 Double_t phiTrigger = p12.Phi() ;
2379 Double_t etaTrigger = p12.Eta() ;
2380 Double_t pt = p12.Pt() ;
9aed5fcd 2381
7deabbe7 2382 if ( GetMEExists(ipid) )
2383 {
42800acf 2384 if ( pt >= GetMEPt(ipid) )
7deabbe7 2385 {
2386 SetMEPt(ipid,pt);
2387 SetMEPhi(ipid, phiTrigger);
2388 SetMEEta(ipid, etaTrigger);
2389 SetMEModCase(ipid, modCase);
2390 }
2391 }
2392 else
2393 {
2394 SetMEPt(ipid,pt);
2395 SetMEPhi(ipid, phiTrigger);
2396 SetMEEta(ipid, etaTrigger);
2397 SetMEModCase(ipid, modCase);
2398 SetMEExists(ipid);
2399 }
9aed5fcd 2400}
42800acf 2401
9aed5fcd 2402//_____________________________________________________________________________
42800acf 2403void AliPHOSCorrelations::ZeroingVariables()
2404{
7deabbe7 2405 // Set Phi, Eta, pT, modNumber andtrigger variable of moust energetic trigger particle to zero.
9aed5fcd 2406 for (int i = 0; i < 4; ++i)
2407 {
7deabbe7 2408 fMEExists[i] = false;
9aed5fcd 2409 fMEPhi[i] = fMEEta[i] = fMEPt[i] = -99;
2410 fMEModCase[i] = 1;
9aed5fcd 2411 }
2412}
2413
42800acf 2414//_____________________________________________________________________________
08a327bf 2415void AliPHOSCorrelations::SetMassMeanParametrs(const Double_t par[2])
42800acf 2416{
2417 for (int i = 0; i < 2; ++i)
2418 {
2419 fMassMean[i] = par[i] ;
2420 }
2421}
2422
2423//_____________________________________________________________________________
08a327bf 2424void AliPHOSCorrelations::SetMassSigmaParametrs(const Double_t par[3])
42800acf 2425{
08a327bf 2426 for (int i = 0; i < 3; ++i)
42800acf 2427 {
2428 fMassSigma[i] = par[i] ;
2429 }
2430}
08a327bf 2431
2432//_____________________________________________________________________________
2433void AliPHOSCorrelations::FillEventBiningProperties() const
2434{
2435 // Fill fCentBin, fEMRPBin, fVtxBin.
2436 if(fPHOSEvent || fMBEvent)
2437 {
2438 FillHistogram( "hCentralityBining", fCentBin) ;
2439 FillHistogram( "phiRPflatBining", fEMRPBin) ;
2440 FillHistogram( "hVertexZBining", fVtxBin) ;
2441 if(fPHOSEvent)
2442 {
2443 FillHistogram( "hCentralityBiningTrigger", fCentBin) ;
2444 FillHistogram( "phiRPflatBiningTrigger", fEMRPBin) ;
2445 FillHistogram( "hVertexZBiningTrigger", fVtxBin) ;
2446 }
2447 if(fMBEvent)
2448 {
2449 FillHistogram( "hCentralityBiningMB", fCentBin) ;
2450 FillHistogram( "phiRPflatBiningMB", fEMRPBin) ;
2451 FillHistogram( "hVertexZBiningMB", fVtxBin) ;
2452 }
2453 }
2454}
2455