coverity fix
[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"
39//#include "AliSpectraBothHistoManager.h"
40#include <iostream>
41
42using namespace std;
43
44ClassImp(AliSpectraBothEventCuts)
45
46AliSpectraBothEventCuts::AliSpectraBothEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0),fAODEvent(AliSpectraBothTrackCuts::kAODobject), fTrackBits(0),fIsMC(0),fCentFromV0(0), fUseCentPatchAOD049(0), fUseSDDPatchforLHC11a(kDoNotCheckforSDD),fTrackCuts(0),
47fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(0),fMaxChi2perNDFforVertex(0),
48fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
49,fHistoEP(0)
50{
51 // Constructor
52 fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
53 fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",500,-15,15);
54 fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",500,-15,15);
55 fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
56 fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
57 fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
58 //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
59 //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
60 fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
61 fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
62 fCentralityCutMin = 0.0; // default value of centrality cut minimum, 0 ~ no cut
63 fCentralityCutMax = 10000.0; // default value of centrality cut maximum, ~ no cut
64 // fQVectorPosCutMin=0.0;
65 // fQVectorPosCutMax=10000.0;
66 // fQVectorNegCutMin=0.0;
67 // fQVectorNegCutMax=10000.0;
68 fQVectorCutMin=0.0;
69 fQVectorCutMax=10000.0;
70 fVertexCutMin=-10.0;
71 fVertexCutMax=10.0;
72 fMultiplicityCutMin=-1.0;
73 fMultiplicityCutMax=-1.0;
74 fTrackBits=1;
75 fCentFromV0=kFALSE;
76 fMaxChi2perNDFforVertex=-1;
77}
78
79//______________________________________________________
80Bool_t AliSpectraBothEventCuts::IsSelected(AliVEvent * aod,AliSpectraBothTrackCuts *trackcuts)
81{
82 // Returns true if Event Cuts are selected and applied
83 fAOD = aod;
84 fTrackCuts = trackcuts;
85 fHistoCuts->Fill(kProcessedEvents);
86 Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
87 if(!IsPhysSel)return IsPhysSel;
88 //loop on tracks, before event selection, filling QA histos
89 AliESDEvent* esdevent=0x0;
90 AliAODEvent* aodevent=0x0;
91 Bool_t isSDD=kFALSE;
92 TString nameoftrack(fAOD->ClassName());
93 if(!nameoftrack.CompareTo("AliESDEvent"))
94 {
95 fAODEvent=AliSpectraBothTrackCuts::kESDobject;
96
97 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
98 {
99 esdevent=dynamic_cast<AliESDEvent*>(fAOD);
737ad8ac 100 if(!esdevent)
94a8ed71 101 return kFALSE;
239a080a 102 if(esdevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
103 isSDD=kTRUE;
104 }
105 }
106 else if(!nameoftrack.CompareTo("AliAODEvent"))
107 {
108 fAODEvent=AliSpectraBothTrackCuts::kAODobject;
109 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
110 {
111 aodevent=dynamic_cast<AliAODEvent*>(fAOD);
737ad8ac 112 if(!aodevent)
94a8ed71 113 return kFALSE;
239a080a 114 if(aodevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
115 isSDD=kTRUE;
116 }
117 }
118 else
119 return false;
120 if(fUseSDDPatchforLHC11a==kwithSDD&&isSDD==kFALSE)
121 return false;
122 if(fUseSDDPatchforLHC11a==kwithoutSDD&&isSDD==kTRUE)
123 return false;
124
125
126
127 fHistoCuts->Fill(kPhysSelEvents);
128
129
130
131
132 const AliVVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
133 if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
134 fIsSelected =kFALSE;
135 if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut() && CheckVtxChi2perNDF()){ //selection on vertex and Centrality
136 fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
137 }
138 if(fIsSelected){
139 fHistoCuts->Fill(kAcceptedEvents);
140 if(vertex)
141 fHistoVtxAftSel->Fill(vertex->GetZ());
142 }
143 Int_t Nch=0;
144 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
145 AliVTrack* track =dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
146 /* if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
147 track=dynamic_cast<AliVTrack*>(esdevent->GetTrack(iTracks));
148 else if (fAODEvent==AliSpectraBothTrackCuts::kAODobject)
149 track=dynamic_cast<AliVTrack*>(aodevent->GetTrack(iTracks));
150 else return false;*/
151
152 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
153 fHistoEtaBefSel->Fill(track->Eta());
154 if(fIsSelected){
155 fHistoEtaAftSel->Fill(track->Eta());
156 Nch++;
157 }
158 }
159 //Printf("NCHARGED_EvSel : %d",Nch);
160 if(fIsSelected)fHistoNChAftSel->Fill(Nch);
161 return fIsSelected;
162}
163
164//______________________________________________________
165Bool_t AliSpectraBothEventCuts::CheckVtxRange()
166{
167 // reject events outside of range
168 const AliVVertex * vertex = fAOD->GetPrimaryVertex();
169 //when moving to 2011 w├Čone has to add a cut using SPD vertex.
170 //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
171 if (!vertex)
172 {
173 fHistoCuts->Fill(kVtxNoEvent);
174 return kFALSE;
175 }
176 if(vertex->GetNContributors()<1)
177 {
178 fHistoCuts->Fill(kVtxNoEvent);
179 return kFALSE;
180
181 }
182 TString tmp(vertex->GetTitle());
183 if(tmp.Contains("TPC"))
184 {
185 fHistoCuts->Fill(kVtxNoEvent);
186 return kFALSE;
187 }
188 if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
189 {
190 return kTRUE;
191 }
192 fHistoCuts->Fill(kVtxRange);
193 return kFALSE;
194}
195
196//______________________________________________________
197Bool_t AliSpectraBothEventCuts::CheckCentralityCut()
198{
199 // Check centrality cut
200 if ( fCentralityCutMax<0.0 && fCentralityCutMin<0.0 ) return kTRUE;
201 Double_t cent=0;
202 TString CentEstimator="TRK";
203 if(fCentFromV0)CentEstimator="V0M";
204 if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile(CentEstimator.Data());
205 else cent=ApplyCentralityPatchAOD049();
206
207 if ( (cent <= fCentralityCutMax) && (cent >= fCentralityCutMin) ) return kTRUE;
208 fHistoCuts->Fill(kVtxCentral);
209
210 return kFALSE;
211}
212
213//______________________________________________________
214Bool_t AliSpectraBothEventCuts::CheckMultiplicityCut()
215{
216 // Check multiplicity cut
217if(fMultiplicityCutMin<0.0 && fMultiplicityCutMax<0.0)
218 return kTRUE;
219 Int_t Ncharged=0;
220 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
221 AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
222
223 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
224
225 Ncharged++;
226 }
227 //Printf("NCHARGED_cut : %d",Ncharged);
228 if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
229
230 return kFALSE;
231}
232
233//______________________________________________________
234Bool_t AliSpectraBothEventCuts::CheckQVectorCut()
235{
236 if(fQVectorCutMin<0.0 && fQVectorCutMax<0.0)return kTRUE;
237 // Check qvector cut
238 /// FIXME: Q vector
239 // //Selection on QVector, before ANY other selection on the event
240 // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
241 // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
242 // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
243
244 // Int_t multPos = 0;
245 // Int_t multNeg = 0;
246 // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
247 // AliAODTrack* aodTrack = fAOD->GetTrack(iT);
248 // if (!aodTrack->TestFilterBit(fTrackBits)) continue;
249 // if (aodTrack->Eta() >= 0){
250 // multPos++;
251 // Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
252 // Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
253 // } else {
254 // multNeg++;
255 // Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
256 // Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
257 // }
258 // }
259 // Double_t qPos=-999;
260 // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
261 // Double_t qNeg=-999;
262 // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
263 //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
264
265 Double_t qxEPVZERO = 0, qyEPVZERO = 0;
266 Double_t qVZERO = -999;
267 Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
268
269 qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
270 if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
271 fHistoQVector->Fill(qVZERO);
272 fHistoEP->Fill(psi);
273
274 fHistoCuts->Fill(kQVector);
275 // fHistoQVectorPos->Fill(qPos);
276 // fHistoQVectorNeg->Fill(qNeg);
277 return kTRUE;
278}
279//____________________________________________________________
280 Bool_t AliSpectraBothEventCuts::CheckVtxChi2perNDF()
281 {
282 if(fMaxChi2perNDFforVertex<0)
283 return kTRUE;
284 const AliVVertex * vertex = fAOD->GetPrimaryVertex();
285 if(TMath::Abs(vertex->GetChi2perNDF())>fMaxChi2perNDFforVertex)
286 return kFALSE;
287 return kTRUE;
288 }
289
290
291
292//______________________________________________________
293void AliSpectraBothEventCuts::PrintCuts()
294{
295 // print info about event cuts
296 cout << "Event Stats" << endl;
297 cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
298 cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
299 cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
300 cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
301 cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
302 cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
303 cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
304 cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
305 // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
306 // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
307 cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
308 cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
309 cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
310 cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
311}
312//______________________________________________________
313
314Double_t AliSpectraBothEventCuts::ApplyCentralityPatchAOD049()
315{
316 //
317 //Apply centrality patch for AOD049 outliers
318 //
319 // if (fCentralityType!="V0M") {
320 // AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
321 // return -999.0;
322 // }
323 AliCentrality *centrality = fAOD->GetCentrality();
324 if (!centrality) {
325 AliWarning("Cannot get centrality from AOD event.");
326 return -999.0;
327 }
328
329 Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
330 /*
331 Bool_t isSelRun = kFALSE;
332 Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
333 if(cent<0){
334 Int_t quality = centrality->GetQuality();
335 if(quality<=1){
336 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
337 } else {
338 Int_t runnum=aodEvent->GetRunNumber();
339 for(Int_t ir=0;ir<5;ir++){
340 if(runnum==selRun[ir]){
341 isSelRun=kTRUE;
342 break;
343 }
344 }
345 if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
346 }
347 }
348 */
349 if(cent>=0.0) {
350 Float_t v0 = 0.0;
351 AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
352 AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
353 v0+=aodV0->GetMTotV0A();
354 v0+=aodV0->GetMTotV0C();
355 if ( (cent==0) && (v0<19500) ) {
356 AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
357 return -999.0;
358 }
359 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
360 Float_t val = 1.30552 + 0.147931 * v0;
361
362 Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
363 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
364 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
365 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
366 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
367 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
368 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
369 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
370 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
371 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
372 };
373
374 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
375 AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
376 return -999.0;
377 }
378 } else {
379 //force it to be -999. whatever the negative value was
380 cent = -999.;
381 }
382 return cent;
383}
384
385//______________________________________________________
386
387
388Long64_t AliSpectraBothEventCuts::Merge(TCollection* list)
389{
390 // Merge a list of AliSpectraBothEventCuts objects with this.
391 // Returns the number of merged objects (including this).
392
393 // AliInfo("Merging");
394
395 if (!list)
396 return 0;
397
398 if (list->IsEmpty())
399 return 1;
400
401 TIterator* iter = list->MakeIterator();
402 TObject* obj;
403
404 // collections of all histograms
405 TList collections;//FIXME we should use only 1 collection
406 TList collections_histoVtxBefSel;
407 TList collections_histoVtxAftSel;
408 TList collections_histoEtaBefSel;
409 TList collections_histoEtaAftSel;
410 TList collections_histoNChAftSel;
411 // TList collections_histoQVectorPos;
412 // TList collections_histoQVectorNeg;
413 TList collections_histoQVector;
414 TList collections_histoEP;
415
416 Int_t count = 0;
417
418 while ((obj = iter->Next())) {
419 AliSpectraBothEventCuts* entry = dynamic_cast<AliSpectraBothEventCuts*> (obj);
420 if (entry == 0)
421 continue;
422
423 TH1I * histo = entry->GetHistoCuts();
424 collections.Add(histo);
425 TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
426 collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
427 TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
428 collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
429 TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
430 collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
431 TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
432 collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
433 TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
434 collections_histoNChAftSel.Add(histo_histoNChAftSel);
435 // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
436 // collections_histoQVectorPos.Add(histo_histoQVectorPos);
437 // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
438 // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
439 TH1F * histo_histoQVector = entry->GetHistoQVector();
440 collections_histoQVector.Add(histo_histoQVector);
441 TH1F * histo_histoEP = entry->GetHistoEP();
442 collections_histoEP.Add(histo_histoEP);
443 count++;
444 }
445
446 fHistoCuts->Merge(&collections);
447 fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
448 fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
449 fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
450 fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
451 fHistoNChAftSel->Merge(&collections_histoNChAftSel);
452 // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
453 // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
454 fHistoQVector->Merge(&collections_histoQVector);
455 fHistoEP->Merge(&collections_histoEP);
456
457
458 delete iter;
459
460 return count+1;
461}
462