]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothEventCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding the multiplcity lodaer class to TestAOD directory and removing second copy...
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraBothEventCuts.cxx
CommitLineData
239a080a 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// AliSpectraBothEventCuts class
19//-----------------------------------------------------------------
20
21#include "TChain.h"
22#include "TTree.h"
23#include "TLegend.h"
24#include "TH1F.h"
25#include "TH1I.h"
26#include "TH2F.h"
27#include "TCanvas.h"
28#include "AliAnalysisTask.h"
29#include "AliAnalysisManager.h"
30#include "AliAODTrack.h"
31#include "AliAODMCParticle.h"
32#include "AliAODEvent.h"
33#include "AliESDEvent.h"
34#include "AliAODInputHandler.h"
35#include "AliAnalysisTaskESDfilter.h"
36#include "AliAnalysisDataContainer.h"
37#include "AliSpectraBothEventCuts.h"
38#include "AliSpectraBothTrackCuts.h"
49ed9c10 39#include "AliAnalysisUtils.h"
40#include "AliPPVsMultUtils.h"
21e4a3cd 41#include "AliGenPythiaEventHeader.h"
239a080a 42//#include "AliSpectraBothHistoManager.h"
43#include <iostream>
c239ee06 44#include "TParameter.h"
239a080a 45
46using namespace std;
47
48ClassImp(AliSpectraBothEventCuts)
49
b3ea73e1 50AliSpectraBothEventCuts::AliSpectraBothEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0),fAODEvent(AliSpectraBothTrackCuts::kAODobject), fTrackBits(0),fIsMC(0),fCentEstimator(""), fUseCentPatchAOD049(0), fUseSDDPatchforLHC11a(kDoNotCheckforSDD),fTriggerSettings(AliVEvent::kMB),fTrackCuts(0),
239a080a 51fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(0),fMaxChi2perNDFforVertex(0),
49ed9c10 52fMinRun(0),fMaxRun(0),fetarangeofmultiplicitycut(0.0),fUseAliPPVsMultUtils(false),
21e4a3cd 53fNMCProcessType(0),fEventMCProcessType(0),fEventMCProcessTypeIncluded(0),
239a080a 54fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
0ab8c127 55,fHistoEP(0),fHistoVtxAftSelwithoutZvertexCut(0),fHistoVtxalltriggerEventswithMCz(0),fHistoVtxAftSelwithoutZvertexCutusingMCz(0),fHistoRunNumbers(0),
49ed9c10 56fHistoCentrality(0),fHistoMultiplicty(0),fAnalysisUtils(0),fAliPPVsMultUtils(0)
0ab8c127 57
239a080a 58{
8fda510f 59 // Constructori
2d98dd91 60 // Bool_t oldStatus = TH1::AddDirectoryStatus();
61 //TH1::AddDirectory(kFALSE);
62 // fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
63 // fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",300,-15,15);
64// fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",300,-15,15);
65// fHistoVtxAftSelwithoutZvertexCut=new TH1F("fHistoVtxAftSelwithoutZvertexcut", "Vtx distr after event selection without Z vertex cut",300,-15,15);
66// fHistoVtxalltriggerEventswithMCz=new TH1F("fHistoVtxalltriggerEventswithMCz", "generated z vertex position",300,-15,15);
67// fHistoVtxAftSelwithoutZvertexCutusingMCz=new TH1F("fHistoVtxAftSelwithoutZvertexCutusingMCz", "Vtx distr after event selection without Z vertex cut using MC z",300,-15,15);
68// fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
69// fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
70 // fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
239a080a 71 //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
72 //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
2d98dd91 73 // fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
74 // fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
239a080a 75 fCentralityCutMin = 0.0; // default value of centrality cut minimum, 0 ~ no cut
76 fCentralityCutMax = 10000.0; // default value of centrality cut maximum, ~ no cut
77 // fQVectorPosCutMin=0.0;
78 // fQVectorPosCutMax=10000.0;
79 // fQVectorNegCutMin=0.0;
80 // fQVectorNegCutMax=10000.0;
81 fQVectorCutMin=0.0;
82 fQVectorCutMax=10000.0;
83 fVertexCutMin=-10.0;
84 fVertexCutMax=10.0;
85 fMultiplicityCutMin=-1.0;
86 fMultiplicityCutMax=-1.0;
87 fTrackBits=1;
b3ea73e1 88 fCentEstimator="V0M";
21e4a3cd 89 fMaxChi2perNDFforVertex=-1;
2d98dd91 90 // TH1::AddDirectory(oldStatus);
8fda510f 91}
92//______________________________________________________
93
94AliSpectraBothEventCuts::~AliSpectraBothEventCuts()
95{
96 if(fHistoCuts)
97 delete fHistoCuts;
98 if(fHistoVtxBefSel)
99 delete fHistoVtxBefSel;
100 if(fHistoVtxAftSel)
101 delete fHistoVtxAftSel;
102 if(fHistoVtxAftSelwithoutZvertexCut)
103 delete fHistoVtxAftSelwithoutZvertexCut;
104 if(fHistoVtxalltriggerEventswithMCz)
105 delete fHistoVtxalltriggerEventswithMCz;
106 if(fHistoVtxAftSelwithoutZvertexCutusingMCz)
107 delete fHistoVtxAftSelwithoutZvertexCutusingMCz;
108 if(fHistoEtaBefSel)
109 delete fHistoEtaBefSel;
110 if(fHistoEtaAftSel)
111 delete fHistoEtaAftSel ;
112 if(fHistoNChAftSel)
113 delete fHistoNChAftSel;
114 if(fHistoQVector)
115 delete fHistoQVector;
116 if(fHistoEP)
117 delete fHistoEP;
2d98dd91 118 if(fHistoRunNumbers)
119 delete fHistoRunNumbers;
0ab8c127 120 if(fHistoCentrality)
121 delete fHistoCentrality;
122 if(fHistoMultiplicty)
123 delete fHistoMultiplicty;
4954491e 124 if(fAnalysisUtils)
125 delete fAnalysisUtils;
49ed9c10 126 if(fAliPPVsMultUtils)
127 delete fAliPPVsMultUtils;
21e4a3cd 128 if(fEventMCProcessType)
129 delete [] fEventMCProcessType;
0ab8c127 130
239a080a 131}
2d98dd91 132//______________________________________________________
133void AliSpectraBothEventCuts::InitHisto()
134{
135 Bool_t oldStatus = TH1::AddDirectoryStatus();
136 TH1::AddDirectory(kFALSE);
137 if(!fHistoCuts)
138 fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
139 if(!fHistoVtxBefSel )
140 fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",300,-15,15);
141 if(!fHistoVtxAftSel)
142 fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",300,-15,15);
143 if(!fHistoVtxAftSelwithoutZvertexCut)
144 fHistoVtxAftSelwithoutZvertexCut=new TH1F("fHistoVtxAftSelwithoutZvertexcut", "Vtx distr after event selection without Z vertex cut",300,-15,15);
145 if(!fHistoVtxalltriggerEventswithMCz)
146 fHistoVtxalltriggerEventswithMCz=new TH1F("fHistoVtxalltriggerEventswithMCz", "generated z vertex position",300,-15,15);
147 if(!fHistoVtxAftSelwithoutZvertexCutusingMCz)
148 fHistoVtxAftSelwithoutZvertexCutusingMCz=new TH1F("fHistoVtxAftSelwithoutZvertexCutusingMCz", "Vtx distr after event selection without Z vertex cut using MC z",300,-15,15);
149 if(!fHistoEtaBefSel)
150 fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
151 if(!fHistoEtaAftSel)
152 fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
153 if(!fHistoNChAftSel)
154 fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
155 //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
156 //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
157 if(!fHistoQVector)
158 fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
159 if(!fHistoEP)
160 fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
161 if(!fHistoRunNumbers)
162 {
163 if(fMaxRun>fMinRun&&fMinRun>=0)
164 fHistoRunNumbers=new TH1F("fHistoRunNumbers","Run numbers",fMaxRun-fMinRun+1,fMinRun-0.5,fMaxRun+0.5);
165 else
166 fHistoRunNumbers=new TH1F("fHistoRunNumbers","Run numbers",1001,120000-.5,121000+0.5);
239a080a 167
2d98dd91 168 }
0ab8c127 169 if(!fHistoCentrality)
170 fHistoCentrality = new TH2F("fHistoCentrality", "centrality",2,0,2,100,0.0,100);
171
172 if(!fHistoMultiplicty)
49ed9c10 173 fHistoMultiplicty= new TH2F("fHistoMultiplicty", "multiplicty estimator",2,0,2,155,-4.5,150.5);
0ab8c127 174
2d98dd91 175 TH1::AddDirectory(oldStatus);
176}
239a080a 177//______________________________________________________
34347952 178Bool_t AliSpectraBothEventCuts::IsSelected(AliVEvent * aod,AliSpectraBothTrackCuts* trackcuts,Bool_t isMC,Double_t mcZ,TH1F* managerhisteventcuts)
239a080a 179{
180 // Returns true if Event Cuts are selected and applied
181 fAOD = aod;
182 fTrackCuts = trackcuts;
183 fHistoCuts->Fill(kProcessedEvents);
34347952 184 if(managerhisteventcuts)
185 managerhisteventcuts->Fill(0);
2d98dd91 186 fHistoRunNumbers->Fill(aod->GetRunNumber());
1c23d41e 187 Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerSettings);//FIXME we can add the trigger mask here
d0592deb 188 if(!IsPhysSel)
189 return IsPhysSel;
190 if(fAnalysisUtils) // we check for pile-up
4954491e 191 IsPhysSel = (!fAnalysisUtils->IsPileUpEvent(fAOD));
192 if(!IsPhysSel)
193 return IsPhysSel;
194
8bb435b0 195 if(isMC)
196 fHistoVtxalltriggerEventswithMCz->Fill(mcZ);
239a080a 197 //loop on tracks, before event selection, filling QA histos
4954491e 198 AliESDEvent* esdevent=0x0;
239a080a 199 AliAODEvent* aodevent=0x0;
200 Bool_t isSDD=kFALSE;
4954491e 201 TString nameoftrack(fAOD->ClassName());
239a080a 202 if(!nameoftrack.CompareTo("AliESDEvent"))
203 {
204 fAODEvent=AliSpectraBothTrackCuts::kESDobject;
205
206 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
207 {
208 esdevent=dynamic_cast<AliESDEvent*>(fAOD);
737ad8ac 209 if(!esdevent)
94a8ed71 210 return kFALSE;
239a080a 211 if(esdevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
212 isSDD=kTRUE;
213 }
214 }
215 else if(!nameoftrack.CompareTo("AliAODEvent"))
216 {
217 fAODEvent=AliSpectraBothTrackCuts::kAODobject;
218 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
219 {
220 aodevent=dynamic_cast<AliAODEvent*>(fAOD);
737ad8ac 221 if(!aodevent)
94a8ed71 222 return kFALSE;
239a080a 223 if(aodevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
224 isSDD=kTRUE;
225 }
226 }
227 else
228 return false;
229 if(fUseSDDPatchforLHC11a==kwithSDD&&isSDD==kFALSE)
230 return false;
231 if(fUseSDDPatchforLHC11a==kwithoutSDD&&isSDD==kTRUE)
232 return false;
233
239a080a 234 fHistoCuts->Fill(kPhysSelEvents);
34347952 235 if(managerhisteventcuts)
236 managerhisteventcuts->Fill(1);
239a080a 237
238
239
240
bca0c28e 241 const AliVVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
242
4954491e 243 if(vertex)
244 fHistoVtxBefSel->Fill(vertex->GetZ());
239a080a 245 fIsSelected =kFALSE;
49ed9c10 246 if(CheckVtx())
247 { //selection on vertex
547e52a0 248
4954491e 249 fIsSelected=kTRUE;
239a080a 250 }
8bb435b0 251 if(fIsSelected&&vertex)
49ed9c10 252 {
253 fHistoVtxAftSelwithoutZvertexCut->Fill(vertex->GetZ());
254 if(isMC)
255 fHistoVtxAftSelwithoutZvertexCutusingMCz->Fill(mcZ);
256 if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
257 {
547e52a0 258 fIsSelected=kTRUE;
259 fHistoVtxAftSel->Fill(vertex->GetZ());
49ed9c10 260 }
261 else
262 {
547e52a0 263 fIsSelected=kFALSE;
49ed9c10 264 }
239a080a 265 }
4954491e 266 if(fIsSelected)
49ed9c10 267 {
268 if( CheckCentralityCut() && CheckMultiplicityCut() && CheckQVectorCut())
269 fIsSelected=kTRUE;
77898649 270 else
271 fIsSelected=kFALSE;
49ed9c10 272 }
4954491e 273
274
275
239a080a 276 Int_t Nch=0;
49ed9c10 277 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
278 {
279 AliVTrack* track =dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
239a080a 280 /* if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
281 track=dynamic_cast<AliVTrack*>(esdevent->GetTrack(iTracks));
282 else if (fAODEvent==AliSpectraBothTrackCuts::kAODobject)
283 track=dynamic_cast<AliVTrack*>(aodevent->GetTrack(iTracks));
284 else return false;*/
285
49ed9c10 286 if (!fTrackCuts->IsSelected(track,kFALSE))
287 continue;
288 fHistoEtaBefSel->Fill(track->Eta());
289 if(fIsSelected)
290 {
291 fHistoEtaAftSel->Fill(track->Eta());
292 Nch++;
293 }
239a080a 294 }
4954491e 295 if(fIsSelected)
77898649 296 {
297 fHistoNChAftSel->Fill(Nch);
298 fHistoCuts->Fill(kAcceptedEvents);
34347952 299 if(managerhisteventcuts)
300 managerhisteventcuts->Fill(2);
301
77898649 302 }
239a080a 303 return fIsSelected;
304}
305
306//______________________________________________________
547e52a0 307Bool_t AliSpectraBothEventCuts::CheckVtx()
239a080a 308{
309 // reject events outside of range
310 const AliVVertex * vertex = fAOD->GetPrimaryVertex();
311 //when moving to 2011 wìone has to add a cut using SPD vertex.
312 //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
313 if (!vertex)
314 {
315 fHistoCuts->Fill(kVtxNoEvent);
316 return kFALSE;
317 }
bca0c28e 318 if(vertex->GetNContributors()<1)
239a080a 319 {
bca0c28e 320
321 fHistoCuts->Fill(kZeroCont);
239a080a 322 return kFALSE;
323
324 }
bca0c28e 325
326 TString tmp(vertex->GetTitle());
327 if(tmp.Contains("NoConstraint"))
239a080a 328 {
bca0c28e 329 fHistoCuts->Fill(kTPCasPV);
330 return kFALSE;
331 }
332
333
547e52a0 334 // if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
335 // {
336 // return kTRUE;
337 // }
4954491e 338 if(!CheckVtxChi2perNDF())
339 return kFALSE;
340
341 fHistoCuts->Fill(kGoodVtx);
547e52a0 342 //return kFALSE;
343 return kTRUE;
239a080a 344}
345
346//______________________________________________________
347Bool_t AliSpectraBothEventCuts::CheckCentralityCut()
348{
349 // Check centrality cut
49ed9c10 350 if ( fCentralityCutMax<0.0 && fCentralityCutMin<0.0 )
351 return kTRUE;
c239ee06 352 Double_t cent=0;
353 Bool_t validcent=kFALSE;
354 if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
49ed9c10 355 {
c239ee06 356
357 TString par1name("cent");
358 par1name+=fCentEstimator;
359 cout<<par1name.Data()<<endl;
360 TParameter<Double_t>* par= dynamic_cast<TParameter<Double_t>*>(fAOD->FindListObject(par1name.Data()));
361 if(par)
362 {
363 validcent=kTRUE;
364 cent=par->GetVal();
365 }
366
367 }
368 if(!validcent)
369 {
370 if(fUseAliPPVsMultUtils)
371 {
372 if(!fAliPPVsMultUtils)
373 fAliPPVsMultUtils=new AliPPVsMultUtils();
374 cent=fAliPPVsMultUtils->GetMultiplicityPercentile(fAOD,fCentEstimator.Data());
375 }
376 else
377 {
378 if(!fUseCentPatchAOD049)
379 cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentEstimator.Data());
380 else
381 cent=ApplyCentralityPatchAOD049();
382 }
49ed9c10 383 }
0ab8c127 384 fHistoCentrality->Fill(0.5,cent);
77898649 385 if ( (cent < fCentralityCutMax) && (cent >= fCentralityCutMin) )
0ab8c127 386 {
387 fHistoCentrality->Fill(1.5,cent);
388 return kTRUE;
389 }
239a080a 390 fHistoCuts->Fill(kVtxCentral);
391
392 return kFALSE;
393}
394
395//______________________________________________________
396Bool_t AliSpectraBothEventCuts::CheckMultiplicityCut()
397{
398 // Check multiplicity cut
0ab8c127 399 if(fMultiplicityCutMin<0 && fMultiplicityCutMax<0)
400 return kTRUE;
401 Int_t Ncharged=-1;
402 if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
c239ee06 403 {
404
405 TParameter<Int_t>* par=0;
406 if(TMath::Abs(0.8-fetarangeofmultiplicitycut)<0.1)
407 par= dynamic_cast<TParameter<Int_t>*>(fAOD->FindListObject("Ncheta0dot8"));
408 else
409 par= dynamic_cast<TParameter<Int_t>*>(fAOD->FindListObject("Ncheta0dot5"));
410 if(par)
411 {
412 Ncharged=par->GetVal();
413 }
414 else
415 {
416 AliESDEvent* esdevent=dynamic_cast<AliESDEvent*>(fAOD);
417 AliESDtrackCuts::MultEstTrackType estType = esdevent->GetPrimaryVertexTracks()->GetStatus() ? AliESDtrackCuts::kTrackletsITSTPC : AliESDtrackCuts::kTracklets;
418 Ncharged=AliESDtrackCuts::GetReferenceMultiplicity(esdevent,estType,fetarangeofmultiplicitycut);
419 }
0ab8c127 420 }
421 else if(fAODEvent==AliSpectraBothTrackCuts::kAODobject)
422 {
423 AliAODEvent* aodevent=0x0;
424 aodevent=dynamic_cast<AliAODEvent*>(fAOD);
0a918d8d 425 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
426 if(!header) AliFatal("Not a standard AOD");
427
0ab8c127 428 if(TMath::Abs(0.8-fetarangeofmultiplicitycut)<0.1)
0a918d8d 429 Ncharged=header->GetRefMultiplicityComb08();
0ab8c127 430 else if (TMath::Abs(0.5-fetarangeofmultiplicitycut)<0.1)
0a918d8d 431 Ncharged=header->GetRefMultiplicityComb05();
0ab8c127 432 else
433 Ncharged=-1;
434 }
435 else
436 return kFALSE;
239a080a 437
0ab8c127 438 fHistoMultiplicty->Fill(0.5,Ncharged);
21e4a3cd 439 if(Ncharged>=fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)
0ab8c127 440 {
441 fHistoMultiplicty->Fill(1.5,Ncharged);
442 return kTRUE;
443 }
444 fHistoCuts->Fill(kVtxCentral);
239a080a 445 return kFALSE;
446}
447
448//______________________________________________________
449Bool_t AliSpectraBothEventCuts::CheckQVectorCut()
450{
451 if(fQVectorCutMin<0.0 && fQVectorCutMax<0.0)return kTRUE;
452 // Check qvector cut
453 /// FIXME: Q vector
454 // //Selection on QVector, before ANY other selection on the event
455 // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
456 // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
457 // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
458
459 // Int_t multPos = 0;
460 // Int_t multNeg = 0;
461 // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
462 // AliAODTrack* aodTrack = fAOD->GetTrack(iT);
463 // if (!aodTrack->TestFilterBit(fTrackBits)) continue;
464 // if (aodTrack->Eta() >= 0){
465 // multPos++;
466 // Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
467 // Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
468 // } else {
469 // multNeg++;
470 // Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
471 // Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
472 // }
473 // }
474 // Double_t qPos=-999;
475 // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
476 // Double_t qNeg=-999;
477 // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
478 //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
479
480 Double_t qxEPVZERO = 0, qyEPVZERO = 0;
481 Double_t qVZERO = -999;
482 Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
483
484 qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
485 if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
486 fHistoQVector->Fill(qVZERO);
487 fHistoEP->Fill(psi);
488
489 fHistoCuts->Fill(kQVector);
490 // fHistoQVectorPos->Fill(qPos);
491 // fHistoQVectorNeg->Fill(qNeg);
492 return kTRUE;
493}
494//____________________________________________________________
495 Bool_t AliSpectraBothEventCuts::CheckVtxChi2perNDF()
496 {
497 if(fMaxChi2perNDFforVertex<0)
498 return kTRUE;
499 const AliVVertex * vertex = fAOD->GetPrimaryVertex();
500 if(TMath::Abs(vertex->GetChi2perNDF())>fMaxChi2perNDFforVertex)
501 return kFALSE;
502 return kTRUE;
503 }
504
505
506
507//______________________________________________________
508void AliSpectraBothEventCuts::PrintCuts()
509{
510 // print info about event cuts
511 cout << "Event Stats" << endl;
2d98dd91 512 if(fHistoCuts)
513 {
514 cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
515 cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
516 cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
4954491e 517 cout << " > With good veretx: " << fHistoCuts->GetBinContent(kGoodVtx + 1) << endl;
2d98dd91 518 cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
519 cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
520 cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
521 }
239a080a 522 cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
523 // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
524 // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
525 cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
526 cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
527 cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
528 cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
529}
530//______________________________________________________
531
532Double_t AliSpectraBothEventCuts::ApplyCentralityPatchAOD049()
533{
534 //
535 //Apply centrality patch for AOD049 outliers
536 //
537 // if (fCentralityType!="V0M") {
538 // AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
539 // return -999.0;
540 // }
541 AliCentrality *centrality = fAOD->GetCentrality();
542 if (!centrality) {
543 AliWarning("Cannot get centrality from AOD event.");
544 return -999.0;
545 }
546
547 Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
548 /*
549 Bool_t isSelRun = kFALSE;
550 Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
551 if(cent<0){
552 Int_t quality = centrality->GetQuality();
553 if(quality<=1){
554 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
555 } else {
556 Int_t runnum=aodEvent->GetRunNumber();
557 for(Int_t ir=0;ir<5;ir++){
558 if(runnum==selRun[ir]){
559 isSelRun=kTRUE;
560 break;
561 }
562 }
563 if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
564 }
565 }
566 */
567 if(cent>=0.0) {
568 Float_t v0 = 0.0;
569 AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
570 AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
571 v0+=aodV0->GetMTotV0A();
572 v0+=aodV0->GetMTotV0C();
573 if ( (cent==0) && (v0<19500) ) {
574 AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
575 return -999.0;
576 }
577 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
578 Float_t val = 1.30552 + 0.147931 * v0;
579
580 Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
581 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
582 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
583 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
584 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
585 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
586 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
587 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
588 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
589 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
590 };
591
592 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
593 AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
594 return -999.0;
595 }
596 } else {
597 //force it to be -999. whatever the negative value was
598 cent = -999.;
599 }
600 return cent;
601}
602
603//______________________________________________________
604
605
606Long64_t AliSpectraBothEventCuts::Merge(TCollection* list)
607{
608 // Merge a list of AliSpectraBothEventCuts objects with this.
609 // Returns the number of merged objects (including this).
610
611 // AliInfo("Merging");
612
613 if (!list)
614 return 0;
615
616 if (list->IsEmpty())
617 return 1;
618
619 TIterator* iter = list->MakeIterator();
620 TObject* obj;
621
622 // collections of all histograms
623 TList collections;//FIXME we should use only 1 collection
624 TList collections_histoVtxBefSel;
625 TList collections_histoVtxAftSel;
626 TList collections_histoEtaBefSel;
627 TList collections_histoEtaAftSel;
628 TList collections_histoNChAftSel;
629 // TList collections_histoQVectorPos;
630 // TList collections_histoQVectorNeg;
631 TList collections_histoQVector;
632 TList collections_histoEP;
3a71f081 633 TList collections_histoVtxAftSelwithoutZvertexCut;
634 TList collections_histoVtxalltriggerEventswithMCz;
635 TList collections_histoVtxAftSelwithoutZvertexCutusingMCz;
2d98dd91 636 TList collections_histoRunNumbers;
0ab8c127 637 TList collections_histoCentrality;
638 TList collections_histoMultiplicty;
639
239a080a 640 Int_t count = 0;
641
642 while ((obj = iter->Next())) {
643 AliSpectraBothEventCuts* entry = dynamic_cast<AliSpectraBothEventCuts*> (obj);
644 if (entry == 0)
645 continue;
646
82c27a75 647 TH1I * histo = entry->GetHistoCuts();
648 if(histo)
649 collections.Add(histo);
239a080a 650 TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
82c27a75 651 if(histo_histoVtxBefSel)
652 collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
239a080a 653 TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
82c27a75 654 if(histo_histoVtxAftSel)
655 collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
656 TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
657 if(histo_histoEtaBefSel)
658 collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
659 TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
660 if(histo_histoEtaAftSel)
661 collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
662 TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
663 if(histo_histoNChAftSel)
664 collections_histoNChAftSel.Add(histo_histoNChAftSel);
239a080a 665 // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
666 // collections_histoQVectorPos.Add(histo_histoQVectorPos);
667 // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
82c27a75 668 // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
669 TH1F * histo_histoQVector = entry->GetHistoQVector();
670 if(histo_histoQVector)
671 collections_histoQVector.Add(histo_histoQVector);
239a080a 672 TH1F * histo_histoEP = entry->GetHistoEP();
82c27a75 673 if(histo_histoEP)
674 collections_histoEP.Add(histo_histoEP);
3a71f081 675 TH1F* histo_histoVtxAftSelwithoutZvertexCut=entry->GetHistoVtxAftSelwithoutZvertexCut();
82c27a75 676 if(histo_histoVtxAftSelwithoutZvertexCut)
677 collections_histoVtxAftSelwithoutZvertexCut.Add(histo_histoVtxAftSelwithoutZvertexCut);
3a71f081 678 TH1F* histo_histoVtxalltriggerEventswithMCz=entry->GetHistoVtxGenerated();
82c27a75 679 if(histo_histoVtxalltriggerEventswithMCz)
680 collections_histoVtxalltriggerEventswithMCz.Add(histo_histoVtxalltriggerEventswithMCz);
2d98dd91 681
682 TH1F* histo_histoVtxAftSelwithoutZvertexCutusingMCz=entry->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
82c27a75 683 if(histo_histoVtxAftSelwithoutZvertexCutusingMCz)
684 collections_histoVtxAftSelwithoutZvertexCutusingMCz.Add(histo_histoVtxAftSelwithoutZvertexCutusingMCz);
2d98dd91 685
686 TH1F* histo_histoRunNumbers=entry->GetHistoRunNumbers();
82c27a75 687 if(histo_histoRunNumbers)
688 collections_histoRunNumbers.Add(histo_histoRunNumbers);
0ab8c127 689
690 TH2F* histo_histoCentrality=entry->GetHistoCentrality();
82c27a75 691 if(histo_histoCentrality)
692 collections_histoCentrality.Add(histo_histoCentrality);
0ab8c127 693
694TH2F* histo_histoMultiplicty=entry->GetHistoMultiplicty();
82c27a75 695 if(histo_histoMultiplicty)
696 collections_histoMultiplicty.Add(histo_histoMultiplicty);
0ab8c127 697
698
239a080a 699 count++;
700 }
82c27a75 701 if(fHistoCuts)
702 fHistoCuts->Merge(&collections);
703 if(fHistoVtxBefSel)
704 fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
705 if(fHistoVtxAftSel)
706 fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
707 if(fHistoEtaBefSel)
708 fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
709 if(fHistoEtaAftSel)
710 fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
711 if(fHistoNChAftSel)
712 fHistoNChAftSel->Merge(&collections_histoNChAftSel);
239a080a 713 // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
714 // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
82c27a75 715 if(fHistoQVector)
716 fHistoQVector->Merge(&collections_histoQVector);
717 if(fHistoEP)
718 fHistoEP->Merge(&collections_histoEP);
719 if(fHistoVtxAftSelwithoutZvertexCut)
720 fHistoVtxAftSelwithoutZvertexCut->Merge(&collections_histoVtxAftSelwithoutZvertexCut);
721 if(fHistoVtxalltriggerEventswithMCz)
722 fHistoVtxalltriggerEventswithMCz->Merge(&collections_histoVtxalltriggerEventswithMCz);
723 if(fHistoVtxAftSelwithoutZvertexCutusingMCz)
724 fHistoVtxAftSelwithoutZvertexCutusingMCz->Merge(&collections_histoVtxAftSelwithoutZvertexCutusingMCz);
725 if(fHistoRunNumbers)
726 fHistoRunNumbers->Merge(&collections_histoRunNumbers);
727 if(fHistoCentrality)
728 fHistoCentrality->Merge(&collections_histoCentrality);
729 if(fHistoMultiplicty)
730 fHistoMultiplicty->Merge(&collections_histoMultiplicty);
0ab8c127 731
239a080a 732
733 delete iter;
734
735 return count+1;
736}
2d98dd91 737//__________________________________________________________________________________________________________
738void AliSpectraBothEventCuts::SetRunNumberRange(Int_t min, Int_t max)
739{
740 if(max>min&&min>=0)
741 {
742 fMinRun=min;
743 fMaxRun=max;
744 }
745}
21e4a3cd 746//__________________________________________________________________________________________________________
747Bool_t AliSpectraBothEventCuts::CheckMCProcessType(AliMCEvent* mcevent)
748{
749 if(fNMCProcessType<0)
750 return kTRUE;
751 if(!mcevent)
752 return kFALSE;
753 AliHeader* aHeader=mcevent->Header();
754 if(!aHeader)
755 return kFALSE;
756 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
757 if(!pythiaGenHeader)
758 return kFALSE;
759 Int_t processtype=pythiaGenHeader->ProcessType();
760 for(int i=0;i<fNMCProcessType;i++)
761 {
762 if(fEventMCProcessType[i]<0)
763 continue;
764 if (processtype==fEventMCProcessType[i])
765 return fEventMCProcessTypeIncluded;
766 }
767 return !fEventMCProcessTypeIncluded;
768}
769//_____________________________________________________________________________________________________________
770void AliSpectraBothEventCuts::SetNMCProcessType(Int_t flag)
771{
772 fNMCProcessType=flag;
773 if(fEventMCProcessType)
774 delete [] fEventMCProcessType;
775
776 fEventMCProcessType= new Int_t[fNMCProcessType];
777 if(!fEventMCProcessType)
778 {
779 fNMCProcessType=-1;
780 return;
781 }
782 for(int i=0;i<fNMCProcessType;i++)
783 fEventMCProcessType[i]=-1.0;
784 return;
785}
786//________________________________________________________________________________________________________________________
787void AliSpectraBothEventCuts::AddMCProcessType(Int_t type,Int_t index)
788{
789 if(index<fNMCProcessType)
790 fEventMCProcessType[index]=type;
791}
792