1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "TObjArray.h"
25 #include "TParticle.h"
29 #include "THashList.h"
31 #include "AliAnalysisManager.h"
32 #include "AliMCEventHandler.h"
33 #include "AliMCEvent.h"
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisTaskPi0Flow.h"
37 #include "AliCaloPhoton.h"
38 #include "AliPHOSGeometry.h"
39 #include "TGeoManager.h"
40 #include "AliPHOSEsdCluster.h"
41 #include "AliPHOSCalibData.h"
42 #include "AliESDEvent.h"
43 #include "AliESDCaloCells.h"
44 #include "AliESDVertex.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliAODEvent.h"
49 #include "AliCDBManager.h"
50 #include <AliAODCaloCluster.h>
51 #include "AliCentrality.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliEventplane.h"
55 #include "AliOADBContainer.h"
56 #include "AliEPFlattener.h"
58 // Analysis task to fill histograms with PHOS ESD or AOD clusters and cells
59 // Authors : Dmitri Peressounko
61 // Modified: 03.08.2012 Henrik Qvigstad
64 ClassImp(AliAnalysisTaskPi0Flow);
66 //________________________________________________________________________
67 Double_t rnlin(Double_t *x, Double_t * /*par*/)
69 //a = par[0], b = par[1].
72 // return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
73 return 1.015*(0.0241+1.0504*x[0]+0.000249*x[0]*x[0]) ;
77 //________________________________________________________________________
78 AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
79 : AliAnalysisTaskSE(name),
85 fManualV0EPCalc(false),
86 fOutputContainer(0x0),
93 fInternalRunNumber(0),
96 fV0Cpol(0.),fV0Apol(0.),
99 fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"),
112 fCaloPhotonsPHOS(0x0),
113 fCaloPhotonsPHOSLists(0x0)
116 Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
117 TArrayD centEdges(nbins+1, edges);
118 Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
119 TArrayI centNMixed(nbins, nMixed);
120 SetCentralityBinning(centEdges, centNMixed);
122 for(int mod=1; mod <= 5; ++mod)
123 fModuleEnabled[mod-1] = kFALSE;
125 for(Int_t i=0;i<kNCenBins;i++){
126 for(Int_t j=0;j<2; j++)
127 for(Int_t k=0; k<2; k++) {
133 fVertex[0]=0; fVertex[1]=0; fVertex[2]=0;
135 // Output slots #0 write into a TH1 container
136 DefineOutput(1,TList::Class());
140 // Set bad channel map
142 for(Int_t i=0; i<6; i++){
143 snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
144 fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
148 // Initialize non-linrarity correction
149 fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
153 //___________________________________________________________________________
154 AliAnalysisTaskPi0Flow::~AliAnalysisTaskPi0Flow()
157 delete fESDtrackCuts;
158 delete fPHOSCalibData;
159 delete fCaloPhotonsPHOSLists;
160 if(fTPCFlat)delete fTPCFlat; fTPCFlat=0x0;
161 if(fV0AFlat)delete fV0AFlat; fV0AFlat=0x0;
162 if(fV0CFlat)delete fV0CFlat; fV0CFlat=0x0;
165 //________________________________________________________________________
166 void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
170 const Int_t nRuns=200 ;
173 if(fOutputContainer != NULL){
174 delete fOutputContainer;
176 fOutputContainer = new THashList();
177 fOutputContainer->SetOwner(kTRUE);
179 //========QA histograms=======
182 fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 12,0.,13.,nRuns,0.,float(nRuns))) ;
183 fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 12,0.,12.)) ;
185 //vertex distribution
186 fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
189 fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
190 fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
191 fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
192 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
193 fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,10.,100,0.,100.)) ;
197 fOutputContainer->Add(new TH3F("hPHOSphi","cos" ,10,0.,100.,20,0.,10.,100,-TMath::Pi(),TMath::Pi()));
199 fOutputContainer->Add(new TH2F("cos2AC","RP correlation between TPC subs", 100,-1.,1.,20,0.,100.)) ;
200 fOutputContainer->Add(new TH2F("cos2V0AC","RP correlation between VO A and C sides", 100,-1.,1.,20,0.,100.)) ;
201 fOutputContainer->Add(new TH2F("cos2V0ATPC","RP correlation between TPC and V0A", 100,-1.,1.,20,0.,100.)) ;
202 fOutputContainer->Add(new TH2F("cos2V0CTPC","RP correlation between TPC and V0C", 100,-1.,1.,20,0.,100.)) ;
204 fOutputContainer->Add(new TH2F("phiRP","RP distribution with TPC", 100,0.,TMath::Pi(),20,0.,100.)) ;
205 fOutputContainer->Add(new TH2F("phiRPflat","RP distribution with TPC flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
206 fOutputContainer->Add(new TH2F("phiRPV0A","RP distribution with V0A", 100,0.,TMath::Pi(),20,0.,100.)) ;
207 fOutputContainer->Add(new TH2F("phiRPV0C","RP distribution with V0C", 100,0.,TMath::Pi(),20,0.,100.)) ;
208 fOutputContainer->Add(new TH3F("phiRPV0AC","RP distribution with V0A and V0C", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
209 fOutputContainer->Add(new TH2F("phiRPV0Aflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
210 fOutputContainer->Add(new TH2F("phiRPV0Cflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
211 fOutputContainer->Add(new TH3F("phiRPV0ATPC","RP distribution with V0A + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
212 fOutputContainer->Add(new TH3F("phiRPV0CTPC","RP distribution with V0C + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
216 fOutputContainer->Add(new TH1I("hCellMultEvent" ,"PHOS cell multiplicity per event" ,2000,0,2000));
217 fOutputContainer->Add(new TH1I("hCellMultEventM1","PHOS cell multiplicity per event, M1",2000,0,2000));
218 fOutputContainer->Add(new TH1I("hCellMultEventM2","PHOS cell multiplicity per event, M2",2000,0,2000));
219 fOutputContainer->Add(new TH1I("hCellMultEventM3","PHOS cell multiplicity per event, M3",2000,0,2000));
221 fOutputContainer->Add(new TH1F("hCellEnergy" ,"Cell energy" ,3000,0.,30.));
222 fOutputContainer->Add(new TH1F("hCellEnergyM1","Cell energy in module 1",3000,0.,30.));
223 fOutputContainer->Add(new TH1F("hCellEnergyM2","Cell energy in module 2",3000,0.,30.));
224 fOutputContainer->Add(new TH1F("hCellEnergyM3","Cell energy in module 3",3000,0.,30.));
226 fOutputContainer->Add(new TH2F("hCellNXZM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
227 fOutputContainer->Add(new TH2F("hCellNXZM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
228 fOutputContainer->Add(new TH2F("hCellNXZM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
229 fOutputContainer->Add(new TH2F("hCellEXZM1","Cell E(X,Z), M1",64,0.5,64.5, 56,0.5,56.5));
230 fOutputContainer->Add(new TH2F("hCellEXZM2","Cell E(X,Z), M2",64,0.5,64.5, 56,0.5,56.5));
231 fOutputContainer->Add(new TH2F("hCellEXZM3","Cell E(X,Z), M3",64,0.5,64.5, 56,0.5,56.5));
234 fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
235 fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
236 fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
238 fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
239 fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
240 fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
242 fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
243 fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
244 fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
246 fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
247 fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
248 fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
251 //Single photon and pi0 spectrum
252 const Int_t nPtPhot = 400 ;
253 const Double_t ptPhotMax = 40 ;
254 const Int_t nM = 500;
255 const Double_t mMin = 0.0;
256 const Double_t mMax = 1.0;
258 //PHOS calibration QA
259 fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
260 fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
261 fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
262 fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
263 fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
264 fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
266 // Histograms for different centralities
267 const int kNPID = 16;
268 const char* pidNames[kNPID] = {"All", "Allcore", "Allwou", "Disp", "Disp2", "Dispcore", "Disp2core", "Dispwou", "CPV", "CPVcore", "CPV2", "CPV2core", "Both", "Bothcore", "Both2", "Both2core"};
271 for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
272 for(Int_t ipid=0; ipid < kNPID; ipid++){
273 name = Form("hPhot%s_cen%i", pidNames[ipid], cent );
274 title = Form("%s clusters", pidNames[ipid]);
275 fOutputContainer->Add(new TH1F(name.Data(), title.Data(), nPtPhot,0.,ptPhotMax));
277 name = Form("hPi0%s_cen%i", pidNames[ipid], cent );
278 title = Form("%s clusters", pidNames[ipid]);
279 fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
281 name = Form("hSingle%s_cen%i", pidNames[ipid], cent );
282 title = Form("%s clusters", pidNames[ipid]);
283 fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
285 name = Form("hMiPi0%s_cen%i", pidNames[ipid], cent );
286 title = Form("%s clusters", pidNames[ipid]);
287 fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
289 name = Form("hMiSingle%s_cen%i", pidNames[ipid], cent );
290 title = Form("%s clusters", pidNames[ipid]);
291 fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
293 // PhotPhi histograms
294 const Int_t nPt = 20;
295 const Double_t xPt[21]={0.6,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.} ;
296 const Int_t nPhi=10 ;
297 Double_t xPhi[nPhi+1] ;
298 for(Int_t i=0;i<=nPhi;i++)
299 xPhi[i]=i*TMath::Pi() /nPhi ;
300 const Int_t nMm=200 ;
302 for(Int_t i=0;i<=nMm;i++)
304 const Int_t kNPhiTitles = 3;
305 const char* phiTitles[kNPhiTitles] = {"TPC", "V0A", "V0C"};
306 for(Int_t iRP=0; iRP<3; iRP++){
307 name = Form("hPhotPhi%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
308 title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
309 fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nPt,xPt,nPhi,xPhi));
311 name = Form("hMassPt%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
312 title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
313 fOutputContainer->Add(new TH3F(name.Data(), title.Data(), nMm,xM,nPt,xPt,nPhi,xPhi));
315 name = Form("hMiMassPt%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
316 title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
317 fOutputContainer->Add(new TH3F(name.Data(), title.Data(), nMm,xM,nPt,xPt,nPhi,xPhi));
322 snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
323 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
324 snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
325 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
326 snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
327 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
328 snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
329 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
330 snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
331 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
333 snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
334 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
335 snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
336 fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
337 snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
338 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
339 snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
340 fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
341 snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
342 fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
345 // Setup photon lists
346 Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
347 fCaloPhotonsPHOSLists = new TObjArray(kapacity);
348 fCaloPhotonsPHOSLists->SetOwner();
350 PostData(1, fOutputContainer);
353 //________________________________________________________________________
354 void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
356 // Main loop, called for each event
360 // Step 0: Event Objects
362 fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
363 fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
364 fMCStack = GetMCStack();
368 // Step 1: Run Number, Misalignment Matrix, and Calibration
369 // fRunNumber, fInternalRunNumber, fMultV0, fV0Cpol, fV0Apol, fMeanQ, fWidthQ
370 if( fRunNumber != fEvent->GetRunNumber()) { // Check run number
371 // this should run only at first call of UserExec(),
372 // or if task runs over multiple runs, which should not occur in normal use.
374 // if run number has changed, set run variables
375 fRunNumber = fEvent->GetRunNumber();
376 fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
377 // then set misalignment and V0 calibration
386 LogSelection(0, fInternalRunNumber);
390 // fVertex, fVertexVector, fVtxBin
392 if( RejectEventVertex() ) {
393 PostData(1, fOutputContainer);
399 // if(event->IsPileupFromSPD()){
400 // PostData(1, fOutputContainer);
401 // return; // Reject!
406 // Step 4: Centrality
407 // fCentralityV0M, fCentBin
409 if( RejectCentrality() ){
410 PostData(1, fOutputContainer);
416 // Step 5: Reaction Plane
417 // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
418 EvalReactionPlane(); //TODO: uncomment this, or at least deal with it
419 EvalV0ReactionPlane(); //TODO: uncomment this, or at least deal with it
420 fEMRPBin = GetRPBin(); //TODO: uncomment this, or at least deal with it
429 // Step 7: QA PHOS cells
430 FillPHOSCellQAHists();
434 // Step 8: Event Photons (PHOS Clusters) selection
435 SelectPhotonClusters();
436 FillSelectedClusterHistograms();
439 if( ! fCaloPhotonsPHOS->GetEntriesFast() )
442 LogSelection(6, fInternalRunNumber);
445 // Step 9: Consider pi0 (photon/cluster) pairs.
453 // Step 11: Update lists
459 PostData(1, fOutputContainer);
462 //________________________________________________________________________
463 // void AliAnalysisTaskPi0Flow::Terminate(Option_t *)
465 // // Draw result to the screen
466 // // Called once at the end of the query
468 // // TH1 * hTotSelEvents = dynamic_cast<TH1*>(fOutputContainer->FindObject("hTotSelEvents"));
469 // // hTotSelEvents->Draw();
471 //________________________________________________________________________
472 void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed)
474 // Define centrality bins by their edges
475 if( edges.At(0) < 0.) AliFatal("lower edge less then 0");
476 if( 90. < edges.At(edges.GetSize()-1) ) AliFatal("upper edge larger then 90.");
477 for(int i=0; i<edges.GetSize()-1; ++i)
478 if(edges.At(i) > edges.At(i+1)) AliFatal("edges are not sorted");
479 if( edges.GetSize() != nMixed.GetSize()+1) AliFatal("edges and nMixed don't have appropriate relative sizes");
482 fCentNMixed = nMixed;
485 //_____________________________________________________________________________
486 void AliAnalysisTaskPi0Flow::SetEnablePHOSModule(int module, Bool_t enable)
488 if( module < 1 || 5 < module )
489 AliFatal("PHOS Module must be between 1 and 5");
490 fModuleEnabled[module-1] = enable;
494 //________________________________________________________________________
495 void AliAnalysisTaskPi0Flow::SetPHOSBadMap(Int_t mod, TH2I* badMapHist)
498 delete fPHOSBadMap[mod];
500 fPHOSBadMap[mod]=new TH2I(*badMapHist);
502 AliInfo(Form("Setting Bad Map Histogram %s",fPHOSBadMap[mod]->GetName()));
505 //________________________________________________________________________
506 Bool_t AliAnalysisTaskPi0Flow::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
508 //Check if this channel belogs to the good ones
510 if(strcmp(det,"PHOS")==0){
512 AliError(Form("No bad map for PHOS module %d",mod)) ;
515 if(!fPHOSBadMap[mod]){
516 AliError(Form("No Bad map for PHOS module %d, !fPHOSBadMap[mod]",mod)) ;
519 if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
525 AliError(Form("Can not find bad channels for detector %s ",det)) ;
528 //Remove 6 noisy channels in run 139036, LHC10h
529 if( 139036 == fRunNumber
531 && (ix==9||ix==10||ix==11)
532 && (iz==45 || iz==46))
537 //_____________________________________________________________________________
538 void AliAnalysisTaskPi0Flow::FillPHOSCellQAHists()
540 // Fill cell occupancy per module
542 AliVCaloCells * cells = fEvent->GetPHOSCells();
544 FillHistogram("hCenPHOSCells",fCentralityV0M,cells->GetNumberOfCells()) ;
545 FillHistogram("hCenTrack",fCentralityV0M,fEvent->GetNumberOfTracks()) ;
548 Int_t nCellModule[3] = {0,0,0};
549 for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
550 Int_t cellAbsId = cells->GetCellNumber(iCell);
551 Int_t relId[4] = {0,0,0,0};
552 fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
553 Int_t mod1 = relId[0];
554 Int_t cellX = relId[2];
555 Int_t cellZ = relId[3] ;
556 Float_t energy = cells->GetAmplitude(iCell);
557 FillHistogram("hCellEnergy",energy);
560 FillHistogram("hCellEnergyM1",cells->GetAmplitude(iCell));
561 FillHistogram("hCellNXZM1",cellX,cellZ,1.);
562 FillHistogram("hCellEXZM1",cellX,cellZ,energy);
566 FillHistogram("hCellEnergyM2",cells->GetAmplitude(iCell));
567 FillHistogram("hCellNXZM2",cellX,cellZ,1.);
568 FillHistogram("hCellEXZM2",cellX,cellZ,energy);
572 FillHistogram("hCellEnergyM3",cells->GetAmplitude(iCell));
573 FillHistogram("hCellNXZM3",cellX,cellZ,1.);
574 FillHistogram("hCellEXZM3",cellX,cellZ,energy);
577 FillHistogram("hCellMultEvent",nCellModule[0]+nCellModule[1]+nCellModule[2]);
578 FillHistogram("hCellMultEventM1",nCellModule[0]);
579 FillHistogram("hCellMultEventM2",nCellModule[1]);
580 FillHistogram("hCellMultEventM3",nCellModule[2]);
583 //_____________________________________________________________________________
584 void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
586 // clear (or create) array for holding events photons/clusters
588 fCaloPhotonsPHOS->Clear();
590 fCaloPhotonsPHOS = new TObjArray(200);
591 fCaloPhotonsPHOS->SetOwner();
595 AliVCaloCells* cells = dynamic_cast<AliVCaloCells*> (fEvent->GetPHOSCells());
596 for (Int_t i=0; i<fEvent->GetNumberOfCaloClusters(); i++) {
597 AliVCluster *clu = fEvent->GetCaloCluster(i);
599 if ( !clu->IsPHOS() || clu->E()< kMinClusterEnergy) continue; // reject cluster
603 // check if cell/channel is good.
605 clu->GetPosition(position);
606 TVector3 global(position) ;
608 fPHOSGeo->GlobalPos2RelId(global,relId) ;
609 Int_t mod = relId[0] ;
610 Int_t cellX = relId[2];
611 Int_t cellZ = relId[3] ;
612 if ( ! fModuleEnabled[mod-1] )
614 if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
615 continue ; // reject if not.
617 Double_t distBC=clu->GetDistanceToBadChannel();
618 if(distBC<kMinBCDistance)
621 FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
623 if(clu->GetNCells() < kMinNCells) continue ;
624 if(clu->GetM02() < kMinM02) continue ;
626 TLorentzVector lorentzMomentum;
628 ecore = CoreEnergy(clu,cells);
630 //if ESD, Apply re-Calibreation
631 Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0)
633 AliPHOSEsdCluster cluPHOS1( *(AliESDCaloCluster*) (clu) );
634 cluPHOS1.Recalibrate(fPHOSCalibData, static_cast<AliESDCaloCells*> (cells)); // modify the cell energies
635 Reclusterize(&cluPHOS1) ;
636 cluPHOS1.EvalAll(kLogWeight, fVertexVector); // recalculate the cluster parameters
637 cluPHOS1.SetE(fNonLinCorr->Eval(cluPHOS1.E()));// Users's nonlinearity
639 if(cluPHOS1.E()<0.3) continue; // check energy again
641 //correct misalignment
643 const Float_t shiftX[6]={0.,-2.3,-2.11,-1.53,0.,0.} ;
644 const Float_t shiftZ[6]={0.,-0.4, 0.52, 0.8,0.,0.} ;
645 fPHOSGeo->Global2Local(localPos,global,mod) ;
646 fPHOSGeo->Local2Global(mod,localPos.X()+shiftX[mod],localPos.Z()+shiftZ[mod],global);
647 position[0]=global.X() ;
648 position[1]=global.Y() ;
649 position[2]=global.Z() ;
650 cluPHOS1.SetPosition(position);
652 cluPHOS1.GetMomentum(lorentzMomentum ,origo);
654 //TODO: Check, this may be LHC10h specific:
655 if(mod==2) lorentzMomentum*=135.5/134.0 ;
656 if(mod==3) lorentzMomentum*=135.5/137.2 ;
657 if(mod==2) ecore*=135.5/134.0 ;
658 if(mod==3) ecore*=135.5/137.2 ;
661 else if (fEventAOD) { // is ! ESD, AOD.
662 AliESDCaloCluster* aodCluster = (AliESDCaloCluster*) (clu);
663 aodCluster->GetMomentum(lorentzMomentum ,origo);
666 AliError("(Calo)Cluster is neither ESD nor AOD");
671 FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.);
672 if(lorentzMomentum.E()>1.5){
673 FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.);
676 fCaloPhotonsPHOS->Add(new AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
677 AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
680 lorentzMomentum*= ecore/lorentzMomentum.E() ;
681 ph->SetMomV2(&lorentzMomentum) ;
682 ph->SetNCells(clu->GetNCells());
683 ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
684 //Evaluate CoreDispersion
685 Double_t m02=0.,m20=0. ;
686 EvalCoreLambdas(clu, cells, m02, m20) ;
687 ph->SetDisp2Bit(TestCoreLambda(clu->E(),m20,m02)) ; //Correct order m20,m02
688 // ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
690 FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
694 Double_t dx=clu->GetTrackDx() ;
695 Double_t dz=clu->GetTrackDz() ;
696 Bool_t cpvBit=kTRUE ; //No track matched by default. True means: not from charged, according to veto.
697 Bool_t cpvBit2=kTRUE ; //More Strict criterion
700 TArrayI * itracks = static_cast<AliESDCaloCluster*> (clu)->GetTracksMatched() ;
701 if(itracks->GetSize()>0){
702 Int_t iTr = itracks->At(0);
703 if(iTr>=0 && iTr<fEvent->GetNumberOfTracks()){
704 AliVParticle* track = fEvent->GetTrack(iTr);
705 Double_t pt = track->Pt() ;
706 Short_t charge = track->Charge() ;
707 Double_t r=TestCPV(dx, dz, pt, charge) ;
713 else if ( fEventAOD ) {
714 int nTracksMatched = clu->GetNTracksMatched();
715 if(nTracksMatched > 0) {
716 AliVTrack* track = dynamic_cast<AliVTrack*> (clu->GetTrackMatched(0));
718 Double_t pt = track->Pt();
719 Short_t charge = track->Charge();
720 Double_t r = TestCPV(dx, dz, pt, charge) ;
726 ph->SetCPVBit(cpvBit) ;
727 ph->SetCPV2Bit(cpvBit2) ;
729 FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.);
731 ph->SetEMCx(float(cellX)) ;
732 ph->SetEMCz(float(cellZ)) ;
733 // ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
734 ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolde
737 FillHistogram("hCenPHOS",fCentralityV0M, fCaloPhotonsPHOS->GetEntriesFast()) ;
739 //_____________________________________________________________________________
740 void AliAnalysisTaskPi0Flow::FillSelectedClusterHistograms()
742 for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
743 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
745 Double_t dphiA=ph1->Phi()-fRPV0A ;
746 while(dphiA<0)dphiA+=TMath::Pi() ;
747 while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
749 Double_t dphiC=ph1->Phi()-fRPV0C ;
750 while(dphiC<0)dphiC+=TMath::Pi() ;
751 while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
753 Double_t dphiT=ph1->Phi()-fRP ;
754 while(dphiT<0)dphiT+=TMath::Pi() ;
755 while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
757 Double_t pt = ph1->Pt() ;
758 Double_t ptcore = ph1->GetMomV2()->Pt() ;
760 FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),pt,dphiA) ;
761 FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),pt,dphiC) ;
763 FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),pt,dphiT) ;
764 FillHistogram(Form("hPhotPhiV0AAllcore_cen%d",fCentBin),ptcore,dphiA) ;
765 FillHistogram(Form("hPhotPhiV0CAllcore_cen%d",fCentBin),ptcore,dphiC) ;
767 FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCentBin),ptcore,dphiT) ;
769 FillHistogram(Form("hPhotAll_cen%d",fCentBin),pt) ;
770 FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ptcore) ;
771 if(ph1->IsntUnfolded()){
772 FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),pt) ;
773 FillHistogram(Form("hPhotPhiV0AAllwou_cen%d",fCentBin),pt,dphiA) ;
774 FillHistogram(Form("hPhotPhiV0CAllwou_cen%d",fCentBin),pt,dphiC) ;
776 FillHistogram(Form("hPhotPhiTPCAllwou_cen%d",fCentBin),pt,dphiT) ;
779 FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),pt,dphiA) ;
780 FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),pt,dphiC) ;
782 FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),pt,dphiT) ;
784 FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ptcore,dphiA) ;
785 FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ptcore,dphiC) ;
787 FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ptcore,dphiT) ;
789 FillHistogram(Form("hPhotCPV_cen%d",fCentBin),pt) ;
790 FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ptcore) ;
793 FillHistogram(Form("hPhotPhiV0ACPV2_cen%d",fCentBin),pt,dphiA) ;
794 FillHistogram(Form("hPhotPhiV0CCPV2_cen%d",fCentBin),pt,dphiC) ;
796 FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCentBin),pt,dphiT) ;
798 FillHistogram(Form("hPhotPhiV0ACPV2core_cen%d",fCentBin),ptcore,dphiA) ;
799 FillHistogram(Form("hPhotPhiV0CCPV2core_cen%d",fCentBin),ptcore,dphiC) ;
801 FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCentBin),ptcore,dphiT) ;
802 FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),pt) ;
803 FillHistogram(Form("hPhotCPV2core_cen%d",fCentBin),ptcore) ;
806 FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),pt,dphiA) ;
807 FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),pt,dphiC) ;
809 FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),pt,dphiT) ;
811 FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ptcore,dphiA) ;
812 FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ptcore,dphiC) ;
814 FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ptcore,dphiT) ;
816 if(ph1->IsntUnfolded()){
817 FillHistogram(Form("hPhotPhiV0ADispwou_cen%d",fCentBin),pt,dphiA) ;
818 FillHistogram(Form("hPhotPhiV0CDispwou_cen%d",fCentBin),pt,dphiC) ;
820 FillHistogram(Form("hPhotPhiTPCDispwou_cen%d",fCentBin),pt,dphiT) ;
823 FillHistogram(Form("hPhotDisp_cen%d",fCentBin),pt) ;
824 FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ptcore) ;
825 if(ph1->IsntUnfolded()){
826 FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),pt) ;
829 FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),pt,dphiA) ;
830 FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),pt,dphiC) ;
832 FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),pt,dphiT) ;
834 FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ptcore,dphiA) ;
835 FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ptcore,dphiC) ;
837 FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ptcore,dphiT) ;
839 FillHistogram(Form("hPhotBoth_cen%d",fCentBin),pt) ;
840 FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ptcore) ;
843 if(ph1->IsDisp2OK()){
844 FillHistogram(Form("hPhotPhiV0ADisp2_cen%d",fCentBin),pt,dphiA) ;
845 FillHistogram(Form("hPhotPhiV0CDisp2_cen%d",fCentBin),pt,dphiC) ;
847 FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCentBin),pt,dphiT) ;
848 FillHistogram(Form("hPhotPhiV0ADisp2core_cen%d",fCentBin),ptcore,dphiA) ;
849 FillHistogram(Form("hPhotPhiV0CDisp2core_cen%d",fCentBin),ptcore,dphiC) ;
851 FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCentBin),ptcore,dphiT) ;
853 FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),pt) ;
854 FillHistogram(Form("hPhotDisp2core_cen%d",fCentBin),ptcore) ;
856 FillHistogram(Form("hPhotPhiV0ABoth2_cen%d",fCentBin),pt,dphiA) ;
857 FillHistogram(Form("hPhotPhiV0CBoth2_cen%d",fCentBin),pt,dphiC) ;
859 FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCentBin),pt,dphiT) ;
861 FillHistogram(Form("hPhotPhiV0ABoth2core_cen%d",fCentBin),ptcore,dphiA) ;
862 FillHistogram(Form("hPhotPhiV0CBoth2core_cen%d",fCentBin),ptcore,dphiC) ;
864 FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCentBin),ptcore,dphiT) ;
866 FillHistogram(Form("hPhotBoth2_cen%d",fCentBin),pt) ;
867 FillHistogram(Form("hPhotBoth2core_cen%d",fCentBin),ptcore) ;
872 //_____________________________________________________________________________
873 void AliAnalysisTaskPi0Flow::ConsiderPi0s()
876 for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
877 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
878 for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
879 AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
880 TLorentzVector p12 = *ph1 + *ph2;
881 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
882 FillHistogram("hPHOSphi",fCentralityV0M,p12.Pt(),p12.Phi());
883 Double_t dphiA=p12.Phi()-fRPV0A ;
884 while(dphiA<0)dphiA+=TMath::Pi() ;
885 while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
887 Double_t dphiC=p12.Phi()-fRPV0C ;
888 while(dphiC<0)dphiC+=TMath::Pi() ;
889 while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
891 Double_t dphiT=p12.Phi()-fRP ;
892 while(dphiT<0)dphiT+=TMath::Pi() ;
893 while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
895 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
897 Double_t mcore=pv12.M() ;
898 Double_t pt=p12.Pt() ;
899 Double_t ptcore=pv12.Pt() ;
900 Double_t pt1=ph1->Pt() ;
901 Double_t pt2=ph2->Pt() ;
902 Double_t ptcore1=ph1->GetMomV2()->Pt() ;
903 Double_t ptcore2=ph2->GetMomV2()->Pt() ;
905 FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA) ;
906 FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC) ;
908 FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
910 FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
911 FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
913 FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
916 FillHistogram(Form("hPi0All_cen%d",fCentBin),m,pt) ;
917 FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
918 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
919 FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),m,pt) ;
920 FillHistogram(Form("hMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA) ;
921 FillHistogram(Form("hMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC) ;
923 FillHistogram(Form("hMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT) ;
926 FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt1) ;
927 FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt2) ;
928 FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore1) ;
929 FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore2) ;
930 if(ph1->IsntUnfolded())
931 FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt1) ;
932 if(ph2->IsntUnfolded())
933 FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt2) ;
935 FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt1) ;
936 FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
939 FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt2) ;
940 FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
943 FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt1) ;
944 FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
947 FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt2) ;
948 FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
951 FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt1) ;
952 if(ph1->IsntUnfolded()){
953 FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt1) ;
955 FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
958 FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt2) ;
959 if(ph1->IsntUnfolded()){
960 FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt2) ;
962 FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
964 if(ph1->IsDisp2OK()){
965 FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt1) ;
966 FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
968 if(ph2->IsDisp2OK()){
969 FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt2) ;
970 FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
972 if(ph1->IsDispOK() && ph1->IsCPVOK()){
973 FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt1) ;
974 FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore1) ;
976 if(ph2->IsDispOK() && ph2->IsCPVOK()){
977 FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt2) ;
978 FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore2) ;
980 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
981 FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt1) ;
982 FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore1) ;
984 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
985 FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt2) ;
986 FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore2) ;
991 FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),m,pt) ;
994 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
995 snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
996 FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
997 FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
999 FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
1001 FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1002 FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1004 FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1006 FillHistogram(Form("hPi0CPV_cen%d",fCentBin),m,pt) ;
1007 FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
1010 FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),m,pt) ;
1013 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1014 FillHistogram(Form("hMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA) ;
1015 FillHistogram(Form("hMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC) ;
1017 FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT) ;
1018 FillHistogram(Form("hMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1019 FillHistogram(Form("hMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1021 FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1023 FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),m,pt) ;
1024 FillHistogram(Form("hPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
1026 FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
1029 if(ph1->IsDispOK() && ph2->IsDispOK()){
1030 snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
1031 FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1032 FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
1034 FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
1036 FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1037 FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1039 FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1041 FillHistogram(Form("hPi0Disp_cen%d",fCentBin),m,pt) ;
1042 FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),mcore, ptcore) ;
1044 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1045 FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),m,pt) ;
1047 FillHistogram(Form("hMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA) ;
1048 FillHistogram(Form("hMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC) ;
1050 FillHistogram(Form("hMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT) ;
1054 FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),m,pt) ;
1056 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1057 FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1058 FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
1060 FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
1062 FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1063 FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1065 FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1067 FillHistogram(Form("hPi0Both_cen%d",fCentBin),m,pt) ;
1068 FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),mcore,ptcore) ;
1071 snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
1072 FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),m,pt) ;
1074 if(ph1->Module()==1 && ph2->Module()==1)
1075 FillHistogram("hPi0M11",m,pt );
1076 else if(ph1->Module()==2 && ph2->Module()==2)
1077 FillHistogram("hPi0M22",m,pt );
1078 else if(ph1->Module()==3 && ph2->Module()==3)
1079 FillHistogram("hPi0M33",m,pt );
1080 else if(ph1->Module()==1 && ph2->Module()==2)
1081 FillHistogram("hPi0M12",m,pt );
1082 else if(ph1->Module()==1 && ph2->Module()==3)
1083 FillHistogram("hPi0M13",m,pt );
1084 else if(ph1->Module()==2 && ph2->Module()==3)
1085 FillHistogram("hPi0M23",m,pt );
1092 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1093 FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),m,pt) ;
1094 FillHistogram(Form("hPi0Disp2core_cen%d",fCentBin),mcore, ptcore) ;
1096 FillHistogram(Form("hMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA) ;
1097 FillHistogram(Form("hMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC) ;
1099 FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT) ;
1101 FillHistogram(Form("hMassPtV0ADisp2core_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1102 FillHistogram(Form("hMassPtV0CDisp2core_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1104 FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1106 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1107 FillHistogram(Form("hMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA) ;
1108 FillHistogram(Form("hMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC) ;
1110 FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT) ;
1112 FillHistogram(Form("hMassPtV0ABoth2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1113 FillHistogram(Form("hMassPtV0CBoth2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1115 FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1117 FillHistogram(Form("hPi0Both2_cen%d",fCentBin),m,pt) ;
1118 FillHistogram(Form("hPi0Both2core_cen%d",fCentBin),mcore,ptcore) ;
1125 //_____________________________________________________________________________
1126 void AliAnalysisTaskPi0Flow::ConsiderPi0sMix()
1130 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1132 for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
1133 AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1134 for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1135 TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
1136 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1137 AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1138 TLorentzVector p12 = *ph1 + *ph2;
1139 TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1141 Double_t dphiA=p12.Phi()-fRPV0A ;
1142 while(dphiA<0)dphiA+=TMath::Pi() ;
1143 while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
1145 Double_t dphiC=p12.Phi()-fRPV0C ;
1146 while(dphiC<0)dphiC+=TMath::Pi() ;
1147 while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
1149 Double_t dphiT=p12.Phi()-fRP ;
1150 while(dphiT<0)dphiT+=TMath::Pi() ;
1151 while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
1154 Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1155 Double_t m=p12.M() ;
1156 Double_t mcore=pv12.M() ;
1157 Double_t pt=p12.Pt() ;
1158 Double_t ptcore=pv12.Pt() ;
1159 Double_t pt1=ph1->Pt() ;
1160 Double_t pt2=ph2->Pt() ;
1161 Double_t ptcore1=ph1->GetMomV2()->Pt() ;
1162 Double_t ptcore2=ph2->GetMomV2()->Pt() ;
1165 snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
1166 FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA) ;
1167 FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC) ;
1169 FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
1171 FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),mcore, ptcore, dphiA) ;
1172 FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),mcore, ptcore, dphiC) ;
1174 FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),mcore, ptcore, dphiT) ;
1176 FillHistogram(Form("hMiPi0All_cen%d",fCentBin),m,pt) ;
1177 FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
1178 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1179 FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),m,pt) ;
1180 FillHistogram(Form("hMiMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA) ;
1181 FillHistogram(Form("hMiMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC) ;
1183 FillHistogram(Form("hMiMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT) ;
1186 FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt1) ;
1187 FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt2) ;
1188 FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore1) ;
1189 FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore2) ;
1190 if(ph1->IsntUnfolded())
1191 FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt1) ;
1192 if(ph2->IsntUnfolded())
1193 FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt2) ;
1195 FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt1) ;
1196 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
1199 FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt2) ;
1200 FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
1202 if(ph1->IsCPV2OK()){
1203 FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt1) ;
1204 FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore1) ;
1206 if(ph2->IsCPV2OK()){
1207 FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt2) ;
1208 FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
1210 if(ph1->IsDispOK()){
1211 FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt1) ;
1212 if(ph1->IsntUnfolded()){
1213 FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt1) ;
1215 FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
1217 if(ph2->IsDispOK()){
1218 FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt2) ;
1219 if(ph1->IsntUnfolded()){
1220 FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt2) ;
1222 FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
1224 if(ph1->IsDisp2OK()){
1225 FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt1) ;
1226 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
1228 if(ph2->IsDisp2OK()){
1229 FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt2) ;
1230 FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore2) ;
1232 if(ph1->IsDispOK() && ph1->IsCPVOK()){
1233 snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1234 FillHistogram(key,m,pt1) ;
1235 snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1236 FillHistogram(key,mcore,ptcore1) ;
1238 if(ph2->IsDispOK() && ph2->IsCPVOK()){
1239 snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1240 FillHistogram(key,m,pt2) ;
1241 snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1242 FillHistogram(key,mcore,ptcore2) ;
1244 if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1245 FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt1) ;
1246 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore1) ;
1248 if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1249 FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt2) ;
1250 FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore2) ;
1256 FillHistogram(Form("hMiPi0All_a07_cen%d",fCentBin),m,pt) ;
1258 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1259 FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
1260 FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
1262 FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
1264 FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1265 FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1267 FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1269 FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),m,pt) ;
1270 FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
1273 FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),m,pt) ;
1276 if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1277 FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),m,pt) ;
1278 FillHistogram(Form("hMiPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
1280 FillHistogram(Form("hMiMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA) ;
1281 FillHistogram(Form("hMiMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC) ;
1283 FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT) ;
1284 FillHistogram(Form("hMiMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1285 FillHistogram(Form("hMiMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1287 FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1290 FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
1293 if(ph1->IsDispOK() && ph2->IsDispOK()){
1294 FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1295 FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
1297 FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
1299 FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1300 FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1302 FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1305 FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),m,pt) ;
1306 FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1307 if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1308 FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),m,pt) ;
1309 FillHistogram(Form("hMiMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA) ;
1310 FillHistogram(Form("hMiMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC) ;
1312 FillHistogram(Form("hMiMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT) ;
1316 FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),m,pt) ;
1318 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1319 FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1320 FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
1322 FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
1324 FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1325 FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1327 FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1329 FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),m,pt) ;
1330 FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1333 FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),m,pt) ;
1338 if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1339 FillHistogram(Form("hMiMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA) ;
1340 FillHistogram(Form("hMiMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC) ;
1342 FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT) ;
1344 FillHistogram(Form("hMiMassPtV0ADisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1345 FillHistogram(Form("hMiMassPtV0CDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1347 FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1350 FillHistogram(Form("hMiPi0Disp2_cen%d",fCentBin),m,pt) ;
1351 FillHistogram(Form("hMiPi0Disp2core_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1353 if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1354 FillHistogram(Form("hMiMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA) ;
1355 FillHistogram(Form("hMiMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC) ;
1357 FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT) ;
1359 FillHistogram(Form("hMiMassPtV0ABoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1360 FillHistogram(Form("hMiMassPtV0CBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1362 FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1364 FillHistogram(Form("hMiPi0Both2_cen%d",fCentBin),m,pt) ;
1365 FillHistogram(Form("hMiPi0Both2core_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1373 //_____________________________________________________________________________
1374 void AliAnalysisTaskPi0Flow::UpdateLists()
1376 //Now we either add current events to stack or remove
1377 //If no photons in current event - no need to add it to mixed
1379 TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1381 AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1382 if(fCaloPhotonsPHOS->GetEntriesFast()>0){
1383 arrayList->AddFirst(fCaloPhotonsPHOS) ;
1385 if(arrayList->GetEntries() > fCentNMixed[fCentBin]){ // Remove redundant events
1386 TObjArray * tmp = static_cast<TObjArray*>(arrayList->Last()) ;
1387 arrayList->RemoveLast() ;
1388 delete tmp ; // TODO: may conflict with delete done by list being owner.
1392 fCaloPhotonsPHOS->Clear(); // TODO: redundant???
1394 //_____________________________________________________________________________
1395 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x)const{
1397 TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
1401 AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
1403 //_____________________________________________________________________________
1404 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y)const{
1406 TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
1410 AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
1413 //_____________________________________________________________________________
1414 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1415 //Fills 1D histograms with key
1416 TObject * obj = fOutputContainer->FindObject(key);
1418 TH2 * th2 = dynamic_cast<TH2*> (obj);
1420 th2->Fill(x, y, z) ;
1424 TH3 * th3 = dynamic_cast<TH3*> (obj);
1426 th3->Fill(x, y, z) ;
1430 AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
1433 //_____________________________________________________________________________
1434 AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
1436 fEvent = InputEvent();
1438 AliError("Event could not be retrieved");
1439 PostData(1, fOutputContainer);
1445 //___________________________________________________________________________
1446 AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
1449 AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
1451 AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
1453 fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
1458 //___________________________________________________________________________
1459 Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
1461 /* fCentBin=1+Int_t(centralityV0M/100. *kNCenBins) ;
1462 if(centralityV0M < 5. || fCentBin < 0)
1464 if(fCentBin > kNCenBins-1)
1465 fCentBin = kNCenBins-1 ;
1467 int lastBinUpperIndex = fCentEdges.GetSize() -1;
1468 if( centralityV0M > fCentEdges[lastBinUpperIndex] ) {
1470 AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1471 return lastBinUpperIndex-1;
1473 if( centralityV0M < fCentEdges[0] ) {
1475 AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1479 fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1483 //___________________________________________________________________________
1484 Int_t AliAnalysisTaskPi0Flow::GetRPBin()
1488 averageRP = fRP ;// If possible, it is better to have EP bin from TPC
1489 // to have similar events for miximng (including jets etc) (fRPV0A+fRPV0C+fRP) /3.;
1491 averageRP = (fRPV0A+fRPV0C) /2.;
1493 fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1495 if(fEMRPBin> (Int_t) fNEMRPBins-1)
1496 fEMRPBin=fNEMRPBins-1 ;
1501 AliInfo(Form("Event Mixing Reaction Plane bin is: %d", fEMRPBin));
1507 //_____________________________________________________________________________
1508 void AliAnalysisTaskPi0Flow::LogProgress(int step)
1511 AliInfo(Form("step %d completed", step));
1513 // the +0.5 is not realy neccisarry, but oh well... -henrik
1514 //FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1515 //FillHistogram("hTotSelEvents", step+0.5);
1518 void AliAnalysisTaskPi0Flow::LogSelection(int step, int internalRunNumber)
1521 // AliInfo(Form("step %d completed", step));
1523 // the +0.5 is not realy neccisarry, but oh well... -henrik
1524 FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1525 FillHistogram("hTotSelEvents", step+0.5);
1529 //___________________________________________________________________________
1530 Int_t AliAnalysisTaskPi0Flow::ConvertToInternalRunNumber(Int_t run){
1531 if( kLHC11h == fPeriod ) {
1533 case 170593 : return 179 ;
1534 case 170572 : return 178 ;
1535 case 170556 : return 177 ;
1536 case 170552 : return 176 ;
1537 case 170546 : return 175 ;
1538 case 170390 : return 174 ;
1539 case 170389 : return 173 ;
1540 case 170388 : return 172 ;
1541 case 170387 : return 171 ;
1542 case 170315 : return 170 ;
1543 case 170313 : return 169 ;
1544 case 170312 : return 168 ;
1545 case 170311 : return 167 ;
1546 case 170309 : return 166 ;
1547 case 170308 : return 165 ;
1548 case 170306 : return 164 ;
1549 case 170270 : return 163 ;
1550 case 170269 : return 162 ;
1551 case 170268 : return 161 ;
1552 case 170267 : return 160 ;
1553 case 170264 : return 159 ;
1554 case 170230 : return 158 ;
1555 case 170228 : return 157 ;
1556 case 170208 : return 156 ;
1557 case 170207 : return 155 ;
1558 case 170205 : return 154 ;
1559 case 170204 : return 153 ;
1560 case 170203 : return 152 ;
1561 case 170195 : return 151 ;
1562 case 170193 : return 150 ;
1563 case 170163 : return 149 ;
1564 case 170162 : return 148 ;
1565 case 170159 : return 147 ;
1566 case 170155 : return 146 ;
1567 case 170152 : return 145 ;
1568 case 170091 : return 144 ;
1569 case 170089 : return 143 ;
1570 case 170088 : return 142 ;
1571 case 170085 : return 141 ;
1572 case 170084 : return 140 ;
1573 case 170083 : return 139 ;
1574 case 170081 : return 138 ;
1575 case 170040 : return 137 ;
1576 case 170038 : return 136 ;
1577 case 170036 : return 135 ;
1578 case 170027 : return 134 ;
1579 case 169981 : return 133 ;
1580 case 169975 : return 132 ;
1581 case 169969 : return 131 ;
1582 case 169965 : return 130 ;
1583 case 169961 : return 129 ;
1584 case 169956 : return 128 ;
1585 case 169926 : return 127 ;
1586 case 169924 : return 126 ;
1587 case 169923 : return 125 ;
1588 case 169922 : return 124 ;
1589 case 169919 : return 123 ;
1590 case 169918 : return 122 ;
1591 case 169914 : return 121 ;
1592 case 169859 : return 120 ;
1593 case 169858 : return 119 ;
1594 case 169855 : return 118 ;
1595 case 169846 : return 117 ;
1596 case 169838 : return 116 ;
1597 case 169837 : return 115 ;
1598 case 169835 : return 114 ;
1599 case 169683 : return 113 ;
1600 case 169628 : return 112 ;
1601 case 169591 : return 111 ;
1602 case 169590 : return 110 ;
1603 case 169588 : return 109 ;
1604 case 169587 : return 108 ;
1605 case 169586 : return 107 ;
1606 case 169584 : return 106 ;
1607 case 169557 : return 105 ;
1608 case 169555 : return 104 ;
1609 case 169554 : return 103 ;
1610 case 169553 : return 102 ;
1611 case 169550 : return 101 ;
1612 case 169515 : return 100 ;
1613 case 169512 : return 99 ;
1614 case 169506 : return 98 ;
1615 case 169504 : return 97 ;
1616 case 169498 : return 96 ;
1617 case 169475 : return 95 ;
1618 case 169420 : return 94 ;
1619 case 169419 : return 93 ;
1620 case 169418 : return 92 ;
1621 case 169417 : return 91 ;
1622 case 169415 : return 90 ;
1623 case 169411 : return 89 ;
1624 case 169238 : return 88 ;
1625 case 169236 : return 87 ;
1626 case 169167 : return 86 ;
1627 case 169160 : return 85 ;
1628 case 169156 : return 84 ;
1629 case 169148 : return 83 ;
1630 case 169145 : return 82 ;
1631 case 169144 : return 81 ;
1632 case 169143 : return 80 ;
1633 case 169138 : return 79 ;
1634 case 169099 : return 78 ;
1635 case 169094 : return 77 ;
1636 case 169091 : return 76 ;
1637 case 169045 : return 75 ;
1638 case 169044 : return 74 ;
1639 case 169040 : return 73 ;
1640 case 169035 : return 72 ;
1641 case 168992 : return 71 ;
1642 case 168988 : return 70 ;
1643 case 168984 : return 69 ;
1644 case 168826 : return 68 ;
1645 case 168777 : return 67 ;
1646 case 168514 : return 66 ;
1647 case 168512 : return 65 ;
1648 case 168511 : return 64 ;
1649 case 168467 : return 63 ;
1650 case 168464 : return 62 ;
1651 case 168461 : return 61 ;
1652 case 168460 : return 60 ;
1653 case 168458 : return 59 ;
1654 case 168362 : return 58 ;
1655 case 168361 : return 57 ;
1656 case 168356 : return 56 ;
1657 case 168342 : return 55 ;
1658 case 168341 : return 54 ;
1659 case 168325 : return 53 ;
1660 case 168322 : return 52 ;
1661 case 168318 : return 51 ;
1662 case 168311 : return 50 ;
1663 case 168310 : return 49 ;
1664 case 168213 : return 48 ;
1665 case 168212 : return 47 ;
1666 case 168208 : return 46 ;
1667 case 168207 : return 45 ;
1668 case 168206 : return 44 ;
1669 case 168205 : return 43 ;
1670 case 168204 : return 42 ;
1671 case 168203 : return 41 ;
1672 case 168181 : return 40 ;
1673 case 168177 : return 39 ;
1674 case 168175 : return 38 ;
1675 case 168173 : return 37 ;
1676 case 168172 : return 36 ;
1677 case 168171 : return 35 ;
1678 case 168115 : return 34 ;
1679 case 168108 : return 33 ;
1680 case 168107 : return 32 ;
1681 case 168105 : return 31 ;
1682 case 168104 : return 30 ;
1683 case 168103 : return 29 ;
1684 case 168076 : return 28 ;
1685 case 168069 : return 27 ;
1686 case 168068 : return 26 ;
1687 case 168066 : return 25 ;
1688 case 167988 : return 24 ;
1689 case 167987 : return 23 ;
1690 case 167986 : return 22 ;
1691 case 167985 : return 21 ;
1692 case 167921 : return 20 ;
1693 case 167920 : return 19 ;
1694 case 167915 : return 18 ;
1695 case 167909 : return 17 ;
1696 case 167903 : return 16 ;
1697 case 167902 : return 15 ;
1698 case 167818 : return 14 ;
1699 case 167814 : return 13 ;
1700 case 167813 : return 12 ;
1701 case 167808 : return 11 ;
1702 case 167807 : return 10 ;
1703 case 167806 : return 9 ;
1704 case 167713 : return 8 ;
1705 case 167712 : return 7 ;
1706 case 167711 : return 6 ;
1707 case 167706 : return 5 ;
1708 case 167693 : return 4 ;
1709 case 166532 : return 3 ;
1710 case 166530 : return 2 ;
1711 case 166529 : return 1 ;
1713 default : return 199;
1716 if( kLHC10h == fPeriod ) {
1718 case 139517 : return 137;
1719 case 139514 : return 136;
1720 case 139513 : return 135;
1721 case 139511 : return 134;
1722 case 139510 : return 133;
1723 case 139507 : return 132;
1724 case 139505 : return 131;
1725 case 139504 : return 130;
1726 case 139503 : return 129;
1727 case 139470 : return 128;
1728 case 139467 : return 127;
1729 case 139466 : return 126;
1730 case 139465 : return 125;
1731 case 139440 : return 124;
1732 case 139439 : return 123;
1733 case 139438 : return 122;
1734 case 139437 : return 121;
1735 case 139360 : return 120;
1736 case 139329 : return 119;
1737 case 139328 : return 118;
1738 case 139314 : return 117;
1739 case 139311 : return 116;
1740 case 139310 : return 115;
1741 case 139309 : return 114;
1742 case 139308 : return 113;
1743 case 139173 : return 112;
1744 case 139172 : return 111;
1745 case 139110 : return 110;
1746 case 139107 : return 109;
1747 case 139105 : return 108;
1748 case 139104 : return 107;
1749 case 139042 : return 106;
1750 case 139038 : return 105;
1751 case 139037 : return 104;
1752 case 139036 : return 103;
1753 case 139029 : return 102;
1754 case 139028 : return 101;
1755 case 138983 : return 100;
1756 case 138982 : return 99;
1757 case 138980 : return 98;
1758 case 138979 : return 97;
1759 case 138978 : return 96;
1760 case 138977 : return 95;
1761 case 138976 : return 94;
1762 case 138973 : return 93;
1763 case 138972 : return 92;
1764 case 138965 : return 91;
1765 case 138924 : return 90;
1766 case 138872 : return 89;
1767 case 138871 : return 88;
1768 case 138870 : return 87;
1769 case 138837 : return 86;
1770 case 138830 : return 85;
1771 case 138828 : return 84;
1772 case 138826 : return 83;
1773 case 138796 : return 82;
1774 case 138795 : return 81;
1775 case 138742 : return 80;
1776 case 138732 : return 79;
1777 case 138730 : return 78;
1778 case 138666 : return 77;
1779 case 138662 : return 76;
1780 case 138653 : return 75;
1781 case 138652 : return 74;
1782 case 138638 : return 73;
1783 case 138624 : return 72;
1784 case 138621 : return 71;
1785 case 138583 : return 70;
1786 case 138582 : return 69;
1787 case 138579 : return 68;
1788 case 138578 : return 67;
1789 case 138534 : return 66;
1790 case 138469 : return 65;
1791 case 138442 : return 64;
1792 case 138439 : return 63;
1793 case 138438 : return 62;
1794 case 138396 : return 61;
1795 case 138364 : return 60;
1796 case 138359 : return 59;
1797 case 138275 : return 58;
1798 case 138225 : return 57;
1799 case 138201 : return 56;
1800 case 138200 : return 55;
1801 case 138197 : return 54;
1802 case 138192 : return 53;
1803 case 138190 : return 52;
1804 case 138154 : return 51;
1805 case 138153 : return 50;
1806 case 138151 : return 49;
1807 case 138150 : return 48;
1808 case 138126 : return 47;
1809 case 138125 : return 46;
1810 case 137848 : return 45;
1811 case 137847 : return 44;
1812 case 137844 : return 43;
1813 case 137843 : return 42;
1814 case 137752 : return 41;
1815 case 137751 : return 40;
1816 case 137748 : return 39;
1817 case 137724 : return 38;
1818 case 137722 : return 37;
1819 case 137718 : return 36;
1820 case 137704 : return 35;
1821 case 137693 : return 34;
1822 case 137692 : return 33;
1823 case 137691 : return 32;
1824 case 137689 : return 31;
1825 case 137686 : return 30;
1826 case 137685 : return 29;
1827 case 137639 : return 28;
1828 case 137638 : return 27;
1829 case 137608 : return 26;
1830 case 137595 : return 25;
1831 case 137549 : return 24;
1832 case 137546 : return 23;
1833 case 137544 : return 22;
1834 case 137541 : return 21;
1835 case 137539 : return 20;
1836 case 137531 : return 19;
1837 case 137530 : return 18;
1838 case 137443 : return 17;
1839 case 137441 : return 16;
1840 case 137440 : return 15;
1841 case 137439 : return 14;
1842 case 137434 : return 13;
1843 case 137432 : return 12;
1844 case 137431 : return 11;
1845 case 137430 : return 10;
1846 case 137366 : return 9;
1847 case 137243 : return 8;
1848 case 137236 : return 7;
1849 case 137235 : return 6;
1850 case 137232 : return 5;
1851 case 137231 : return 4;
1852 case 137165 : return 3;
1853 case 137162 : return 2;
1854 case 137161 : return 1;
1855 default : return 199;
1858 if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) ) {
1859 AliWarning("Period not defined");
1863 //_____________________________________________________________________________
1864 Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1866 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1867 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1868 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1869 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1870 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1871 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1872 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1873 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1874 return (R2<2.5*2.5) ;
1877 //_____________________________________________________________________________
1878 Bool_t AliAnalysisTaskPi0Flow::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1880 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1881 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1882 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1883 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1884 Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1885 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1886 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1887 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1888 return (R2<1.5*1.5) ;
1891 //____________________________________________________________________________
1892 TList* AliAnalysisTaskPi0Flow::GetCaloPhotonsPHOSList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin)
1894 int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins
1895 + centBin * fNEMRPBins
1897 if( fCaloPhotonsPHOSLists->At(offset) ) { // list exists
1898 TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1900 AliError("object in fCaloPhotonsPHOSLists at %i did not cast");
1903 else {// no list for this bin has been created, yet
1904 TList* list = new TList();
1906 fCaloPhotonsPHOSLists->AddAt(list, offset);
1911 //____________________________________________________________________________
1912 Double_t AliAnalysisTaskPi0Flow::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1913 //Parameterization of LHC10h period
1918 Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1919 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1920 Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1921 Double_t mf = 0.; //
1922 if(fEventAOD) mf = fEventAOD->GetMagneticField(); //Positive for ++ and negative for --
1923 else if(fEventESD) mf = fEventESD->GetMagneticField(); //Positive for ++ and negative for --
1926 if(mf<0.){ //field --
1929 meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1931 meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1936 meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1938 meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1941 Double_t rz=(dz-meanZ)/sz ;
1942 Double_t rx=(dx-meanX)/sx ;
1943 return TMath::Sqrt(rx*rx+rz*rz) ;
1945 //____________________________________________________________________________
1946 void AliAnalysisTaskPi0Flow::SetFlatteningData(){
1947 //Read objects with flattening parameters
1948 AliOADBContainer flatContainer("phosFlat");
1949 flatContainer.InitFromFile(fEPcalibFileName.Data(),"phosFlat");
1950 TObjArray *maps = (TObjArray*)flatContainer.GetObject(fRunNumber,"phosFlat");
1952 AliError(Form("Can not read Flattening for run %d. \n From file >%s<\n",fRunNumber,fEPcalibFileName.Data())) ;
1955 AliInfo(Form("Setting PHOS flattening with name %s \n",maps->GetName())) ;
1956 AliEPFlattener * h = (AliEPFlattener*)maps->At(0) ;
1957 if(fTPCFlat) delete fTPCFlat ;
1958 fTPCFlat = new AliEPFlattener() ;
1960 h = (AliEPFlattener*)maps->At(1) ;
1961 if(fV0AFlat) delete fV0AFlat ;
1962 fV0AFlat = new AliEPFlattener() ;
1964 h = (AliEPFlattener*)maps->At(2) ;
1965 if(fV0CFlat) delete fV0CFlat ;
1966 fV0CFlat = new AliEPFlattener() ;
1971 //____________________________________________________________________________
1972 Double_t AliAnalysisTaskPi0Flow::ApplyFlattening(Double_t phi, Double_t c){
1975 return fTPCFlat->MakeFlat(phi,c);
1979 //____________________________________________________________________________
1980 Double_t AliAnalysisTaskPi0Flow::ApplyFlatteningV0A(Double_t phi, Double_t c){
1983 return fV0AFlat->MakeFlat(phi,c);
1987 //____________________________________________________________________________
1988 Double_t AliAnalysisTaskPi0Flow::ApplyFlatteningV0C(Double_t phi, Double_t c){
1991 return fV0CFlat->MakeFlat(phi,c);
1995 //____________________________________________________________________________
1996 Double_t AliAnalysisTaskPi0Flow::CoreEnergy(AliVCluster * clu, AliVCaloCells * cells)
1998 //calculate energy of the cluster in the circle with radius distanceCut around the maximum
2000 //Can not use already calculated coordinates?
2001 //They have incidence correction...
2002 const Double_t distanceCut =3.5 ;
2003 const Double_t logWeight=4.5 ;
2005 const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2006 // Calculates the center of gravity in the local PHOS-module coordinates
2008 const Int_t mulDigit=clu->GetNCells() ;
2009 Double_t xc[mulDigit] ;
2010 Double_t zc[mulDigit] ;
2011 Double_t ei[mulDigit] ;
2014 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2018 fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
2019 fPHOSGeo->RelPosInModule(relid, xi, zi);
2022 ei[iDigit]=elist[iDigit]*cells->GetCellAmplitude(clu->GetCellsAbsId()[iDigit]);
2024 printf("%f ",ei[iDigit]);
2025 if (clu->E()>0 && ei[iDigit]>0) {
2026 Float_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
2027 x += xc[iDigit] * w ;
2028 z += zc[iDigit] * w ;
2037 for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
2038 Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
2039 if(distance < distanceCut)
2040 coreE += ei[iDigit] ;
2042 //Apply non-linearity correction
2043 return fNonLinCorr->Eval(coreE) ;
2045 //____________________________________________________________________________
2046 Bool_t AliAnalysisTaskPi0Flow::AreNeibors(Int_t id1,Int_t id2){
2047 // return true if absId are "Neighbors" (adjacent, including diagornaly,)
2051 fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
2054 fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
2056 // if inside the same PHOS module
2057 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) {
2058 const Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
2059 const Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
2061 // and if diff in both direction is 1 or less
2062 if (( coldiff <= 1 ) && ( rowdiff <= 1 ))
2063 return true; // are neighbors
2069 //____________________________________________________________________________
2070 void AliAnalysisTaskPi0Flow::Reclusterize(AliVCluster * clu){
2071 //Re-clusterize to make continues cluster
2073 const Int_t oldMulDigit=clu->GetNCells() ;
2074 Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2075 UShort_t * dlist = clu->GetCellsAbsId();
2077 Int_t index[oldMulDigit] ;
2078 Bool_t used[oldMulDigit] ;
2079 for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
2083 for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
2084 if(eMax<elist[iDigit]){
2090 if(inClu==0){ //empty cluster
2093 used[index[0]]=kTRUE ; //mark as used
2094 for(Int_t i=0; i<inClu; i++){
2095 for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
2096 if(used[iDigit]) //already used
2098 if(AreNeibors(dlist[index[i]],dlist[iDigit])){
2099 index[inClu]= iDigit ;
2101 used[iDigit]=kTRUE ;
2106 if(inClu==oldMulDigit) //no need to modify
2109 clu->SetNCells(inClu);
2111 UShort_t tmpD[oldMulDigit] ;
2112 Double_t tmpE[oldMulDigit] ;
2113 for(Int_t i=0; i<oldMulDigit; i++){
2117 //change order of digits in list so that
2118 //first inClu cells were true ones
2119 for(Int_t i=0; i<inClu; i++){
2120 dlist[i]=tmpD[index[i]] ;
2121 elist[i]=tmpE[index[i]] ;
2127 //_____________________________________________________________________________
2128 void AliAnalysisTaskPi0Flow::SetMisalignment(){
2129 // sets the misalignment vertex if ESD
2131 for(Int_t mod=0; mod<5; mod++) {
2132 const TGeoHMatrix* modMatrix = fEvent->GetPHOSMatrix(mod);
2135 AliInfo(Form("no PHOS Geometric Misalignment Matrix for module %d", mod));
2139 fPHOSGeo->SetMisalMatrix(modMatrix, mod);
2141 AliInfo(Form("PHOS Geometric Misalignment Matrix set for module %d", mod));
2147 //_____________________________________________________________________________
2148 void AliAnalysisTaskPi0Flow::SetV0Calibration(){
2149 // assigns: fMultV0, fV0Cpol, fV0Apol, fMeanQ, and fWidthQ
2151 if ( ! fManualV0EPCalc ) {
2153 AliInfo("Not setting V0Calibration, only needed for manual V0 EP Calculation");
2157 int runNumber = this->fRunNumber;
2159 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
2160 TFile *foadb = TFile::Open(oadbfilename.Data());
2163 AliError(Form("OADB file %s cannot be opened\n", oadbfilename.Data()));
2164 AliError("V0 Calibration not set !\n");
2168 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
2170 AliError("OADB object hMultV0BefCorr is not available in the file");
2171 AliError("V0 Calibration not set!\n");
2175 if(!(cont->GetObject(runNumber))){
2176 AliError(Form("OADB object hMultV0BefCorr is not available for run %i, trying 137366)",runNumber));
2179 if(!(cont->GetObject(runNumber))){
2180 AliError(Form("OADB object hMultV0BefCorr is not available for run %i ",runNumber));
2181 AliError("V0 Calibration not set!\n");
2185 if( fDebug ) AliInfo("Setting V0 calibration") ;
2186 fMultV0 = ((TH2F *) cont->GetObject(runNumber))->ProfileX();
2188 TF1 *fpol0 = new TF1("fpol0","pol0");
2189 fMultV0->Fit(fpol0,"Q0","",0,31);
2190 fV0Cpol = fpol0->GetParameter(0);
2191 fMultV0->Fit(fpol0,"Q0","",32,64);
2192 fV0Apol = fpol0->GetParameter(0);
2194 for(Int_t iside=0;iside<2;iside++){
2195 for(Int_t icoord=0;icoord<2;icoord++){
2196 for(Int_t i=0;i < kNCenBins;i++){
2198 if(iside==0 && icoord==0)
2199 snprintf(namecont,100,"hQxc2_%i",i);
2200 else if(iside==1 && icoord==0)
2201 snprintf(namecont,100,"hQxa2_%i",i);
2202 else if(iside==0 && icoord==1)
2203 snprintf(namecont,100,"hQyc2_%i",i);
2204 else if(iside==1 && icoord==1)
2205 snprintf(namecont,100,"hQya2_%i",i);
2207 cont = (AliOADBContainer*) foadb->Get(namecont);
2209 AliError(Form("OADB object %s is not available in the file %s", namecont, oadbfilename.Data()));
2210 AliError("V0 Calibration not fully set!\n");
2214 if(!(cont->GetObject(runNumber))){
2215 AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2218 if(!(cont->GetObject(runNumber))){
2219 AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2220 AliError("V0 Calibration not fully set!\n");
2223 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetMean();
2224 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetRMS();
2227 // if(iside==0 && icoord==0)
2228 // snprintf(namecont,100,"hQxc3_%i",i);
2229 // else if(iside==1 && icoord==0)
2230 // snprintf(namecont,100,"hQxa3_%i",i);
2231 // else if(iside==0 && icoord==1)
2232 // snprintf(namecont,100,"hQyc3_%i",i);
2233 // else if(iside==1 && icoord==1)
2234 // snprintf(namecont,100,"hQya3_%i",i);
2236 // cont = (AliOADBContainer*) foadb->Get(namecont);
2238 // AliError(Form("OADB object %s is not available in the file",namecont));
2239 // AliError("V0 Calibration not fully set!\n");
2243 // if(!(cont->GetObject(runNumber))){
2244 // AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2245 // runNumber = 137366;
2247 // if(!(cont->GetObject(runNumber))){
2248 // AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2249 // AliError("V0 Calibration not fully set!\n");
2252 // fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2253 // fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2259 delete fpol0; fpol0=0;
2262 //_____________________________________________________________________________
2263 void AliAnalysisTaskPi0Flow::SetESDTrackCuts()
2266 // Create ESD track cut
2267 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
2268 //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
2269 fESDtrackCuts->SetRequireTPCRefit(kTRUE);
2273 //_____________________________________________________________________________
2274 void AliAnalysisTaskPi0Flow::SetGeometry()
2276 // Initialize the PHOS geometry
2277 if( kLHC10h == fPeriod && fEventESD ) {
2278 TGeoManager::Import("geometry.root"); //TODO: should perhaps not be done
2279 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
2281 AliError("geometry (fPHOSGeo) not initialised");
2286 AliOADBContainer geomContainer("phosGeo");
2287 geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
2288 TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
2289 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
2290 for(Int_t mod=0; mod<5; mod++) {
2291 if(!matrixes->At(mod)) {
2293 AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2297 fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
2299 AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2305 //_____________________________________________________________________________
2306 void AliAnalysisTaskPi0Flow::SetPHOSCalibData()
2308 if( fPHOSCalibData )
2309 delete fPHOSCalibData;
2312 // Calibration only needed for ESD
2313 if( fEventESD /*&& */ ) {
2314 if( kLHC10h == fPeriod && fEventESD ) {
2315 //We have to apply re-calibration for pass1 LCH10h
2316 // Initialize decalibration factors in the form of the OCDB object
2317 AliCDBManager * man = AliCDBManager::Instance();
2318 man->SetRun(140000) ; //TODO; revise, this should probably not b done.
2319 man->SetDefaultStorage("local://OCDB");
2321 fPHOSCalibData = new AliPHOSCalibData();
2325 //_____________________________________________________________________________
2326 void AliAnalysisTaskPi0Flow::SetVertex()
2328 const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
2329 if( primaryVertex ) {
2330 fVertex[0] = primaryVertex->GetX();
2331 fVertex[1] = primaryVertex->GetY();
2332 fVertex[2] = primaryVertex->GetZ();
2335 AliError("Event has 0x0 Primary Vertex, defaulting to origo");
2340 fVertexVector = TVector3(fVertex);
2341 FillHistogram("hZvertex", fVertexVector.z(), fInternalRunNumber-0.5);
2344 AliInfo(Form("Vertex is set to (%.1f,%.1f,%.1f)", fVertex[0], fVertex[1], fVertex[2]));
2346 fVtxBin=0 ;// No support for vtx binning implemented.
2349 //_____________________________________________________________________________
2350 Bool_t AliAnalysisTaskPi0Flow::RejectEventVertex()
2352 if( ! fEvent->GetPrimaryVertex() )
2353 return true; // reject
2354 LogSelection(1, fInternalRunNumber);
2356 if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ )
2357 return true; // reject
2358 LogSelection(2, fInternalRunNumber);
2360 return false; // accept event.
2363 //_____________________________________________________________________________
2364 void AliAnalysisTaskPi0Flow::SetCentrality()
2366 AliCentrality *centrality = fEvent->GetCentrality();
2368 fCentralityV0M=centrality->GetCentralityPercentile("V0M");
2370 AliError("Event has 0x0 centrality");
2371 fCentralityV0M = -1.;
2373 FillHistogram("hCentrality",fCentralityV0M,fInternalRunNumber-0.5) ;
2375 fCentBin = GetCentralityBin(fCentralityV0M);
2378 AliInfo(Form("Centrality (bin) is: %f (%d)", fCentralityV0M, fCentBin));
2381 //_____________________________________________________________________________
2382 Bool_t AliAnalysisTaskPi0Flow::RejectCentrality()
2384 if( ! fEvent->GetCentrality() )
2385 return true; // reject
2386 LogSelection(3, fInternalRunNumber);
2388 // if( fCentralityV0M <= 0. || fCentralityV0M>80. )
2389 // return true; // reject
2391 int lastBinUpperIndex = fCentEdges.GetSize() -1;
2392 if( fCentralityV0M > fCentEdges[lastBinUpperIndex] ) {
2394 AliInfo("Rejecting due to centrality outside of binning.");
2395 return true; // reject
2397 LogSelection(4, fInternalRunNumber);
2399 if( fCentralityV0M < fCentEdges[0] ) {
2401 AliInfo("Rejecting due to centrality outside of binning.");
2402 return true; // reject
2404 LogSelection(5, fInternalRunNumber);
2410 //_____________________________________________________________________________
2411 void AliAnalysisTaskPi0Flow::EvalReactionPlane()
2413 // assigns: fHaveTPCRP and fRP
2414 // also does a few histogram fills
2416 AliEventplane *eventPlane = fEvent->GetEventplane();
2417 if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2419 Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
2420 FillHistogram("phiRP",reactionPlaneQ,fCentralityV0M) ;
2422 if(reactionPlaneQ==999 || reactionPlaneQ < 0.){ //reaction plain was not defined
2423 if( fDebug ) AliInfo(Form("No Q Reaction Plane, value is %f", reactionPlaneQ));
2424 fHaveTPCRP = kFALSE;
2427 if( fDebug >= 2 ) AliInfo(Form("Q Reaction Plane is %f", reactionPlaneQ));
2432 fRP = ApplyFlattening(reactionPlaneQ, fCentralityV0M) ;
2434 while(fRP<0) fRP+=TMath::Pi();
2435 while(fRP>TMath::Pi()) fRP-=TMath::Pi();
2436 FillHistogram("phiRPflat",fRP,fCentralityV0M) ;
2437 Double_t dPsi = eventPlane->GetQsubRes() ;
2438 FillHistogram("cos2AC",TMath::Cos(2.*dPsi),fCentralityV0M) ;
2445 //____________________________________________________________________________
2446 void AliAnalysisTaskPi0Flow::EvalV0ReactionPlane(){
2447 // set: fRPV0A and fRPV0C
2449 // Do Manual V0 EP Calculation
2450 if ( fManualV0EPCalc )
2453 AliVVZERO* v0 = fEvent->GetVZEROData();
2455 //reset Q vector info
2456 Double_t Qxa2 = 0, Qya2 = 0;
2457 Double_t Qxc2 = 0, Qyc2 = 0;
2459 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
2460 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
2461 Float_t multv0 = v0->GetMultiplicity(iv0);
2462 if (iv0 < 32){ // V0C
2463 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2464 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2466 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2467 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2473 if(fCentralityV0M < 5) iC = 0;
2474 else if(fCentralityV0M < 10) iC = 1;
2475 else if(fCentralityV0M < 20) iC = 2;
2476 else if(fCentralityV0M < 30) iC = 3;
2477 else if(fCentralityV0M < 40) iC = 4;
2478 else if(fCentralityV0M < 50) iC = 5;
2479 else if(fCentralityV0M < 60) iC = 6;
2480 else if(fCentralityV0M < 70) iC = 7;
2483 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
2484 Double_t Qxamean2 = fMeanQ[iC][1][0];
2485 Double_t Qxarms2 = fWidthQ[iC][1][0];
2486 Double_t Qyamean2 = fMeanQ[iC][1][1];
2487 Double_t Qyarms2 = fWidthQ[iC][1][1];
2489 Double_t Qxcmean2 = fMeanQ[iC][0][0];
2490 Double_t Qxcrms2 = fWidthQ[iC][0][0];
2491 Double_t Qycmean2 = fMeanQ[iC][0][1];
2492 Double_t Qycrms2 = fWidthQ[iC][0][1];
2494 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
2495 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
2496 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
2497 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
2499 fRPV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
2500 fRPV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
2502 else // Use Official V0 EP Calculation.
2504 AliEventplane *eventPlane = fEvent->GetEventplane();
2505 if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2506 fRPV0A = eventPlane->GetEventplane("V0A", fEvent);
2507 fRPV0C = eventPlane->GetEventplane("V0C", fEvent);
2510 // Check that the A&C RP are within allowed range.
2511 if( fDebug >= 3 && (fRPV0A<0 || fRPV0A>TMath::Pi() ) )
2512 AliInfo(Form("RPV0A outside of permited range [0,pi]: %f, correcting", fRPV0A));
2513 if( fDebug >= 3 && (fRPV0C<0 || fRPV0C>TMath::Pi() ) )
2514 AliInfo(Form("RPV0C outside of permited range [0,pi]: %f, correcting", fRPV0C));
2515 while (fRPV0A<0 ) fRPV0A+=TMath::Pi() ;
2516 while (fRPV0A>TMath::Pi()) fRPV0A-=TMath::Pi() ;
2517 while (fRPV0C<0 ) fRPV0C+=TMath::Pi() ;
2518 while (fRPV0C>TMath::Pi()) fRPV0C-=TMath::Pi() ;
2520 // Reaction plane histograms before flattening
2522 AliInfo(Form("V0 Reaction Plane before flattening: A side: %f, C side: %f", fRPV0A, fRPV0C));
2524 FillHistogram("phiRPV0A" ,fRPV0A,fCentralityV0M);
2525 FillHistogram("phiRPV0C" ,fRPV0C,fCentralityV0M);
2526 FillHistogram("phiRPV0AC",fRPV0A,fRPV0C,fCentralityV0M) ;
2529 fRPV0A=ApplyFlatteningV0A(fRPV0A,fCentralityV0M) ;
2530 while (fRPV0A<0 ) fRPV0A+=TMath::Pi() ;
2531 while (fRPV0A>TMath::Pi()) fRPV0A-=TMath::Pi() ;
2533 fRPV0C=ApplyFlatteningV0C(fRPV0C,fCentralityV0M) ;
2534 while (fRPV0C<0 ) fRPV0C+=TMath::Pi() ;
2535 while (fRPV0C>TMath::Pi()) fRPV0C-=TMath::Pi() ;
2538 AliInfo(Form("V0 Reaction Plane after flattening: A side: %f, C side: %f", fRPV0A, fRPV0C));
2540 FillHistogram("phiRPV0Aflat",fRPV0A,fCentralityV0M) ;
2541 FillHistogram("cos2V0AC",TMath::Cos(2.*(fRPV0A-fRPV0C)),fCentralityV0M) ;
2543 FillHistogram("phiRPV0ATPC",fRP,fRPV0A,fCentralityV0M) ;
2544 FillHistogram("cos2V0ATPC",TMath::Cos(2.*(fRP-fRPV0A)),fCentralityV0M) ;
2547 FillHistogram("phiRPV0Cflat",fRPV0C,fCentralityV0M) ;
2549 FillHistogram("phiRPV0CTPC",fRP,fRPV0C,fCentralityV0M) ;
2550 FillHistogram("cos2V0CTPC",TMath::Cos(2.*(fRP-fRPV0C)),fCentralityV0M) ;
2553 //____________________________________________________________________________
2554 void AliAnalysisTaskPi0Flow::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20){
2555 //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
2557 const Double_t rCut=4.5 ;
2559 const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2560 // Calculates the center of gravity in the local PHOS-module coordinates
2562 const Int_t mulDigit=clu->GetNCells() ;
2563 Double_t xc[mulDigit] ;
2564 Double_t zc[mulDigit] ;
2565 Double_t wi[mulDigit] ;
2568 const Double_t logWeight=4.5 ;
2569 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2573 Int_t absId = clu->GetCellAbsId(iDigit) ;
2574 fPHOSGeo->AbsToRelNumbering(absId, relid) ;
2575 fPHOSGeo->RelPosInModule(relid, xi, zi);
2578 Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
2580 if (clu->E()>0 && ei>0) {
2581 wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
2582 Double_t w=wi[iDigit];
2583 x += xc[iDigit] * w ;
2584 z += zc[iDigit] * w ;
2597 Double_t xCut = 0. ;
2598 Double_t zCut = 0. ;
2599 for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2600 Double_t w=wi[iDigit];
2602 Double_t xi= xc[iDigit] ;
2603 Double_t zi= zc[iDigit] ;
2604 if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
2607 dxx += w * xi * xi ;
2608 dzz += w * zi * zi ;
2609 dxz += w * xi * zi ;
2621 dxx -= xCut * xCut ;
2622 dzz -= zCut * zCut ;
2623 dxz -= xCut * zCut ;
2625 m02 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
2626 m20 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
2633 //____________________________________________________________________________
2634 Bool_t AliAnalysisTaskPi0Flow::TestCoreLambda(Double_t pt,Double_t l1,Double_t l2){
2635 //Evaluates if lambdas correspond to photon cluster
2636 //Tuned using pp date
2637 //For core radius R=4.5
2638 Double_t l1Mean = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
2639 Double_t l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
2640 Double_t l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
2641 Double_t l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
2642 Double_t c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
2644 Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
2645 0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
2646 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
2647 return (R2<2.5*2.5) ;