some cleanup+ add diagnostic histograms
[u/mrichter/AliRoot.git] / CORRFW / AliCFEventRecCuts.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 "AliESDEvent.h"
30 #include "AliESDVertex.h"
31 #include "AliCFEventRecCuts.h"
32 ClassImp(AliCFEventRecCuts) 
33 //____________________________________________________________________
34 AliCFEventRecCuts::AliCFEventRecCuts() : 
35   AliCFCutBase(),
36   fNTracksMin(-1),
37   fNTracksMax(1000000),
38   fRequireVtxCuts(kFALSE),
39   fVtxXMax(1.e99),
40   fVtxYMax(1.e99),
41   fVtxZMax(1.e99),
42   fVtxXMin(-1.e99),
43   fVtxYMin(-1.e99),
44   fVtxZMin(-1.e99),
45   fVtxXResMax(1.e99),
46   fVtxYResMax(1.e99),
47   fVtxZResMax(1.e99),
48   fBitMap(0x0),
49   fhNBinsNTracks(0),
50   fhBinLimNTracks(0),
51   fhNBinsVtxPosX(0),
52   fhBinLimVtxPosX(0),
53   fhNBinsVtxPosY(0),
54   fhBinLimVtxPosY(0),
55   fhNBinsVtxPosZ(0),
56   fhBinLimVtxPosZ(0),
57   fhNBinsVtxResX(0),
58   fhBinLimVtxResX(0),
59   fhNBinsVtxResY(0),
60   fhBinLimVtxResY(0),
61   fhNBinsVtxResZ(0),
62   fhBinLimVtxResZ(0)
63 {
64   //
65   //ctor
66   //
67   fBitMap=new TBits(0);
68   Initialise();
69 }
70
71 //____________________________________________________________________
72 AliCFEventRecCuts::AliCFEventRecCuts(Char_t* name, Char_t* title) : 
73   AliCFCutBase(name,title),
74   fNTracksMin(-1),
75   fNTracksMax(1000000),
76   fRequireVtxCuts(kFALSE),
77   fVtxXMax(1.e99),
78   fVtxYMax(1.e99),
79   fVtxZMax(1.e99),
80   fVtxXMin(-1.e99),
81   fVtxYMin(-1.e99),
82   fVtxZMin(-1.e99),
83   fVtxXResMax(1.e99),
84   fVtxYResMax(1.e99),
85   fVtxZResMax(1.e99),
86   fBitMap(0x0),
87   fhNBinsNTracks(0),
88   fhBinLimNTracks(0),
89   fhNBinsVtxPosX(0),
90   fhBinLimVtxPosX(0),
91   fhNBinsVtxPosY(0),
92   fhBinLimVtxPosY(0),
93   fhNBinsVtxPosZ(0),
94   fhBinLimVtxPosZ(0),
95   fhNBinsVtxResX(0),
96   fhBinLimVtxResX(0),
97   fhNBinsVtxResY(0),
98   fhBinLimVtxResY(0),
99   fhNBinsVtxResZ(0),
100   fhBinLimVtxResZ(0)
101  {
102   //
103   //ctor
104   //
105   fBitMap=new TBits(0);
106   Initialise();
107  }
108
109 //____________________________________________________________________
110 AliCFEventRecCuts::AliCFEventRecCuts(const AliCFEventRecCuts& c) : 
111   AliCFCutBase(c),
112   fNTracksMin(c.fNTracksMin),
113   fNTracksMax(c.fNTracksMax),
114   fRequireVtxCuts(c.fRequireVtxCuts),
115   fVtxXMax(c.fVtxXMax),
116   fVtxYMax(c.fVtxYMax),
117   fVtxZMax(c.fVtxZMax),
118   fVtxXMin(c.fVtxXMin),
119   fVtxYMin(c.fVtxYMin),
120   fVtxZMin(c.fVtxZMin),
121   fVtxXResMax(c.fVtxXResMax),
122   fVtxYResMax(c.fVtxYResMax),
123   fVtxZResMax(c.fVtxZResMax),
124   fBitMap(c.fBitMap),
125   fhNBinsNTracks(c.fhNBinsNTracks),
126   fhBinLimNTracks(c.fhBinLimNTracks),
127   fhNBinsVtxPosX(c.fhNBinsVtxPosX),
128   fhBinLimVtxPosX(c.fhBinLimVtxPosX),
129   fhNBinsVtxPosY(c.fhNBinsVtxPosY),
130   fhBinLimVtxPosY(c.fhBinLimVtxPosY),
131   fhNBinsVtxPosZ(c.fhNBinsVtxPosZ),
132   fhBinLimVtxPosZ(c.fhBinLimVtxPosZ),
133   fhNBinsVtxResX(c.fhNBinsVtxResX),
134   fhBinLimVtxResX(c.fhBinLimVtxResX),
135   fhNBinsVtxResY(c.fhNBinsVtxResY),
136   fhBinLimVtxResY(c.fhBinLimVtxResY),
137   fhNBinsVtxResZ(c.fhNBinsVtxResZ),
138   fhBinLimVtxResZ(c.fhBinLimVtxResZ) 
139 {
140   //
141   //copy constructor
142   //
143   for (Int_t i=0; i<c.kNCuts; i++){
144     for (Int_t j=0; j<c.kNStepQA; j++){
145       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
146     }
147   }
148
149 }
150
151 //____________________________________________________________________
152 AliCFEventRecCuts::~AliCFEventRecCuts() {
153   //
154   //dtor
155   //
156
157   for (Int_t i=0; i<kNCuts; i++){
158     for (Int_t j=0; j<kNStepQA; j++){
159       if(fhQA[i][j]) delete fhQA[i][j];
160     }
161   }
162
163   if(fhBinLimNTracks)delete fhBinLimNTracks;
164   if(fhBinLimVtxPosX)delete fhBinLimVtxPosX;
165   if(fhBinLimVtxPosY)delete fhBinLimVtxPosY;
166   if(fhBinLimVtxPosZ)delete fhBinLimVtxPosZ;
167   if(fhBinLimVtxResX)delete fhBinLimVtxResX;
168   if(fhBinLimVtxResY)delete fhBinLimVtxResY;
169   if(fhBinLimVtxResZ)delete fhBinLimVtxResZ;
170
171   if(fBitMap)delete fBitMap;
172
173 }
174
175 //_____________________________________________________________________________
176 void AliCFEventRecCuts::Init() {
177   //
178   // initialises all QA histograms
179   //
180   if(fIsQAOn)
181     DefineHistograms();
182 }
183
184 //_____________________________________________________________________________
185 void AliCFEventRecCuts::Initialise()
186 {
187
188   //
189   //initialization
190   //
191
192   //
193   // sets pointers to histos to zero
194   //
195
196   for(Int_t i=0; i<kNCuts; i++){
197     for(Int_t j =0; j<kNStepQA; j++){
198       fhQA[i][j]=0x0;
199     }
200   }
201
202   //set default bin number/ranges for QA histograms
203
204   SetHistogramBins(kNTracks,23,-0.5,22.5);
205   SetHistogramBins(kVtxPosX,100,-5,5);
206   SetHistogramBins(kVtxPosY,100,-5,5);
207   SetHistogramBins(kVtxPosZ,100,-50,50);
208   SetHistogramBins(kVtxResX,100,-1,1);
209   SetHistogramBins(kVtxResY,100,-1,1);
210   SetHistogramBins(kVtxResZ,100,-1,1);
211
212 }
213
214 //____________________________________________________________________
215 AliCFEventRecCuts& AliCFEventRecCuts::operator=(const AliCFEventRecCuts& c)
216 {
217   //
218   // Assignment operator
219   //
220   if (this != &c) {
221     AliCFCutBase::operator=(c) ;
222     fNTracksMin=c.fNTracksMin;
223     fNTracksMax=c.fNTracksMax;
224     fRequireVtxCuts=c.fRequireVtxCuts;
225     fVtxXMax=c.fVtxXMax;
226     fVtxYMax=c.fVtxYMax;
227     fVtxZMax=c.fVtxZMax;
228     fVtxXMin=c.fVtxXMin;
229     fVtxYMin=c.fVtxYMin;
230     fVtxZMin=c.fVtxZMin;
231     fVtxXResMax=c.fVtxXResMax;
232     fVtxYResMax=c.fVtxYResMax;
233     fVtxZResMax=c.fVtxZResMax;
234     fBitMap=c.fBitMap;
235     fhNBinsNTracks=c.fhNBinsNTracks;
236     fhBinLimNTracks=c.fhBinLimNTracks;
237     fhNBinsVtxPosX=c.fhNBinsVtxPosX;
238     fhBinLimVtxPosX=c.fhBinLimVtxPosX;
239     fhNBinsVtxPosY=c.fhNBinsVtxPosY;
240     fhBinLimVtxPosY=c.fhBinLimVtxPosY;
241     fhNBinsVtxPosZ=c.fhNBinsVtxPosZ;
242     fhBinLimVtxPosZ=c.fhBinLimVtxPosZ;
243     fhNBinsVtxResX=c.fhNBinsVtxResX;
244     fhBinLimVtxResX=c.fhBinLimVtxResX;
245     fhNBinsVtxResY=c.fhNBinsVtxResY;
246     fhBinLimVtxResY=c.fhBinLimVtxResY;
247     fhNBinsVtxResZ=c.fhNBinsVtxResZ;
248     fhBinLimVtxResZ=c.fhBinLimVtxResZ;
249   }
250
251   for (Int_t i=0; i<c.kNCuts; i++){
252     for (Int_t j=0; j<c.kNStepQA; j++){
253       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
254     }
255   }
256
257
258   return *this ;
259 }
260
261 //____________________________________________________________________
262 Bool_t AliCFEventRecCuts::IsSelected(TObject* obj) {
263   //
264   //Check if the requested cuts are passed
265   //
266
267   TBits *bitmap = SelectionBitMap(obj);
268
269   Bool_t isSelected = kTRUE;
270
271   for (UInt_t icut=0; icut<bitmap->GetNbits();icut++)
272         if(!bitmap->TestBitNumber(icut)) isSelected = kFALSE;
273
274   return isSelected;
275
276 }
277
278 //____________________________________________________________________
279 TBits *AliCFEventRecCuts::SelectionBitMap(TObject* obj) {
280   //
281   //cut on the number of charged tracks and on the event vertex.
282   //so far specific to AliESDEvents
283   //
284
285   //Check if the requested cuts are passed and return a bitmap
286   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
287   AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
288   if ( !esd ) return fBitMap ;
289
290   //now start checking the cuts,
291   //first assume the event will be accepted: 
292   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
293
294   //Number of charged tracks:
295   Int_t nTracks = esd->GetNumberOfTracks();
296   if(nTracks<fNTracksMin || nTracks>fNTracksMax)
297     fBitMap->SetBitNumber(0,kFALSE); 
298   
299   if(fRequireVtxCuts){
300     const AliESDVertex* vtxESD = esd->GetVertex();
301     if(!vtxESD){
302       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
303       return fBitMap;
304     }
305     // Require the vertex to have been reconstructed successfully
306     if (strcmp(vtxESD->GetName(), "default")==0){
307       AliWarning(Form(" No reconstructed vertex found, skip event"));    
308       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
309       return fBitMap;
310     }    
311     // Pick up the position and uncertainties
312     
313     Double_t vtxPos[3];
314     vtxPos[0] = vtxESD->GetXv();
315     vtxPos[1] = vtxESD->GetYv();
316     vtxPos[2] = vtxESD->GetZv();
317     
318     Double_t vtxRes[3];
319     vtxRes[0] = vtxESD->GetXRes();
320     vtxRes[1] = vtxESD->GetYRes();
321     vtxRes[2] = vtxESD->GetZRes();
322  
323     // Apply the cut
324     
325     if (vtxPos[0]>fVtxXMax || vtxPos[0]<fVtxXMin)
326       fBitMap->SetBitNumber(1,kFALSE); 
327     if (vtxPos[1]>fVtxYMax || vtxPos[1]<fVtxYMin)
328       fBitMap->SetBitNumber(2,kFALSE); 
329     if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
330       fBitMap->SetBitNumber(3,kFALSE); 
331     if (vtxRes[0]==0 || vtxRes[0]>fVtxXResMax)
332       fBitMap->SetBitNumber(4,kFALSE); 
333     if (vtxRes[1]==0 || vtxRes[1]>fVtxYResMax)
334       fBitMap->SetBitNumber(5,kFALSE); 
335     if (vtxRes[2]==0 || vtxRes[2]>fVtxZResMax)
336       fBitMap->SetBitNumber(6,kFALSE); 
337   }  
338   return fBitMap;
339 }
340
341 //_____________________________________________________________________________
342 void AliCFEventRecCuts::GetBitMap(TObject* obj, TBits *bitmap) {
343   //
344   // retrieve the pointer to the bitmap
345   //
346   bitmap = SelectionBitMap(obj);
347
348 }
349
350 //_____________________________________________________________________________
351 void AliCFEventRecCuts::FillHistograms(TObject* obj, Bool_t b)
352 {
353   //
354   // fill the QA histograms
355   //
356
357   if(!fIsQAOn) return;
358   // cast TObject into VParticle
359   AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
360   if (!esd ) return  ;
361
362   // index = 0: fill histograms before cuts
363   // index = 1: fill histograms after cuts
364   Int_t index = -1;
365   index = ((b) ? 1 : 0);
366
367
368   //number of charged tracks:
369   Int_t nTracks = esd->GetNumberOfTracks();
370   fhQA[kNTracks][index]->Fill(nTracks);
371
372   //look at vertex parameters:
373   const AliESDVertex* vtxESD = esd->GetVertex();
374   if(!vtxESD)return;
375   // Require the vertex to have been reconstructed successfully
376   if (strcmp(vtxESD->GetName(), "default")==0)return;
377   // vertex position and uncertainties
378   fhQA[kVtxPosX][index]->Fill(vtxESD->GetXv());
379   fhQA[kVtxPosY][index]->Fill(vtxESD->GetYv());
380   fhQA[kVtxPosZ][index]->Fill(vtxESD->GetZv());
381   fhQA[kVtxResX][index]->Fill(vtxESD->GetXRes());
382   fhQA[kVtxResY][index]->Fill(vtxESD->GetYRes());
383   fhQA[kVtxResZ][index]->Fill(vtxESD->GetZRes());
384   
385 }
386
387 //_____________________________________________________________________________
388 void AliCFEventRecCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
389 {
390   //
391   // QA histogram axis parameters
392   // variable bin size:user inputs nbins and the vector of bin limits
393   //
394
395   switch(index){
396   case kNTracks:
397     fhNBinsNTracks=nbins;
398     fhBinLimNTracks=new Double_t[nbins+1];
399     for(Int_t i=0;i<nbins+1;i++)fhBinLimNTracks[i]=bins[i];
400     break;
401   case kVtxPosX:
402     fhNBinsVtxPosX=nbins;
403     fhBinLimVtxPosX=new Double_t[nbins+1];
404     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosX[i]=bins[i];
405     break;
406   case kVtxPosY:
407     fhNBinsVtxPosY=nbins;
408     fhBinLimVtxPosY=new Double_t[nbins+1];
409     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosY[i]=bins[i];
410     break;
411   case kVtxPosZ:
412     fhNBinsVtxPosZ=nbins;
413     fhBinLimVtxPosZ=new Double_t[nbins+1];
414     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosZ[i]=bins[i];
415     break;
416   case kVtxResX:
417     fhNBinsVtxResX=nbins;
418     fhBinLimVtxResX=new Double_t[nbins+1];
419     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResX[i]=bins[i];
420     break;
421   case kVtxResY:
422     fhNBinsVtxResY=nbins;
423     fhBinLimVtxResY=new Double_t[nbins+1];
424     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResY[i]=bins[i];
425     break;
426   case kVtxResZ:
427     fhNBinsVtxResZ=nbins;
428     fhBinLimVtxResZ=new Double_t[nbins+1];
429     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResZ[i]=bins[i];
430     break;
431   }
432
433 }
434
435 //_____________________________________________________________________________
436 void AliCFEventRecCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
437 {
438   //
439   // QA histogram axis parameters
440   // fixed bin size: user inputs nbins, xmin and xmax
441   //
442   switch(index){
443   case kNTracks:
444     fhNBinsNTracks=nbins;
445     fhBinLimNTracks=new Double_t[nbins+1];
446     for(Int_t i=0;i<nbins+1;i++)fhBinLimNTracks[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
447     break;
448   case kVtxPosX:
449     fhNBinsVtxPosX=nbins;
450     fhBinLimVtxPosX=new Double_t[nbins+1];
451     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosX[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
452     break;
453   case kVtxPosY:
454     fhNBinsVtxPosY=nbins;
455     fhBinLimVtxPosY=new Double_t[nbins+1];
456     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosY[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
457     break;
458   case kVtxPosZ:
459     fhNBinsVtxPosZ=nbins;
460     fhBinLimVtxPosZ=new Double_t[nbins+1];
461     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxPosZ[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
462     break;
463   case kVtxResX:
464     fhNBinsVtxResX=nbins;
465     fhBinLimVtxResX=new Double_t[nbins+1];
466     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResX[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
467     break;
468   case kVtxResY:
469     fhNBinsVtxResY=nbins;
470     fhBinLimVtxResY=new Double_t[nbins+1];
471     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResY[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
472     break;
473   case kVtxResZ:
474     fhNBinsVtxResZ=nbins;
475     fhBinLimVtxResZ=new Double_t[nbins+1];
476     for(Int_t i=0;i<nbins+1;i++)fhBinLimVtxResZ[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
477     break;
478   }
479 }
480
481 //_____________________________________________________________________________
482  void AliCFEventRecCuts::DefineHistograms() {
483   //
484   // histograms for cut variables
485   //
486   Int_t color = 2;
487
488   if(!fIsQAOn) {
489     AliInfo(Form("No QA histos requested, Please first set the QA flag on!"));
490     return;
491   }  
492   
493   // book QA histograms
494
495   Char_t str[256];
496   for (Int_t i=0; i<kNStepQA; i++) {
497     if (i==0) sprintf(str," ");
498     else sprintf(str,"_cut");
499
500     fhQA[kNTracks][i]   = new  TH1F(Form("%s_NTracks%s",GetName(),str),                 "",fhNBinsNTracks,fhBinLimNTracks);
501     fhQA[kVtxPosX][i]   = new  TH1F(Form("%s_Vtx_Pos_X%s",GetName(),str),               "",fhNBinsVtxPosX,fhBinLimVtxPosX);
502     fhQA[kVtxPosY][i]   = new  TH1F(Form("%s_Vtx_Pos_Y%s",GetName(),str),               "",fhNBinsVtxPosY,fhBinLimVtxPosY);
503     fhQA[kVtxPosZ][i]   = new  TH1F(Form("%s_Vtx_Pos_Z%s",GetName(),str),               "",fhNBinsVtxPosZ,fhBinLimVtxPosZ);
504
505     fhQA[kVtxResX][i]   = new  TH1F(Form("%s_Vtx_Res_X%s",GetName(),str),               "",fhNBinsVtxResX,fhBinLimVtxResX);
506     fhQA[kVtxResY][i]   = new  TH1F(Form("%s_Vtx_Res_Y%s",GetName(),str),               "",fhNBinsVtxResY,fhBinLimVtxResY);
507     fhQA[kVtxResZ][i]   = new  TH1F(Form("%s_Vtx_Res_Z%s",GetName(),str),               "",fhNBinsVtxResZ,fhBinLimVtxResZ);
508  
509     fhQA[kNTracks][i]   ->SetXTitle("Number of ESD tracks");
510     fhQA[kVtxPosX][i]   ->SetXTitle("Vertex Position X (cm)");
511     fhQA[kVtxPosY][i]   ->SetXTitle("Vertex Position Y (cm)");
512     fhQA[kVtxPosZ][i]   ->SetXTitle("Vertex Position Z (cm)");
513     fhQA[kVtxResX][i]   ->SetXTitle("Vertex Resolution X (cm)");
514     fhQA[kVtxResY][i]   ->SetXTitle("Vertex Resolution Y (cm)");
515     fhQA[kVtxResZ][i]   ->SetXTitle("Vertex Resolution Z (cm)");
516
517   }
518
519   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
520
521 }
522
523 //_____________________________________________________________________________
524 void AliCFEventRecCuts::AddQAHistograms(TList *list) const {
525   //
526   // saves the histograms in a TList
527   //
528   if(!fIsQAOn) return;  
529
530   for (Int_t j=0; j<kNStepQA; j++) {
531     for(Int_t i=0; i<kNCuts; i++)
532         list->Add(fhQA[i][j]);
533   }
534 }