]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODEventCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
q vector calibration for generated data
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraAODEventCuts.cxx
CommitLineData
c88234ad 1
2/**************************************************************************
3 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17//-----------------------------------------------------------------
18// AliSpectraAODEventCuts class
19//-----------------------------------------------------------------
20
21#include "TChain.h"
22#include "TTree.h"
23#include "TLegend.h"
24#include "TH1F.h"
829b5a81 25#include "TH1I.h"
c88234ad 26#include "TH2F.h"
27#include "TCanvas.h"
93db93de 28#include "TProfile.h"
7b364a00 29#include "TSpline.h"
c88234ad 30#include "AliAnalysisTask.h"
31#include "AliAnalysisManager.h"
32#include "AliAODTrack.h"
33#include "AliAODMCParticle.h"
34#include "AliAODEvent.h"
35#include "AliAODInputHandler.h"
36#include "AliAnalysisTaskESDfilter.h"
37#include "AliAnalysisDataContainer.h"
93db93de 38#include "AliESDVZERO.h"
39#include "AliAODVZERO.h"
c88234ad 40#include "AliSpectraAODEventCuts.h"
f6a38178 41#include "AliSpectraAODTrackCuts.h"
c88234ad 42#include <iostream>
43
44using namespace std;
45
46ClassImp(AliSpectraAODEventCuts)
47
93db93de 48AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) :
49 TNamed(name, "AOD Event Cuts"),
50 fAOD(0),
51 fSelectBit(AliVEvent::kMB),
52 fCentralityMethod("V0M"),
53 fTrackBits(1),
54 fIsMC(0),
48c4d28e 55 fIsLHC10h(1),
93db93de 56 fTrackCuts(0),
57 fIsSelected(0),
58 fCentralityCutMin(0.),
59 fCentralityCutMax(999),
60 fQVectorCutMin(-999.),
61 fQVectorCutMax(999.),
a162b3e2 62 fVertexCutMin(-10.),
63 fVertexCutMax(10.),
93db93de 64 fMultiplicityCutMin(-999.),
65 fMultiplicityCutMax(99999.),
48c4d28e 66 fqV0C(-999.),
67 fqV0A(-999.),
4fd2f0ad
LM		 68 fqV0Cx(-999.),
69 fqV0Ax(-999.),
70 fqV0Cy(-999.),
71 fqV0Ay(-999.),
fdbff5a1 72 fPsiV0C(-999.),
73 fPsiV0A(-999.),
48c4d28e 74 fCent(-999.),
93db93de 75 fOutput(0),
76 fCalib(0),
77 fRun(-1),
78 fMultV0(0),
79 fV0Cpol1(-1),
80 fV0Cpol2(-1),
81 fV0Cpol3(-1),
82 fV0Cpol4(-1),
83 fV0Apol1(-1),
84 fV0Apol2(-1),
85 fV0Apol3(-1),
7b364a00 86 fV0Apol4(-1),
87 fQvecIntList(0),
88 fQvecIntegral(0),
89 fSplineArrayV0A(0),
91083d79 90 fSplineArrayV0C(0),
cafe9ad1 91 fQgenIntegral(0),
92 fSplineArrayV0Agen(0),
93 fSplineArrayV0Cgen(0),
91083d79 94 fNch(0),
95 fQvecCalibType(0)
c88234ad 96{
f6a38178 97 // Constructor
93db93de 98 fOutput=new TList();
99 fOutput->SetOwner();
100 fOutput->SetName("fOutput");
101
102 fCalib=new TList();
103 fCalib->SetOwner();
104 fCalib->SetName("fCalib");
105
7b364a00 106 fQvecIntList=new TList();
107 fQvecIntList->SetOwner();
108 fQvecIntList->SetName("fQvecIntList");
109
93db93de 110 TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
23f51e78 111 TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
112 TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
113 TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
114 TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
115 TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
116 TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
117 TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
118 TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
119 TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
120 TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
121 TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
122 TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
123 TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
93db93de 124
7b364a00 125 fSplineArrayV0A = new TObjArray();
126 fSplineArrayV0A->SetOwner();
127 fSplineArrayV0C = new TObjArray();
128 fSplineArrayV0C->SetOwner();
cafe9ad1 129 fSplineArrayV0Agen = new TObjArray();
130 fSplineArrayV0Agen->SetOwner();
131 fSplineArrayV0Cgen = new TObjArray();
132 fSplineArrayV0Cgen->SetOwner();
7b364a00 133
93db93de 134 fOutput->Add(fHistoCuts);
135 fOutput->Add(fHistoVtxBefSel);
136 fOutput->Add(fHistoVtxAftSel);
137 fOutput->Add(fHistoEtaBefSel);
138 fOutput->Add(fHistoEtaAftSel);
139 fOutput->Add(fHistoNChAftSel);
140 fOutput->Add(fHistoQVector);
141 fOutput->Add(fHistoEP);
142 fOutput->Add(fPsiACor);
143 fOutput->Add(fPsiCCor);
144 fOutput->Add(fQVecACor);
145 fOutput->Add(fQVecCCor);
23f51e78 146 fOutput->Add(fV0M);
147 fOutput->Add(fV0MCor);
93db93de 148
149 for (Int_t i = 0; i<10; i++){
150 fMeanQxa2[i] = -1;
151 fMeanQya2[i] = -1;
152 fMeanQxc2[i] = -1;
153 fMeanQyc2[i] = -1;
154 }
c88234ad 155}
156
157//______________________________________________________
f6a38178 158Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts *trackcuts)
c88234ad 159{
f6a38178 160 // Returns true if Event Cuts are selected and applied
161 fAOD = aod;
31579941 162 fTrackCuts = trackcuts; // FIXME: if track cuts is 0, do not set (use the pre-set member). Do we need to pass this here at all??
163 // FIXME: all those references by name are slow.
93db93de 164 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
165 Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
decf69d9 166 if(!IsPhysSel)return IsPhysSel;
93db93de 167 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
f6a38178 168 //loop on tracks, before event selection, filling QA histos
169 AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
93db93de 170 if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
decf69d9 171 fIsSelected =kFALSE;
10a8ccbe 172 if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
00493191 173 fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
decf69d9 174 }
f6a38178 175 if(fIsSelected){
93db93de 176 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
177 if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
f6a38178 178 }
ae0fdd7d 179 Int_t Nch=0;
f6a38178 180 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
181 AliAODTrack* track = fAOD->GetTrack(iTracks);
decf69d9 182 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
93db93de 183 ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
ae0fdd7d 184 if(fIsSelected){
93db93de 185 ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
ae0fdd7d 186 Nch++;
187 }
f6a38178 188 }
91083d79 189 fNch = Nch;
93db93de 190 if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
f6a38178 191 return fIsSelected;
c88234ad 192}
193
194//______________________________________________________
195Bool_t AliSpectraAODEventCuts::CheckVtxRange()
196{
197 // reject events outside of range
f6a38178 198 AliAODVertex * vertex = fAOD->GetPrimaryVertex();
264de30e 199 //when moving to 2011 wìone has to add a cut using SPD vertex.
200 //The point is that for events with |z|>20 the vertexer tracks is not working (only 2011!). One has to put a safety cut using SPD vertex large e.g. 15cm
f6a38178 201 if (!vertex)
202 {
93db93de 203 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
c88234ad 204 return kFALSE;
f6a38178 205 }
10a8ccbe 206 if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
f6a38178 207 {
c88234ad 208 return kTRUE;
f6a38178 209 }
93db93de 210 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
f6a38178 211 return kFALSE;
c88234ad 212}
213
214//______________________________________________________
215Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
216{
f6a38178 217 // Check centrality cut
48c4d28e 218 fCent=-999.;
219 fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
220 if ( (fCent <= fCentralityCutMax) && (fCent >= fCentralityCutMin) ) return kTRUE;
93db93de 221 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
f6a38178 222 return kFALSE;
c88234ad 223}
224
10a8ccbe 225//______________________________________________________
226Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
227{
228 // Check multiplicity cut
31579941 229 // FIXME: why this is not tracket in the track stats histos?
10a8ccbe 230 Int_t Ncharged=0;
231 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
232 AliAODTrack* track = fAOD->GetTrack(iTracks);
233 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
234 Ncharged++;
235 }
10a8ccbe 236 if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
237
238 return kFALSE;
239}
240
decf69d9 241//______________________________________________________
242Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
23f51e78 243{
fdbff5a1 244 Double_t qVZERO = -999.;
245 Double_t psi=-999.;
246
247 if(fIsLHC10h){
248 qVZERO=CalculateQVectorLHC10h();
249 psi=fPsiV0A;
250 }else{
2e20a62a 251 qVZERO=CalculateQVector();
252 psi=fPsiV0A;
fdbff5a1 253 }
254
255 //cut on q vector
16ee3540 256 if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
23f51e78 257 Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
258 ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
259 ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
93db93de 260 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
23f51e78 261
fdbff5a1 262
decf69d9 263 return kTRUE;
decf69d9 264}
265
93db93de 266//______________________________________________________
93db93de 267Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
268
269 Int_t run = fAOD->GetRunNumber();
270 if(run != fRun){
271 // Load the calibrations run dependent
272 if(OpenInfoCalbration(run))fRun=run;
48c4d28e 273 else{
274 fqV0C=-999.;
275 fqV0A=-999.;
276 return -999.;
277 }
93db93de 278 }
279
280 //V0 info
281 Double_t Qxa2 = 0, Qya2 = 0;
282 Double_t Qxc2 = 0, Qyc2 = 0;
283 Double_t sumMc = 0, sumMa = 0;
284
285 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
286
287 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
288
289 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
290
291 Float_t multv0 = aodV0->GetMultiplicity(iv0);
23f51e78 292 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
293
93db93de 294 if (iv0 < 32){
295
296 Double_t multCorC = -10;
297
298 if (iv0 < 8)
299 multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
300 else if (iv0 >= 8 && iv0 < 16)
301 multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
302 else if (iv0 >= 16 && iv0 < 24)
303 multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
304 else if (iv0 >= 24 && iv0 < 32)
305 multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
306
307 if (multCorC < 0){
308 cout<<"Problem with multiplicity in V0C"<<endl;
48c4d28e 309 fqV0C=-999.;
310 fqV0A=-999.;
311 return -999.;
93db93de 312 }
313
314 Qxc2 += TMath::Cos(2*phiV0) * multCorC;
315 Qyc2 += TMath::Sin(2*phiV0) * multCorC;
316
317 sumMc += multCorC;
23f51e78 318 ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
93db93de 319
320 } else {
321
322 Double_t multCorA = -10;
323
324 if (iv0 >= 32 && iv0 < 40)
325 multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
326 else if (iv0 >= 40 && iv0 < 48)
327 multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
328 else if (iv0 >= 48 && iv0 < 56)
329 multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
330 else if (iv0 >= 56 && iv0 < 64)
331 multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
332
333 if (multCorA < 0){
334 cout<<"Problem with multiplicity in V0A"<<endl;
48c4d28e 335 fqV0C=-999.;
336 fqV0A=-999.;
337 return -999.;
93db93de 338 }
339
340 Qxa2 += TMath::Cos(2*phiV0) * multCorA;
341 Qya2 += TMath::Sin(2*phiV0) * multCorA;
342
343 sumMa += multCorA;
23f51e78 344 ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
93db93de 345 }
346 }
347
348 Short_t centrV0 = GetCentrCode(fAOD);
349
6a6d2363 350 Double_t Qxamean2 = 0.;
351 Double_t Qyamean2 = 0.;
352 Double_t Qxcmean2 = 0.;
353 Double_t Qycmean2 = 0.;
354
355 if(centrV0!=-1){
356 Qxamean2 = fMeanQxa2[centrV0];
357 Qyamean2 = fMeanQya2[centrV0];
358 Qxcmean2 = fMeanQxc2[centrV0];
359 Qycmean2 = fMeanQyc2[centrV0];
360 }
93db93de 361
362 Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
363 Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
364 Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
365 Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
366
fdbff5a1 367 fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
368 fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
93db93de 369
fdbff5a1 370 ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
371 ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
93db93de 372
48c4d28e 373 fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
fdbff5a1 374 fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
4fd2f0ad
LM		 375 fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
376 fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
377 fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
378 fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
48c4d28e 379
380 ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
381 ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
93db93de 382
fdbff5a1 383 return fqV0A; //FIXME we have to combine VZERO-A and C
93db93de 384}
385
2e20a62a 386//______________________________________________________
387Double_t AliSpectraAODEventCuts::CalculateQVector(){
388
389 //V0 info
390 Double_t Qxa2 = 0, Qya2 = 0;
391 Double_t Qxc2 = 0, Qyc2 = 0;
2e20a62a 392
393 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
394
395 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
396
397 Float_t multv0 = aodV0->GetMultiplicity(iv0);
16e0a56f 398
2e20a62a 399 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
400
2e20a62a 401 }
402
403 fPsiV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,8,2,Qxa2,Qya2); // V0A
404 fPsiV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,9,2,Qxc2,Qyc2); // V0C
405
406 ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
407 ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
408
16e0a56f 409 fqV0A = TMath::Sqrt((Qxa2*Qxa2 + Qya2*Qya2));
410 fqV0C = TMath::Sqrt((Qxc2*Qxc2 + Qyc2*Qyc2));
411 fqV0Ax = Qxa2;
412 fqV0Cx = Qxc2;
413 fqV0Ay = Qya2;
414 fqV0Cy = Qyc2;
2e20a62a 415
416 ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
417 ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
418
419 return fqV0A; //FIXME we have to combine VZERO-A and C
420
421}
422
fdbff5a1 423//______________________________________________________
93db93de 424Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
425{
426
427 Short_t centrCode = -1;
428
429 AliCentrality* centrality = 0;
430 AliAODEvent* aod = (AliAODEvent*)ev;
431 centrality = aod->GetHeader()->GetCentralityP();
432
433 Float_t centV0 = centrality->GetCentralityPercentile("V0M");
434 Float_t centTrk = centrality->GetCentralityPercentile("TRK");
435
436
437 if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){
438
439 if ((centV0 > 0) && (centV0 <= 5.0))
440 centrCode = 0;
441 else if ((centV0 > 5.0) && (centV0 <= 10.0))
442 centrCode = 1;
443 else if ((centV0 > 10.0) && (centV0 <= 20.0))
444 centrCode = 2;
445 else if ((centV0 > 20.0) && (centV0 <= 30.0))
446 centrCode = 3;
447 else if ((centV0 > 30.0) && (centV0 <= 40.0))
448 centrCode = 4;
449 else if ((centV0 > 40.0) && (centV0 <= 50.0))
450 centrCode = 5;
451 else if ((centV0 > 50.0) && (centV0 <= 60.0))
452 centrCode = 6;
453 else if ((centV0 > 60.0) && (centV0 <= 70.0))
454 centrCode = 7;
455 else if ((centV0 > 70.0) && (centV0 <= 80.0))
456 centrCode = 8;
457 }
458
459 return centrCode;
460
461}
462
c88234ad 463//______________________________________________________
464void AliSpectraAODEventCuts::PrintCuts()
465{
f6a38178 466 // print info about event cuts
467 cout << "Event Stats" << endl;
93db93de 468 cout << " > Trigger Selection: " << fSelectBit << endl;
469 cout << " > Centrality estimator: " << fCentralityMethod << endl;
470 cout << " > Number of accepted events: " << NumberOfEvents() << endl;
471 cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
472 cout << " > Number of PhysSel events: " << NumberOfPhysSelEvents() << endl;
473 cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
474 cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
475 cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
476 cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
829b5a81 477 cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
16ee3540 478 cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
10a8ccbe 479 cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
480 cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
829b5a81 481 cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
f6a38178 482}
c88234ad 483
93db93de 484//_____________________________________________________________________________
485Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
972a21ad 486{
93db93de 487
488 AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
489 if(!cont){
490 printf("OADB object hMultV0BefCorr is not available in the file\n");
491 return 0;
492 }
493
494 if(!(cont->GetObject(run))){
495 printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
496 return 0;
497 }
498 fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
499
500 TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7);
501 fMultV0->Fit(fpolc1, "RN");
502 fV0Cpol1 = fpolc1->GetParameter(0);
503
504 TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15);
505 fMultV0->Fit(fpolc2, "RN");
506 fV0Cpol2 = fpolc2->GetParameter(0);
507
508 TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23);
509 fMultV0->Fit(fpolc3, "RN");
510 fV0Cpol3 = fpolc3->GetParameter(0);
511
512 TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31);
513 fMultV0->Fit(fpolc4, "RN");
514 fV0Cpol4 = fpolc4->GetParameter(0);
515
516 TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
6a6d2363 517 fMultV0->Fit(fpola1, "RN");
93db93de 518 fV0Apol1 = fpola1->GetParameter(0);
519
520 TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
521 fMultV0->Fit(fpola2, "RN");
522 fV0Apol2 = fpola2->GetParameter(0);
a162b3e2 523
93db93de 524 TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
525 fMultV0->Fit(fpola3, "RN");
526 fV0Apol3 = fpola3->GetParameter(0);
527
528 TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
529 fMultV0->Fit(fpola4, "RN");
530 fV0Apol4 = fpola4->GetParameter(0);
531
93db93de 532 for(Int_t i=0; i < 10; i++){
533
534 char nameQxa2[100];
6a6d2363 535 snprintf(nameQxa2,100, "hQxa2m_%i", i);
93db93de 536
537 char nameQya2[100];
6a6d2363 538 snprintf(nameQya2,100, "hQya2m_%i", i);
93db93de 539
540 char nameQxc2[100];
6a6d2363 541 snprintf(nameQxc2,100, "hQxc2m_%i", i);
93db93de 542
543 char nameQyc2[100];
6a6d2363 544 snprintf(nameQyc2,100, "hQyc2m_%i", i);
93db93de 545
546 AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
547 if(!contQxa2){
548 printf("OADB object %s is not available in the file\n", nameQxa2);
549 return 0;
550 }
551
552 if(!(contQxa2->GetObject(run))){
553 printf("OADB object %s is not available for run %i\n", nameQxa2, run);
554 return 0;
555 }
556
557 fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
558
559
560 AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
561 if(!contQya2){
562 printf("OADB object %s is not available in the file\n", nameQya2);
563 return 0;
564 }
565
566 if(!(contQya2->GetObject(run))){
567 printf("OADB object %s is not available for run %i\n", nameQya2, run);
568 return 0;
569 }
570
571 fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
572
573
574 AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
575 if(!contQxc2){
576 printf("OADB object %s is not available in the file\n", nameQxc2);
577 return 0;
578 }
579
580 if(!(contQxc2->GetObject(run))){
581 printf("OADB object %s is not available for run %i\n", nameQxc2, run);
582 return 0;
583 }
584
585 fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
586
587
588 AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
589 if(!contQyc2){
590 printf("OADB object %s is not available in the file\n", nameQyc2);
591 return 0;
592 }
593
594 if(!(contQyc2->GetObject(run))){
595 printf("OADB object %s is not available for run %i\n", nameQyc2, run);
596 return 0;
597 }
598
599 fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
600
601 }
602 return 1;
972a21ad 603}
604
605//______________________________________________________
606
607
c88234ad 608Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
609{
610 // Merge a list of AliSpectraAODEventCuts objects with this.
611 // Returns the number of merged objects (including this).
93db93de 612
613 AliInfo("Merging");
614
c88234ad 615 if (!list)
616 return 0;
617
618 if (list->IsEmpty())
619 return 1;
93db93de 620
c88234ad 621 TIterator* iter = list->MakeIterator();
622 TObject* obj;
93db93de 623
c88234ad 624 // collections of all histograms
93db93de 625 TList collections;
626
c88234ad 627 Int_t count = 0;
628
629 while ((obj = iter->Next())) {
630 AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
631 if (entry == 0)
632 continue;
633
93db93de 634 TList * l = entry->GetOutputList();
635 collections.Add(l);
c88234ad 636 count++;
637 }
638
93db93de 639 fOutput->Merge(&collections);
c88234ad 640
641 delete iter;
642
643 return count+1;
644}
645
7b364a00 646//______________________________________________________
647Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
648
649 //check Qvec and Centrality consistency
650 if(fCent>90.) return -999.;
651 if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
652 if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
653
654 fQvecIntegral = 0x0;
655
656 if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
657 if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
658
91083d79 659
660 Int_t ic = -999;
661
662 if(fQvecCalibType==1){
663 if(fNch<0.) return -999.;
664 ic = GetNchBin();
665 } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
7b364a00 666
667 TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
668
669 TSpline *spline = 0x0;
670
671 if(v0side==0/*V0A*/){
91083d79 672 if( CheckSplineArray(fSplineArrayV0A, ic) ) {
7b364a00 673 spline = (TSpline*)fSplineArrayV0A->At(ic);
d7696a14 674 //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
7b364a00 675 } else {
676 spline = new TSpline3(h1D,"sp3");
677 fSplineArrayV0A->AddAtAndExpand(spline,ic);
d7696a14 678 //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
7b364a00 679 }
680 }
681 else if(v0side==1/*V0C*/){
91083d79 682 if( CheckSplineArray(fSplineArrayV0C, ic) ) {
7b364a00 683 spline = (TSpline*)fSplineArrayV0C->At(ic);
684 } else {
685 spline = new TSpline3(h1D,"sp3");
686 fSplineArrayV0C->AddAtAndExpand(spline,ic);
687 }
688 }
689
690 Double_t percentile = -999.;
691 if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
692 if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
693
694 if(percentile>100. || percentile<0.) return -999.;
695
696 return percentile;
697}
698
ed33b9ef 699//______________________________________________________
700Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side){
701
702 if(!fIsMC) return -999.;
703
704 // 1. get MC array
705 TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
706
707 if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
708
709 Double_t Qx2mc = 0., Qy2mc = 0.;
710 Int_t mult2mc = 0;
711
712 Int_t nMC = arrayMC->GetEntries();
713
714 // 2. loop on generated
715 for (Int_t iMC = 0; iMC < nMC; iMC++){
716 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
717
718 // 3. set VZERO side - FIXME Add TPC!
719 Double_t EtaVZERO[2] = {-999.,-999.};
720 if(v0side==0/*V0A*/){ EtaVZERO[0] = 2.8; EtaVZERO[1] = 5.1; }
721 if(v0side==1/*V0C*/){ EtaVZERO[0] = -3.7; EtaVZERO[1] = -1.7; }
722
723 if(partMC->Eta()<EtaVZERO[0] || partMC->Eta()>EtaVZERO[1]) continue;
724
725 // 4. Calculate Qvec components
726
727 Qx2mc += TMath::Cos(2.*partMC->Phi());
728 Qy2mc += TMath::Sin(2.*partMC->Phi());
729 mult2mc++;
730
731 }
732
733 // 5. return q vector
734 return TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
735
736}
737
cafe9ad1 738//______________________________________________________
739Double_t AliSpectraAODEventCuts::GetQvecPercentileMC(Int_t v0side){
740
741 //check Qvec and Centrality consistency
742 if(fCent>90.) return -999.;
743
744 Double_t qvec = CalculateQVectorMC(v0side);
745 if(qvec==-999.) return -999.;
746
747 fQgenIntegral = 0x0;
748
749 if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen"); }
750 if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen"); }
751 //FIXME you need a check on the TH2D, add AliFatal or a return.
752
753 Int_t ic = -999;
754
755 if(fQvecCalibType==1){
756 if(fNch<0.) return -999.;
757 ic = GetNchBin();
758 } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
759
760 TH1D *h1D = (TH1D*)fQgenIntegral->ProjectionY("h1Dgen",ic+1,ic+1);
761
762 TSpline *spline = 0x0;
763
764 if(v0side==0/*V0A*/){
765 if( CheckSplineArray(fSplineArrayV0Agen, ic) ) {
766 spline = (TSpline*)fSplineArrayV0Agen->At(ic);
767 //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
768 } else {
769 spline = new TSpline3(h1D,"sp3");
770 fSplineArrayV0Agen->AddAtAndExpand(spline,ic);
771 }
772 }
773 else if(v0side==1/*V0C*/){
774 if( CheckSplineArray(fSplineArrayV0Cgen, ic) ) {
775 spline = (TSpline*)fSplineArrayV0Cgen->At(ic);
776 } else {
777 spline = new TSpline3(h1D,"sp3");
778 fSplineArrayV0Cgen->AddAtAndExpand(spline,ic);
779 }
780 }
781
782 Double_t percentile = 100*spline->Eval(qvec);
783
784 if(percentile>100. || percentile<0.) return -999.;
785
786 return percentile;
787
788}
789
7b364a00 790//______________________________________________________
91083d79 791Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr, Int_t n){
7b364a00 792 //check TSpline array consistency
793
91083d79 794// Int_t n = (Int_t)fCent; //FIXME should be ok for icentr and nch
7b364a00 795
796 if(splarr->GetSize()<n) return kFALSE;
797
798 if(!splarr->At(n)) return kFALSE;
799
800 return kTRUE;
801
802}
91083d79 803
804//______________________________________________________
805Int_t AliSpectraAODEventCuts::GetNchBin(){
806
807 Double_t xmax = fQvecIntegral->GetXaxis()->GetXmax();
808
809 if(fNch>xmax) return fQvecIntegral->GetNbinsX();
810
811 return (fNch*fQvecIntegral->GetNbinsX())/fQvecIntegral->GetXaxis()->GetXmax();
812
813}
814