]> 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
TPC percentile calibration
[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) :
a72cae45 49TNamed(name, "AOD Event Cuts"),
50fAOD(0),
51fSelectBit(AliVEvent::kMB),
52fCentralityMethod("V0M"),
53fTrackBits(1),
54fIsMC(0),
55fIsLHC10h(1),
56fTrackCuts(0),
57fIsSelected(0),
58fCentralityCutMin(0.),
59fCentralityCutMax(999),
60fQVectorCutMin(-999.),
61fQVectorCutMax(999.),
62fVertexCutMin(-10.),
63fVertexCutMax(10.),
64fMultiplicityCutMin(-999.),
65fMultiplicityCutMax(99999.),
66fqTPC(-999.),
67fqV0C(-999.),
68fqV0A(-999.),
69fqV0Cx(-999.),
70fqV0Ax(-999.),
71fqV0Cy(-999.),
72fqV0Ay(-999.),
73fPsiV0C(-999.),
74fPsiV0A(-999.),
75fCent(-999.),
76fOutput(0),
77fCalib(0),
78fRun(-1),
79fMultV0(0),
80fV0Cpol1(-1),
81fV0Cpol2(-1),
82fV0Cpol3(-1),
83fV0Cpol4(-1),
84fV0Apol1(-1),
85fV0Apol2(-1),
86fV0Apol3(-1),
87fV0Apol4(-1),
88fQvecIntList(0),
89fQvecIntegral(0),
90fSplineArrayV0A(0),
91fSplineArrayV0C(0),
92fSplineArrayTPC(0),
93fQgenIntegral(0),
94fSplineArrayV0Agen(0),
95fSplineArrayV0Cgen(0),
96fQvecMC(0),
97fNch(0),
98fQvecCalibType(0),
99fV0Aeff(0)
c88234ad 100{
a72cae45 101 // Constructor
102 fOutput=new TList();
103 fOutput->SetOwner();
104 fOutput->SetName("fOutput");
105
106 fCalib=new TList();
107 fCalib->SetOwner();
108 fCalib->SetName("fCalib");
109
110 fQvecIntList=new TList();
111 fQvecIntList->SetOwner();
112 fQvecIntList->SetName("fQvecIntList");
113
114 TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
115 TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
116 TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
117 TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
118 TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
119 TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
120 TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
121 TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
122 TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
123 TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
124 TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
125 TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
126 TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
127 TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
128 TH2F *fV0Mmc = new TH2F("fV0Mmc", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
129
130 fSplineArrayV0A = new TObjArray();
131 fSplineArrayV0A->SetOwner();
132 fSplineArrayV0C = new TObjArray();
133 fSplineArrayV0C->SetOwner();
134 fSplineArrayTPC = new TObjArray();
135 fSplineArrayTPC->SetOwner();
136 fSplineArrayV0Agen = new TObjArray();
137 fSplineArrayV0Agen->SetOwner();
138 fSplineArrayV0Cgen = new TObjArray();
139 fSplineArrayV0Cgen->SetOwner();
140
141 fOutput->Add(fHistoCuts);
142 fOutput->Add(fHistoVtxBefSel);
143 fOutput->Add(fHistoVtxAftSel);
144 fOutput->Add(fHistoEtaBefSel);
145 fOutput->Add(fHistoEtaAftSel);
146 fOutput->Add(fHistoNChAftSel);
147 fOutput->Add(fHistoQVector);
148 fOutput->Add(fHistoEP);
149 fOutput->Add(fPsiACor);
150 fOutput->Add(fPsiCCor);
151 fOutput->Add(fQVecACor);
152 fOutput->Add(fQVecCCor);
153 fOutput->Add(fV0M);
154 fOutput->Add(fV0MCor);
155 fOutput->Add(fV0Mmc);
156
157 for (Int_t i = 0; i<10; i++){
158 fMeanQxa2[i] = -1;
159 fMeanQya2[i] = -1;
160 fMeanQxc2[i] = -1;
161 fMeanQyc2[i] = -1;
162 }
c88234ad 163}
164
165//______________________________________________________
f6a38178 166Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts *trackcuts)
c88234ad 167{
a72cae45 168 // Returns true if Event Cuts are selected and applied
169 fAOD = aod;
170 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??
171 // FIXME: all those references by name are slow.
172 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
173 Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
174 if(!IsPhysSel)return IsPhysSel;
175 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
176 //loop on tracks, before event selection, filling QA histos
177 AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
178 if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
179 fIsSelected =kFALSE;
180 if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
181 fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
182 }
183 if(fIsSelected){
184 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
185 if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
186 }
187 Int_t Nch=0;
188 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
0a918d8d 189 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
190 if(!track) AliFatal("Not a standard AOD");
191 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
a72cae45 192 ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
193 if(fIsSelected){
194 ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
195 Nch++;
196 }
197 }
198 fNch = Nch;
199 if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
200 return fIsSelected;
c88234ad 201}
202
203//______________________________________________________
204Bool_t AliSpectraAODEventCuts::CheckVtxRange()
205{
a72cae45 206 // reject events outside of range
207 AliAODVertex * vertex = fAOD->GetPrimaryVertex();
208 //when moving to 2011 wìone has to add a cut using SPD vertex.
209 //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
210 if (!vertex)
211 {
212 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
213 return kFALSE;
214 }
215 if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
216 {
217 return kTRUE;
218 }
219 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
220 return kFALSE;
c88234ad 221}
222
223//______________________________________________________
224Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
225{
a72cae45 226 // Check centrality cut
227 fCent=-999.;
228 fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
229 if ( (fCent <= fCentralityCutMax) && (fCent >= fCentralityCutMin) ) return kTRUE;
230 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
231 return kFALSE;
c88234ad 232}
233
10a8ccbe 234//______________________________________________________
235Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
236{
a72cae45 237 // Check multiplicity cut
238 // FIXME: why this is not tracket in the track stats histos?
239 Int_t Ncharged=0;
240 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
0a918d8d 241 AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
242 if(!track) AliFatal("Not a standard AOD");
a72cae45 243 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
244 Ncharged++;
245 }
246 if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
247
248 return kFALSE;
10a8ccbe 249}
250
decf69d9 251//______________________________________________________
a72cae45 252
decf69d9 253Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
23f51e78 254{
a72cae45 255 Double_t qVZERO = -999.;
256 Double_t psi=-999.;
257
258 if(fIsLHC10h){
259 qVZERO=CalculateQVectorLHC10h();
260 psi=fPsiV0A;
261 }else{
262 qVZERO=CalculateQVector();
263 psi=fPsiV0A;
264 }
265 //q vector from TPC
266 fqTPC=CalculateQVectorTPC();
267
268 //cut on q vector
269 if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
270 Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
271 ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
272 ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
273 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
274
275
276 return kTRUE;
decf69d9 277}
278
93db93de 279//______________________________________________________
93db93de 280Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
a72cae45 281
282 Int_t run = fAOD->GetRunNumber();
283 if(run != fRun){
284 // Load the calibrations run dependent
285 if(OpenInfoCalbration(run))fRun=run;
286 else{
287 fqV0C=-999.;
288 fqV0A=-999.;
289 return -999.;
290 }
291 }
292
293 //V0 info
294 Double_t Qxa2 = 0, Qya2 = 0;
295 Double_t Qxc2 = 0, Qyc2 = 0;
296 Double_t sumMc = 0, sumMa = 0;
297
298 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
299
300 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
301
302 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
303
304 Float_t multv0 = aodV0->GetMultiplicity(iv0);
305 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
306
307 if (iv0 < 32){
308
309 Double_t multCorC = -10;
310
311 if (iv0 < 8)
312 multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
313 else if (iv0 >= 8 && iv0 < 16)
314 multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
315 else if (iv0 >= 16 && iv0 < 24)
316 multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
317 else if (iv0 >= 24 && iv0 < 32)
318 multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
319
320 if (multCorC < 0){
321 cout<<"Problem with multiplicity in V0C"<<endl;
322 fqV0C=-999.;
323 fqV0A=-999.;
324 return -999.;
325 }
326
327 Qxc2 += TMath::Cos(2*phiV0) * multCorC;
328 Qyc2 += TMath::Sin(2*phiV0) * multCorC;
329
330 sumMc += multCorC;
331 ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
332
333 } else {
334
335 Double_t multCorA = -10;
336
337 if (iv0 >= 32 && iv0 < 40)
338 multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
339 else if (iv0 >= 40 && iv0 < 48)
340 multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
341 else if (iv0 >= 48 && iv0 < 56)
342 multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
343 else if (iv0 >= 56 && iv0 < 64)
344 multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
345
346 if (multCorA < 0){
347 cout<<"Problem with multiplicity in V0A"<<endl;
348 fqV0C=-999.;
349 fqV0A=-999.;
350 return -999.;
351 }
352
353 Qxa2 += TMath::Cos(2*phiV0) * multCorA;
354 Qya2 += TMath::Sin(2*phiV0) * multCorA;
355
356 sumMa += multCorA;
357 ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
358 }
359 }
360
361 Short_t centrV0 = GetCentrCode(fAOD);
362
363 Double_t Qxamean2 = 0.;
364 Double_t Qyamean2 = 0.;
365 Double_t Qxcmean2 = 0.;
366 Double_t Qycmean2 = 0.;
367
368 if(centrV0!=-1){
369 Qxamean2 = fMeanQxa2[centrV0];
370 Qyamean2 = fMeanQya2[centrV0];
371 Qxcmean2 = fMeanQxc2[centrV0];
372 Qycmean2 = fMeanQyc2[centrV0];
373 }
374
375 Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
376 Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
377 Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
378 Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
379
380 fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
381 fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
382
383 ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
384 ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
385
386 fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
387 fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
388 fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
389 fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
390 fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
391 fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
392
393 ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
394 ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
395
396 return fqV0A; //FIXME we have to combine VZERO-A and C
397}
398
399//______________________________________________________
400
401Double_t AliSpectraAODEventCuts::CalculateQVectorTPC(Double_t etaMin,Double_t etaMax){
402
403 Double_t Qx2 = 0, Qy2 = 0;
404 Int_t mult = 0;
405 for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
0a918d8d 406 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iT));
407 if(!aodTrack) AliFatal("Not a standard AOD");
a72cae45 408 if (!aodTrack->TestFilterBit(128)) continue; //FIXME track type hard coded -> TPC only constrained to the vertex
409 if (aodTrack->Eta() <etaMin || aodTrack->Eta() > etaMax)continue;
410 mult++;
411 Qx2 += TMath::Cos(2*aodTrack->Phi());
412 Qy2 += TMath::Sin(2*aodTrack->Phi());
413 }
414 if(mult!=0)fqTPC= TMath::Sqrt((Qx2*Qx2 + Qy2*Qy2)/mult);
415 return fqTPC;
93db93de 416}
417
2e20a62a 418//______________________________________________________
a72cae45 419
2e20a62a 420Double_t AliSpectraAODEventCuts::CalculateQVector(){
a72cae45 421
422 //V0 info
423 Double_t Qxa2 = 0, Qya2 = 0;
424 Double_t Qxc2 = 0, Qyc2 = 0;
425
426 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
427
428 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
429
430 Float_t multv0 = aodV0->GetMultiplicity(iv0);
431
432 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
433
434 }
435
436 fPsiV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,8,2,Qxa2,Qya2); // V0A
437 fPsiV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,9,2,Qxc2,Qyc2); // V0C
438
439 ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
440 ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
441
442 fqV0A = TMath::Sqrt((Qxa2*Qxa2 + Qya2*Qya2));
443 fqV0C = TMath::Sqrt((Qxc2*Qxc2 + Qyc2*Qyc2));
444 fqV0Ax = Qxa2;
445 fqV0Cx = Qxc2;
446 fqV0Ay = Qya2;
447 fqV0Cy = Qyc2;
448
449 ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
450 ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
451
452 return fqV0A; //FIXME we have to combine VZERO-A and C
453
2e20a62a 454}
455
fdbff5a1 456//______________________________________________________
93db93de 457Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
458{
a72cae45 459
460 Short_t centrCode = -1;
461
462 AliCentrality* centrality = 0;
463 AliAODEvent* aod = (AliAODEvent*)ev;
0a918d8d 464 centrality = ((AliVAODHeader*)aod->GetHeader())->GetCentralityP();
a72cae45 465
466 Float_t centV0 = centrality->GetCentralityPercentile("V0M");
467 Float_t centTrk = centrality->GetCentralityPercentile("TRK");
468
469
470 if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){
471
472 if ((centV0 > 0) && (centV0 <= 5.0))
473 centrCode = 0;
474 else if ((centV0 > 5.0) && (centV0 <= 10.0))
475 centrCode = 1;
476 else if ((centV0 > 10.0) && (centV0 <= 20.0))
477 centrCode = 2;
478 else if ((centV0 > 20.0) && (centV0 <= 30.0))
479 centrCode = 3;
480 else if ((centV0 > 30.0) && (centV0 <= 40.0))
481 centrCode = 4;
482 else if ((centV0 > 40.0) && (centV0 <= 50.0))
483 centrCode = 5;
484 else if ((centV0 > 50.0) && (centV0 <= 60.0))
485 centrCode = 6;
486 else if ((centV0 > 60.0) && (centV0 <= 70.0))
487 centrCode = 7;
488 else if ((centV0 > 70.0) && (centV0 <= 80.0))
489 centrCode = 8;
490 }
491
492 return centrCode;
493
93db93de 494}
495
c88234ad 496//______________________________________________________
497void AliSpectraAODEventCuts::PrintCuts()
498{
a72cae45 499 // print info about event cuts
500 cout << "Event Stats" << endl;
501 cout << " > Trigger Selection: " << fSelectBit << endl;
502 cout << " > Centrality estimator: " << fCentralityMethod << endl;
503 cout << " > Number of accepted events: " << NumberOfEvents() << endl;
504 cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
505 cout << " > Number of PhysSel events: " << NumberOfPhysSelEvents() << endl;
506 cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
507 cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
508 cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
509 cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
510 cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
511 cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
512 cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
513 cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
514 cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
f6a38178 515}
c88234ad 516
93db93de 517//_____________________________________________________________________________
518Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
972a21ad 519{
a72cae45 520
521 AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
522 if(!cont){
523 printf("OADB object hMultV0BefCorr is not available in the file\n");
524 return 0;
525 }
526
527 if(!(cont->GetObject(run))){
528 printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
529 return 0;
530 }
531 fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
532
533 TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7);
534 fMultV0->Fit(fpolc1, "RN");
535 fV0Cpol1 = fpolc1->GetParameter(0);
536
537 TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15);
538 fMultV0->Fit(fpolc2, "RN");
539 fV0Cpol2 = fpolc2->GetParameter(0);
540
541 TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23);
542 fMultV0->Fit(fpolc3, "RN");
543 fV0Cpol3 = fpolc3->GetParameter(0);
544
545 TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31);
546 fMultV0->Fit(fpolc4, "RN");
547 fV0Cpol4 = fpolc4->GetParameter(0);
548
549 TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
550 fMultV0->Fit(fpola1, "RN");
551 fV0Apol1 = fpola1->GetParameter(0);
552
553 TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
554 fMultV0->Fit(fpola2, "RN");
555 fV0Apol2 = fpola2->GetParameter(0);
556
557 TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
558 fMultV0->Fit(fpola3, "RN");
559 fV0Apol3 = fpola3->GetParameter(0);
560
561 TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
562 fMultV0->Fit(fpola4, "RN");
563 fV0Apol4 = fpola4->GetParameter(0);
564
565 for(Int_t i=0; i < 10; i++){
566
567 char nameQxa2[100];
568 snprintf(nameQxa2,100, "hQxa2m_%i", i);
569
570 char nameQya2[100];
571 snprintf(nameQya2,100, "hQya2m_%i", i);
572
573 char nameQxc2[100];
574 snprintf(nameQxc2,100, "hQxc2m_%i", i);
575
576 char nameQyc2[100];
577 snprintf(nameQyc2,100, "hQyc2m_%i", i);
578
579 AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
580 if(!contQxa2){
581 printf("OADB object %s is not available in the file\n", nameQxa2);
582 return 0;
583 }
584
585 if(!(contQxa2->GetObject(run))){
586 printf("OADB object %s is not available for run %i\n", nameQxa2, run);
587 return 0;
588 }
589
590 fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
591
592
593 AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
594 if(!contQya2){
595 printf("OADB object %s is not available in the file\n", nameQya2);
596 return 0;
597 }
598
599 if(!(contQya2->GetObject(run))){
600 printf("OADB object %s is not available for run %i\n", nameQya2, run);
601 return 0;
602 }
603
604 fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
605
606
607 AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
608 if(!contQxc2){
609 printf("OADB object %s is not available in the file\n", nameQxc2);
610 return 0;
611 }
612
613 if(!(contQxc2->GetObject(run))){
614 printf("OADB object %s is not available for run %i\n", nameQxc2, run);
615 return 0;
616 }
617
618 fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
619
620
621 AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
622 if(!contQyc2){
623 printf("OADB object %s is not available in the file\n", nameQyc2);
624 return 0;
625 }
626
627 if(!(contQyc2->GetObject(run))){
628 printf("OADB object %s is not available for run %i\n", nameQyc2, run);
629 return 0;
630 }
631
632 fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
633
634 }
635 return 1;
972a21ad 636}
637
638//______________________________________________________
639
640
c88234ad 641Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
642{
a72cae45 643 // Merge a list of AliSpectraAODEventCuts objects with this.
644 // Returns the number of merged objects (including this).
645
646 AliInfo("Merging");
647
648 if (!list)
649 return 0;
650
651 if (list->IsEmpty())
652 return 1;
653
654 TIterator* iter = list->MakeIterator();
655 TObject* obj;
656
657 // collections of all histograms
658 TList collections;
659
660 Int_t count = 0;
661
662 while ((obj = iter->Next())) {
663 AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
664 if (entry == 0)
665 continue;
666
667 TList * l = entry->GetOutputList();
668 collections.Add(l);
669 count++;
670 }
671
672 fOutput->Merge(&collections);
673
674 delete iter;
675
676 return count+1;
c88234ad 677}
678
7b364a00 679//______________________________________________________
680Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
a72cae45 681
682 //check Qvec and Centrality consistency
683 if(fCent>90.) return -999.;
684 if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
685 if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
686 if(v0side==2/*TPC*/){ if(fqTPC== -999.) return -999.; }
687
688 fQvecIntegral = 0x0;
689
690 if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
691 if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
692 if(v0side==2/*TPC*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("TPC"); }
693
694
695 Int_t ic = -999;
696
697 if(fQvecCalibType==1){
698 if(fNch<0.) return -999.;
699 ic = GetNchBin(fQvecIntegral);
700 } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
701
702 TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
703
704 TSpline *spline = 0x0;
705
706 if(v0side==0/*V0A*/){
707 if( CheckSplineArray(fSplineArrayV0A, ic) ) {
708 spline = (TSpline*)fSplineArrayV0A->At(ic);
709 //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
710 } else {
711 spline = new TSpline3(h1D,"sp3");
712 fSplineArrayV0A->AddAtAndExpand(spline,ic);
713 //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
714 }
715 }
716 else if(v0side==1/*V0C*/){
717 if( CheckSplineArray(fSplineArrayV0C, ic) ) {
718 spline = (TSpline*)fSplineArrayV0C->At(ic);
719 } else {
720 spline = new TSpline3(h1D,"sp3");
721 fSplineArrayV0C->AddAtAndExpand(spline,ic);
722 }
723 }
724 else if(v0side==2/*TPC*/){
725 if( CheckSplineArray(fSplineArrayTPC, ic) ) {
726 spline = (TSpline*)fSplineArrayTPC->At(ic);
727 } else {
728 spline = new TSpline3(h1D,"sp3");
729 fSplineArrayTPC->AddAtAndExpand(spline,ic);
730 }
731 }
732
733 Double_t percentile = -999.;
734 if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
735 if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
1d1bb564 736 if(v0side==2/*TPC*/){ percentile = 100*spline->Eval(fqTPC); }
a72cae45 737
738 if(percentile>100. || percentile<0.) return -999.;
739
740 return percentile;
7b364a00 741}
742
ed33b9ef 743//______________________________________________________
725da720 744Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side, Int_t type=1){
a72cae45 745
746 if(!fIsMC) return -999.;
747
748 // V0A efficiecy
749 if(type==1){
750 fV0Aeff = 0x0;
751 fV0Aeff = (TH1F*)fQvecIntList->FindObject("vzeroa_efficiency_prim_plus_sec_over_gen");
752 if(!fV0Aeff) return -999.;
753 }
754
755 // get MC array
756 TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
757
758 if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
759
760 Double_t Qx2mc = 0., Qy2mc = 0.;
761 Double_t mult2mc = 0;
762
763 Int_t nMC = arrayMC->GetEntries();
764
765 if(type==0){ // type==0: q-vec from tracks in vzero acceptance
766
767 // loop on generated
768 for (Int_t iMC = 0; iMC < nMC; iMC++){
769 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
770 if(!partMC->Charge()) continue;//Skip neutrals
771
772 // check vzero side
773 if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
774
775 // Calculate Qvec components
776 Qx2mc += TMath::Cos(2.*partMC->Phi());
777 Qy2mc += TMath::Sin(2.*partMC->Phi());
778 mult2mc++;
779
780 }// end loop on generated.
781
782 }//end if on type==0
783
784 else if(type==1){ // type==1 (default): q-vec from vzero
785
786 // only used in qgen_vzero
787 Double_t multv0mc[64];
788 for(Int_t iCh=0;iCh<64;iCh++) multv0mc[iCh]=0; // initialize multv0mc to zero
789
790 // loop on generated
791 for (Int_t iMC = 0; iMC < nMC; iMC++){
792
793 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
794 if(!partMC->Charge()) continue;//Skip neutrals
795
796 // check vzero side
797 if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
798
799 //get v0 channel of gen track
800 Int_t iv0 = CheckVZEROchannel(v0side, partMC->Eta(), partMC->Phi());
801
802 //use the efficiecy as a weigth
803 multv0mc[iv0] += fV0Aeff->GetFunction("f")->Eval(partMC->Pt());
804
805 //calculate multiplicity for each vzero channell
806 //multv0mc[iv0]++;
807
808 }// end loop on generated.
809
810 Int_t firstCh=-1,lastCh=-1;
811 if (v0side == 0) {
812 firstCh = 32;
813 lastCh = 64;
814 }
815 else if (v0side == 1) {
816 firstCh = 0;
817 lastCh = 32;
818 }
819 for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
820 Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
821 Double_t mult = multv0mc[iCh];
822 Qx2mc += mult*TMath::Cos(2*phi);
823 Qy2mc += mult*TMath::Sin(2*phi);
824 mult2mc += mult;
825
826 ((TH2F*)fOutput->FindObject("fV0Mmc"))->Fill(iCh,multv0mc[iCh]);
827 }
828
829 }//end if on type==1
830
831 // return q vector
832 fQvecMC = TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
833
834 return fQvecMC;
835
ed33b9ef 836}
837
cafe9ad1 838//______________________________________________________
725da720 839Int_t AliSpectraAODEventCuts::CheckVZEROchannel(Int_t vzeroside, Double_t eta, Double_t phi){
a72cae45 840
841 //VZEROA eta acceptance
842 Int_t ring = -1.;
843
844 if ( vzeroside == 0){
845 if (eta > 4.5 && eta <= 5.1 ) ring = 0;
846 if (eta > 3.9 && eta <= 4.5 ) ring = 1;
847 if (eta > 3.4 && eta <= 3.9 ) ring = 2;
848 if (eta > 2.8 && eta <= 3.4 ) ring = 3;
849 } else if (vzeroside == 1){
850 if (eta > -3.7 && eta <= -3.2 ) ring = 0;
851 if (eta > -3.2 && eta <= -2.7 ) ring = 1;
852 if (eta > -2.7 && eta <= -2.2 ) ring = 2;
853 if (eta > -2.2 && eta <= -1.7 ) ring = 3;
854 }
855
856
857 //VZEROAC phi acceptance
858 Int_t phisector= -1;
859
860
861 if ( phi > 0. && phi <= TMath::Pi()/4. ) phisector = 0;
862
863 else if (phi > TMath::Pi()/4. && phi <= TMath::Pi()/2.) phisector = 1;
864
865 else if (phi > TMath::Pi()/2 && phi <= (3./4.)*TMath::Pi() ) phisector =2;
866
867 else if (phi > (3./4.)*TMath::Pi() && phi <= TMath::Pi() ) phisector = 3;
868
869 else if (phi > TMath::Pi() && phi <= (5./4.)*TMath::Pi() ) phisector = 4;
870
871 else if (phi > (5./4.)*TMath::Pi() && phi <= (3./2.)*TMath::Pi() ) phisector = 5;
872
873 else if (phi > (3./2.)*TMath::Pi() && phi <= (7./4.)*TMath::Pi() ) phisector = 6;
874
875 else if (phi > (7./4.)*TMath::Pi() && phi <= TMath::TwoPi() ) phisector = 7;
876
877 if (vzeroside ==0) return Int_t(32+(phisector+(ring*8.))); // iv0 >= 32 -> V0A
878 if (vzeroside ==1) return Int_t(phisector+(ring*8.)); // iv0 < 32 -> V0C
879
880 else return -999.;
725da720 881}
882
883//______________________________________________________
884Int_t AliSpectraAODEventCuts::CheckVZEROacceptance(Double_t eta){
a72cae45 885
886 // eval VZERO side - FIXME Add TPC!
887
888 if(eta > 2.8 && eta < 5.1) return 0; //VZEROA
889
890 else if (eta > -3.7 && eta < -1.7) return 1; //VZEROC
891
892 return -999.;
893
725da720 894}
895
896//______________________________________________________
897Double_t AliSpectraAODEventCuts::GetQvecPercentileMC(Int_t v0side, Int_t type=1){
a72cae45 898
899 //check Qvec and Centrality consistency
900 if(fCent>90.) return -999.;
901
902 Double_t qvec = CalculateQVectorMC(v0side, type);
903 if(qvec==-999.) return -999.;
904
905 fQgenIntegral = 0x0;
906
907 if(type==0){
908 if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen"); }
909 if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen"); }
910 } else if (type==1){
911 if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen_vzero"); }
912 if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen_vzero"); }
913 }
914
915
916 //FIXME you need a check on the TH2D, add AliFatal or a return.
917
918 Int_t ic = -999;
919
920 if(fQvecCalibType==1){
921 if(fNch<0.) return -999.;
922 ic = GetNchBin(fQgenIntegral);
923 } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
924
925 TH1D *h1D = (TH1D*)fQgenIntegral->ProjectionY("h1Dgen",ic+1,ic+1);
926
927 TSpline *spline = 0x0;
928
929 if(v0side==0/*V0A*/){
930 if( CheckSplineArray(fSplineArrayV0Agen, ic) ) {
931 spline = (TSpline*)fSplineArrayV0Agen->At(ic);
932 //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
933 } else {
934 spline = new TSpline3(h1D,"sp3");
935 fSplineArrayV0Agen->AddAtAndExpand(spline,ic);
936 }
937 }
938 else if(v0side==1/*V0C*/){
939 if( CheckSplineArray(fSplineArrayV0Cgen, ic) ) {
940 spline = (TSpline*)fSplineArrayV0Cgen->At(ic);
941 } else {
942 spline = new TSpline3(h1D,"sp3");
943 fSplineArrayV0Cgen->AddAtAndExpand(spline,ic);
944 }
945 }
946
947 Double_t percentile = 100*spline->Eval(qvec);
948
949 if(percentile>100. || percentile<0.) return -999.;
950
951 return percentile;
cafe9ad1 952
953}
954
7b364a00 955//______________________________________________________
91083d79 956Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr, Int_t n){
a72cae45 957 //check TSpline array consistency
958
959 // Int_t n = (Int_t)fCent; //FIXME should be ok for icentr and nch
960
961 if(splarr->GetSize()<n) return kFALSE;
962
963 if(!splarr->At(n)) return kFALSE;
964
965 return kTRUE;
966
7b364a00 967}
91083d79 968
969//______________________________________________________
34f4a5fc 970Int_t AliSpectraAODEventCuts::GetNchBin(TH2D * h){
a72cae45 971
972 Double_t xmax = h->GetXaxis()->GetXmax();
973
974 if(fNch>xmax) return (Int_t)h->GetNbinsX();
975
976 return (fNch*h->GetNbinsX())/h->GetXaxis()->GetXmax();
977
91083d79 978}
979