Fix Coverity defect
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEextraEventCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 // Cut on the Event at reconstructed level: for the moment 
16 // the requirements on the number of charged tracks and on 
17 // the vertex position and resolution are implemented
18 // The argument of IsSelected member function (passed object) is cast into 
19 // an AliESDEvent. In the future may be modified to use AliVEvent interface
20 // and include more cut variables.
21 // The class derives from AliCFCutBase
22 // Author:S.Arcelli Silvia.Arcelli@cern.ch
23 //
24 //
25 #include "TH1F.h"
26 #include "TList.h"
27 #include "TBits.h"
28 #include "AliLog.h"
29 #include "AliAnalysisUtils.h"
30 #include "AliESDEvent.h"
31 #include "AliESDVertex.h"
32 #include "AliAODEvent.h"
33 #include "AliAODVertex.h"
34 #include "AliHFEextraEventCuts.h"
35 ClassImp(AliHFEextraEventCuts) 
36 //____________________________________________________________________
37 AliHFEextraEventCuts::AliHFEextraEventCuts() : 
38   AliCFCutBase(),
39   fAnalysisUtils(NULL),
40   fRequireVtxCuts(kFALSE),
41   fVtxZMax(1.e99),
42   fVtxZMin(-1.e99),
43   fVtxNCtrbMin(0),
44   fVtxMixed(0),
45   fVtxSPD(0),
46   fCheckCorrelationSPDVtx(0),
47   fVtxResolution(kFALSE),
48   fPApileupCut(kFALSE),
49   fBitMap(0x0)
50 {
51   //
52   //ctor
53   //
54 }
55
56 //____________________________________________________________________
57 AliHFEextraEventCuts::AliHFEextraEventCuts(Char_t* name, Char_t* title) : 
58   AliCFCutBase(name,title),
59   fAnalysisUtils(NULL),
60   fRequireVtxCuts(kFALSE),
61   fVtxZMax(1.e99),
62   fVtxZMin(-1.e99),
63   fVtxNCtrbMin(0),
64   fVtxMixed(0),
65   fVtxSPD(0),
66   fCheckCorrelationSPDVtx(0),
67   fVtxResolution(kFALSE),
68   fPApileupCut(kFALSE),
69   fBitMap(0x0)
70  {
71   //
72   //ctor
73   //
74   fBitMap=new TBits(0);
75   Initialise();
76   fAnalysisUtils = new AliAnalysisUtils;
77  }
78
79 //____________________________________________________________________
80 AliHFEextraEventCuts::AliHFEextraEventCuts(const AliHFEextraEventCuts& c) : 
81   AliCFCutBase(c),
82   fAnalysisUtils(NULL),
83   fRequireVtxCuts(c.fRequireVtxCuts),
84   fVtxZMax(c.fVtxZMax),
85   fVtxZMin(c.fVtxZMin),
86   fVtxNCtrbMin(c.fVtxNCtrbMin),
87   fVtxMixed(c.fVtxMixed),
88   fVtxSPD(c.fVtxSPD),
89   fCheckCorrelationSPDVtx(c.fCheckCorrelationSPDVtx),
90   fVtxResolution(c.fVtxResolution),
91   fPApileupCut(c.fPApileupCut),
92   fBitMap(c.fBitMap)
93 {
94   //
95   //copy constructor
96   //
97   for (Int_t i=0; i<c.kNCuts; i++){
98     for (Int_t j=0; j<c.kNStepQA; j++){
99       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
100     }
101   }
102   fAnalysisUtils = new AliAnalysisUtils;
103 }
104
105 //____________________________________________________________________
106 AliHFEextraEventCuts::~AliHFEextraEventCuts() {
107   //
108   //dtor
109   //
110
111   for (Int_t i=0; i<kNCuts; i++){
112     for (Int_t j=0; j<kNStepQA; j++){
113       if(fhQA[i][j]) delete fhQA[i][j];
114     }
115   }
116
117   if(fBitMap) delete fBitMap;
118   if(fAnalysisUtils) delete fAnalysisUtils;
119 }
120 //_____________________________________________________________________________
121 void AliHFEextraEventCuts::Initialise()
122 {
123
124   //
125   //initialization
126   //
127
128   //
129   // sets pointers to histos to zero
130   //
131
132   for(Int_t i=0; i<kNCuts; i++){
133     for(Int_t j =0; j<kNStepQA; j++){
134       fhQA[i][j]=0x0;
135     }
136   }
137 }
138
139 //____________________________________________________________________
140 AliHFEextraEventCuts& AliHFEextraEventCuts::operator=(const AliHFEextraEventCuts& c)
141 {
142   //
143   // Assignment operator
144   //
145   if (this != &c) {
146     AliCFCutBase::operator=(c) ;
147     fRequireVtxCuts=c.fRequireVtxCuts;
148     fVtxZMax=c.fVtxZMax;
149     fVtxZMin=c.fVtxZMin;
150     fVtxNCtrbMin=c.fVtxNCtrbMin;
151     fVtxMixed=c.fVtxMixed;
152     fVtxSPD=c.fVtxSPD;
153     fCheckCorrelationSPDVtx=c.fCheckCorrelationSPDVtx;
154     fVtxResolution = c.fVtxResolution;
155     fBitMap=c.fBitMap;
156     if(!fAnalysisUtils) fAnalysisUtils = new AliAnalysisUtils;
157     fPApileupCut = c.fPApileupCut;
158   }
159
160   for (Int_t i=0; i<c.kNCuts; i++){
161     for (Int_t j=0; j<c.kNStepQA; j++){
162       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
163     }
164   }
165
166
167   return *this ;
168 }
169
170 //____________________________________________________________________
171 Bool_t AliHFEextraEventCuts::IsSelected(TObject* obj) {
172   //
173   //Check if the requested cuts are passed
174   //
175
176
177   SelectionBitMap(obj);
178
179   if (fIsQAOn) FillHistograms(obj,0);
180   Bool_t isSelected = kTRUE;
181
182   for (UInt_t icut=0; icut<fBitMap->GetNbits();icut++)
183         if(!fBitMap->TestBitNumber(icut)) isSelected = kFALSE;
184
185   if (!isSelected) return kFALSE ;
186   if (fIsQAOn) FillHistograms(obj,1);
187   return kTRUE;
188
189 }
190 //____________________________________________________________________
191 void AliHFEextraEventCuts::SelectionBitMap(TObject* obj) {
192   //
193   //cut on the number of charged tracks and on the event vertex.
194   //so far specific to AliESDEvents
195   //
196
197   //Check if the requested cuts are passed and return a bitmap
198
199   AliVEvent *inputEvent = dynamic_cast<AliVEvent *>(obj);
200   if(!inputEvent){
201     AliDebug(1, "Not a virtual event");
202     for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
203     return;
204   }
205   const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
206                    *vtxSPD = GetPrimaryVertexSPD(inputEvent),
207                    *vtxPrim(NULL);
208   
209   //first assume the event will be accepted: 
210   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
211
212   if(fVtxMixed){
213     // Use mixed vertex: Prefer vertex with tracks, in case not available use SPD vertex
214     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
215     else if(vtxSPD && vtxSPD->GetNContributors() > 0) vtxPrim = vtxSPD;
216   } else if(fVtxSPD){
217     if(vtxSPD && vtxSPD->GetNContributors () > 0) vtxPrim = vtxSPD;
218   } else {
219     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
220   }
221   if(!vtxPrim){
222     // No primary vertex: Reject event
223           for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
224           AliWarning("Cannot get vertex, skipping event");
225           return;
226   }
227   // Check standard vertex cuts using the primary vertex 
228   if(fVtxZMin > - 1000. && fVtxZMax < 1000.){
229     // Primary vertex z cut required
230     // Pick up the position and uncertainties
231     Double_t vtxPos[3];
232     vtxPos[0] = vtxPrim->GetX();
233     vtxPos[1] = vtxPrim->GetY();
234     vtxPos[2] = vtxPrim->GetZ();
235     if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
236             fBitMap->SetBitNumber(kVtxPosZ,kFALSE); 
237   }
238   if(fVtxNCtrbMin > 0){  
239     // cut required if the cut value is set to something larger than 0
240     // same effect as setting the min. cut value to 0
241     Int_t nCtrb = vtxPrim->GetNContributors();
242     if (nCtrb<fVtxNCtrbMin)
243             fBitMap->SetBitNumber(kVtxNCtrb,kFALSE);
244   }
245
246   // check vertex correlation cut
247   if(fCheckCorrelationSPDVtx){
248     if(vtxTracks && vtxTracks->GetNContributors() && vtxSPD && vtxSPD->GetNContributors()){
249       // Both vertices available, check correlation
250       if(TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) >= 0.5) fBitMap->SetBitNumber(kCorrelation, kFALSE);
251     } else {
252       // No correlation available: set the cut to false
253       fBitMap->SetBitNumber(kCorrelation, kFALSE); 
254     }
255   }
256   // check cut on the SPD vertex resolution
257   if(fVtxResolution){
258     if(vtxSPD){
259       TString vtxTyp = vtxSPD->GetTitle();
260       Double_t cov[6]={0};
261       vtxSPD->GetCovarianceMatrix(cov);
262       Double_t zRes = TMath::Sqrt(cov[5]);
263       if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fBitMap->SetBitNumber(kResolution, kFALSE);
264     } else{
265       fBitMap->SetBitNumber(kResolution, kFALSE);
266     }
267   }
268   // check pA pileup cut
269   if(fPApileupCut){
270     if(fAnalysisUtils->IsPileUpEvent(inputEvent)) fBitMap->SetBitNumber(kpApileup,kFALSE);
271   }
272 }
273
274 //_____________________________________________________________________________
275 void AliHFEextraEventCuts::FillHistograms(TObject* obj, Bool_t b)
276 {
277   //
278   // fill the QA histograms
279   //
280
281   if(!fIsQAOn) return;
282   // cast TObject into VEvent
283   AliVEvent* inputEvent = dynamic_cast<AliESDEvent *>(obj);
284   if (!inputEvent) return;
285
286   // index = 0: fill histograms before cuts
287   // index = 1: fill histograms after cuts
288   Int_t index = -1;
289   index = ((b) ? 1 : 0);
290
291   // Obtain vertices
292   const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
293                    *vtxSPD = GetPrimaryVertexSPD(inputEvent),
294                    *vtxPrim(NULL);
295
296   //look at vertex parameters:
297   if(fVtxMixed){
298     if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
299     else if(vtxSPD && vtxSPD->GetNContributors()) vtxPrim = vtxSPD;
300   }
301   else {   
302     vtxPrim = vtxTracks; 
303   }
304   if(!vtxPrim)return;
305   // vertex position and uncertainties
306   fhQA[kVtxPosZ] [index]->Fill(vtxPrim->GetZ());
307   fhQA[kVtxNCtrb][index]->Fill(vtxPrim->GetNContributors());
308   // SPD Vertex Position Correlation 
309   if(vtxTracks && vtxSPD){
310     fhQA[kCorrelation][index]->Fill(vtxTracks->GetZ()-vtxSPD->GetZ());
311   }
312   if(vtxSPD){
313     Double_t cov[6]={0};
314     vtxSPD->GetCovarianceMatrix(cov);
315     fhQA[kResolution][index]->Fill(TMath::Sqrt(cov[5]));
316   }
317   double isPileup =0 ;
318   if(fAnalysisUtils->IsPileUpEvent(inputEvent)) isPileup = 1,
319   fhQA[kpApileup][index]->Fill(isPileup);
320 }
321
322 //____________________________________________________________________
323 void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
324 {
325   //
326   //setting x-axis bin limits of QA histogram fhQA[index] 
327   // 
328
329   for(Int_t i=0;i<kNStepQA;i++){
330     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
331     return;
332     }
333     fhQA[index][i]->GetXaxis()->Set(nbins,bins);
334   }
335 }
336 //____________________________________________________________________
337 void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
338 {
339   //
340   //setting x-axis bins and range of QA histogram fhQA[index] 
341   // 
342
343   for(Int_t i=0;i<kNStepQA;i++){
344     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
345     return;
346     }
347     fhQA[index][i]->GetXaxis()->Set(nbins,xmin,xmax);
348   }
349 }
350
351 //_____________________________________________________________________________
352  void AliHFEextraEventCuts::DefineHistograms() {
353   //
354   // histograms for cut variables
355   //
356   Int_t color = 2;
357
358   if(!fIsQAOn) {
359     AliInfo(Form("No QA histos requested, Please first set the QA flag on!"));
360     return;
361   }  
362   
363   // book QA histograms
364   Char_t str[5];
365   for (Int_t i=0; i<kNStepQA; i++) {
366     if (i==0) snprintf(str,5," ");
367     else snprintf(str,5,"_cut");
368
369     fhQA[kVtxPosZ][i]   = new  TH1F(Form("%s_Vtx_Pos_Z%s",GetName(),str),               "",200,-50.,50.);
370     fhQA[kVtxNCtrb][i]  = new  TH1F(Form("%s_Vtx_N_Ctrb%s",GetName(),str),              "",1000,0.,1000.);
371     fhQA[kCorrelation][i] = new TH1F(Form("%s_SPDCorrelation_%s",GetName(),str), "",200, -10., 10);
372     fhQA[kResolution][i] = new TH1F(Form("%s_SPDResolution_%s",GetName(), str), "", 100, 0., 1.); 
373     fhQA[kpApileup][i] = new TH1F(Form("%s_pApileup_%s",GetName(),str), "", 2, 0., 2.);
374  
375     fhQA[kVtxPosZ][i]   ->SetXTitle("Vertex Position Z (cm)");
376     fhQA[kVtxNCtrb][i]  ->SetXTitle("Number of contributors");
377   }
378
379   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
380
381 }
382
383 //_____________________________________________________________________________
384 void AliHFEextraEventCuts::AddQAHistograms(TList *qaList) {
385   //
386   // saves the histograms in a TList
387   //
388
389   DefineHistograms();
390
391   for (Int_t j=0; j<kNStepQA; j++) {
392     for(Int_t i=0; i<kNCuts; i++)
393         qaList->Add(fhQA[i][j]);
394   }
395 }
396
397 //_____________________________________________________________________________
398 const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexSPD(const AliVEvent * const inputEvent){
399   //
400   // Get SPD vertex from event
401   //
402   const AliVVertex *spdvtx(NULL);
403   const AliESDEvent *esd(NULL);
404   const AliAODEvent *aod(NULL);
405   if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
406     spdvtx = esd->GetPrimaryVertexSPD();
407   } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
408     spdvtx = aod->GetPrimaryVertexSPD();
409   }
410   return spdvtx;
411 }
412
413 //_____________________________________________________________________________
414 const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexTracks(const AliVEvent *const inputEvent){
415   // 
416   // Get Primary Vertex from tracks
417   //
418   const AliVVertex *trkvtx(NULL);
419   const AliESDEvent *esd(NULL);
420   const AliAODEvent *aod(NULL);
421   if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
422     trkvtx = esd->GetPrimaryVertexTracks();
423   } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
424     const AliVVertex *vtxTmp = aod->GetPrimaryVertex();
425     // check whether the primary vertex is the vertex from tracks
426     TString vtxTtl = vtxTmp->GetTitle();
427     if(vtxTtl.Contains("VertexerTracks")){
428       trkvtx = vtxTmp;
429     }
430   }
431   return trkvtx;
432 }
433