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