bug fix: data member of AliCFTrackIsPrimaryCuts changed from AliESD to AliESDEvent
[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   fVtxNCtrbMin(0),
49   fVtxNCtrbMax((Int_t)1.e9),
50   fVtxTPC(0),
51   fVtxSPD(0),
52   fBitMap(0x0)
53 {
54   //
55   //ctor
56   //
57   fBitMap=new TBits(0);
58   Initialise();
59 }
60
61 //____________________________________________________________________
62 AliCFEventRecCuts::AliCFEventRecCuts(Char_t* name, Char_t* title) : 
63   AliCFCutBase(name,title),
64   fNTracksMin(-1),
65   fNTracksMax(1000000),
66   fRequireVtxCuts(kFALSE),
67   fVtxXMax(1.e99),
68   fVtxYMax(1.e99),
69   fVtxZMax(1.e99),
70   fVtxXMin(-1.e99),
71   fVtxYMin(-1.e99),
72   fVtxZMin(-1.e99),
73   fVtxXResMax(1.e99),
74   fVtxYResMax(1.e99),
75   fVtxZResMax(1.e99),
76   fVtxNCtrbMin(0),
77   fVtxNCtrbMax((Int_t)1.e9),
78   fVtxTPC(0),
79   fVtxSPD(0),
80   fBitMap(0x0)
81  {
82   //
83   //ctor
84   //
85   fBitMap=new TBits(0);
86   Initialise();
87  }
88
89 //____________________________________________________________________
90 AliCFEventRecCuts::AliCFEventRecCuts(const AliCFEventRecCuts& c) : 
91   AliCFCutBase(c),
92   fNTracksMin(c.fNTracksMin),
93   fNTracksMax(c.fNTracksMax),
94   fRequireVtxCuts(c.fRequireVtxCuts),
95   fVtxXMax(c.fVtxXMax),
96   fVtxYMax(c.fVtxYMax),
97   fVtxZMax(c.fVtxZMax),
98   fVtxXMin(c.fVtxXMin),
99   fVtxYMin(c.fVtxYMin),
100   fVtxZMin(c.fVtxZMin),
101   fVtxXResMax(c.fVtxXResMax),
102   fVtxYResMax(c.fVtxYResMax),
103   fVtxZResMax(c.fVtxZResMax),
104   fVtxNCtrbMin(c.fVtxNCtrbMin),
105   fVtxNCtrbMax(c.fVtxNCtrbMax),
106   fVtxTPC(c.fVtxTPC),
107   fVtxSPD(c.fVtxSPD),
108   fBitMap(c.fBitMap)
109 {
110   //
111   //copy constructor
112   //
113   for (Int_t i=0; i<c.kNCuts; i++){
114     for (Int_t j=0; j<c.kNStepQA; j++){
115       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
116     }
117   }
118
119 }
120
121 //____________________________________________________________________
122 AliCFEventRecCuts::~AliCFEventRecCuts() {
123   //
124   //dtor
125   //
126
127   for (Int_t i=0; i<kNCuts; i++){
128     for (Int_t j=0; j<kNStepQA; j++){
129       if(fhQA[i][j]) delete fhQA[i][j];
130     }
131   }
132
133   if(fBitMap)delete fBitMap;
134
135 }
136 //_____________________________________________________________________________
137 void AliCFEventRecCuts::Initialise()
138 {
139
140   //
141   //initialization
142   //
143
144   //
145   // sets pointers to histos to zero
146   //
147
148   for(Int_t i=0; i<kNCuts; i++){
149     for(Int_t j =0; j<kNStepQA; j++){
150       fhQA[i][j]=0x0;
151     }
152   }
153 }
154
155 //____________________________________________________________________
156 AliCFEventRecCuts& AliCFEventRecCuts::operator=(const AliCFEventRecCuts& c)
157 {
158   //
159   // Assignment operator
160   //
161   if (this != &c) {
162     AliCFCutBase::operator=(c) ;
163     fNTracksMin=c.fNTracksMin;
164     fNTracksMax=c.fNTracksMax;
165     fRequireVtxCuts=c.fRequireVtxCuts;
166     fVtxXMax=c.fVtxXMax;
167     fVtxYMax=c.fVtxYMax;
168     fVtxZMax=c.fVtxZMax;
169     fVtxXMin=c.fVtxXMin;
170     fVtxYMin=c.fVtxYMin;
171     fVtxZMin=c.fVtxZMin;
172     fVtxXResMax=c.fVtxXResMax;
173     fVtxYResMax=c.fVtxYResMax;
174     fVtxZResMax=c.fVtxZResMax;
175     fVtxNCtrbMin=c.fVtxNCtrbMin;
176     fVtxNCtrbMax=c.fVtxNCtrbMax;
177     fVtxTPC=c.fVtxTPC;
178     fVtxSPD=c.fVtxSPD;
179     fBitMap=c.fBitMap;
180   }
181
182   for (Int_t i=0; i<c.kNCuts; i++){
183     for (Int_t j=0; j<c.kNStepQA; j++){
184       if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
185     }
186   }
187
188
189   return *this ;
190 }
191
192 //____________________________________________________________________
193 Bool_t AliCFEventRecCuts::IsSelected(TObject* obj) {
194   //
195   //Check if the requested cuts are passed
196   //
197
198
199   SelectionBitMap(obj);
200
201   if (fIsQAOn) FillHistograms(obj,0);
202   Bool_t isSelected = kTRUE;
203
204   for (UInt_t icut=0; icut<fBitMap->GetNbits();icut++)
205         if(!fBitMap->TestBitNumber(icut)) isSelected = kFALSE;
206
207   if (!isSelected) return kFALSE ;
208   if (fIsQAOn) FillHistograms(obj,1);
209   return kTRUE;
210
211 }
212 //____________________________________________________________________
213 void AliCFEventRecCuts::SelectionBitMap(TObject* obj) {
214   //
215   //cut on the number of charged tracks and on the event vertex.
216   //so far specific to AliESDEvents
217   //
218
219   //Check if the requested cuts are passed and return a bitmap
220   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
221   AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
222   if ( !esd ) return;
223
224   //now start checking the cuts,
225   //first assume the event will be accepted: 
226   for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
227
228   //Number of charged tracks:
229   Int_t nTracks = esd->GetNumberOfTracks();
230   if(nTracks<fNTracksMin || nTracks>fNTracksMax)
231     fBitMap->SetBitNumber(0,kFALSE); 
232   
233   if(fRequireVtxCuts){
234     const AliESDVertex* vtxESD = 0x0;
235     if      (fVtxTPC) vtxESD = esd->GetPrimaryVertexTPC() ;
236     else if (fVtxSPD) vtxESD = esd->GetPrimaryVertexSPD() ;
237     else              vtxESD = esd->GetPrimaryVertexTracks() ;
238     if(!vtxESD){
239       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
240       AliWarning("Cannot get vertex, skipping event");
241       return;
242     }
243     // Require the vertex to have been reconstructed successfully
244     if (strcmp(vtxESD->GetName(), "default")==0){
245       AliWarning(Form(" No reconstructed vertex found, skipping event"));    
246       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
247       return;
248     }    
249     // Pick up the position and uncertainties
250     
251     Double_t vtxPos[3];
252     vtxPos[0] = vtxESD->GetXv();
253     vtxPos[1] = vtxESD->GetYv();
254     vtxPos[2] = vtxESD->GetZv();
255     
256     Double_t vtxRes[3];
257     vtxRes[0] = vtxESD->GetXRes();
258     vtxRes[1] = vtxESD->GetYRes();
259     vtxRes[2] = vtxESD->GetZRes();
260  
261     Int_t nCtrb = vtxESD->GetNContributors();
262
263     // Apply the cut
264     
265     if (vtxPos[0]>fVtxXMax || vtxPos[0]<fVtxXMin)
266       fBitMap->SetBitNumber(1,kFALSE); 
267     if (vtxPos[1]>fVtxYMax || vtxPos[1]<fVtxYMin)
268       fBitMap->SetBitNumber(2,kFALSE); 
269     if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
270       fBitMap->SetBitNumber(3,kFALSE); 
271     if (vtxRes[0]==0 || vtxRes[0]>fVtxXResMax)
272       fBitMap->SetBitNumber(4,kFALSE); 
273     if (vtxRes[1]==0 || vtxRes[1]>fVtxYResMax)
274       fBitMap->SetBitNumber(5,kFALSE); 
275     if (vtxRes[2]==0 || vtxRes[2]>fVtxZResMax)
276       fBitMap->SetBitNumber(6,kFALSE); 
277     if (nCtrb<fVtxNCtrbMin || nCtrb>fVtxNCtrbMax)
278       fBitMap->SetBitNumber(7,kFALSE);
279
280   }  
281   return;
282 }
283
284 //_____________________________________________________________________________
285 void AliCFEventRecCuts::FillHistograms(TObject* obj, Bool_t b)
286 {
287   //
288   // fill the QA histograms
289   //
290
291   if(!fIsQAOn) return;
292   // cast TObject into VParticle
293   AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
294   if (!esd ) return  ;
295
296   // index = 0: fill histograms before cuts
297   // index = 1: fill histograms after cuts
298   Int_t index = -1;
299   index = ((b) ? 1 : 0);
300
301
302   //number of charged tracks:
303   Int_t nTracks = esd->GetNumberOfTracks();
304   fhQA[kNTracks][index]->Fill(nTracks);
305
306   //look at vertex parameters:
307   const AliESDVertex* vtxESD = 0x0;
308   if      (fVtxTPC) vtxESD = esd->GetPrimaryVertexTPC() ;
309   else if (fVtxSPD) vtxESD = esd->GetPrimaryVertexSPD() ;
310   else              vtxESD = esd->GetPrimaryVertexTracks() ;
311   if(!vtxESD)return;
312   // Require the vertex to have been reconstructed successfully
313   if (strcmp(vtxESD->GetName(), "default")==0)return;
314   // vertex position and uncertainties
315   fhQA[kVtxPosX] [index]->Fill(vtxESD->GetXv());
316   fhQA[kVtxPosY] [index]->Fill(vtxESD->GetYv());
317   fhQA[kVtxPosZ] [index]->Fill(vtxESD->GetZv());
318   fhQA[kVtxResX] [index]->Fill(vtxESD->GetXRes());
319   fhQA[kVtxResY] [index]->Fill(vtxESD->GetYRes());
320   fhQA[kVtxResZ] [index]->Fill(vtxESD->GetZRes());
321   fhQA[kVtxNCtrb][index]->Fill(vtxESD->GetNContributors());
322   
323 }
324
325 //____________________________________________________________________
326 void AliCFEventRecCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
327 {
328   //
329   //setting x-axis bin limits of QA histogram fhQA[index] 
330   // 
331
332   for(Int_t i=0;i<kNStepQA;i++){
333     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
334     return;
335     }
336     fhQA[index][i]->GetXaxis()->Set(nbins,bins);
337   }
338 }
339 //____________________________________________________________________
340 void AliCFEventRecCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
341 {
342   //
343   //setting x-axis bins and range of QA histogram fhQA[index] 
344   // 
345
346   for(Int_t i=0;i<kNStepQA;i++){
347     if(!fhQA[index][i]){AliWarning("non-existing histogram!");
348     return;
349     }
350     fhQA[index][i]->GetXaxis()->Set(nbins,xmin,xmax);
351   }
352 }
353
354 //_____________________________________________________________________________
355  void AliCFEventRecCuts::DefineHistograms() {
356   //
357   // histograms for cut variables
358   //
359   Int_t color = 2;
360
361   if(!fIsQAOn) {
362     AliInfo(Form("No QA histos requested, Please first set the QA flag on!"));
363     return;
364   }  
365   
366   // book QA histograms
367
368   Char_t str[256];
369   for (Int_t i=0; i<kNStepQA; i++) {
370     if (i==0) sprintf(str," ");
371     else sprintf(str,"_cut");
372
373     fhQA[kNTracks][i]   = new  TH1F(Form("%s_NTracks%s",GetName(),str),                 "",501,-0.5,500.5);
374     fhQA[kVtxPosX][i]   = new  TH1F(Form("%s_Vtx_Pos_X%s",GetName(),str),               "",100,-5.,5.);
375     fhQA[kVtxPosY][i]   = new  TH1F(Form("%s_Vtx_Pos_Y%s",GetName(),str),               "",100,-5.,5.);
376     fhQA[kVtxPosZ][i]   = new  TH1F(Form("%s_Vtx_Pos_Z%s",GetName(),str),               "",200,-50.,50.);
377     fhQA[kVtxResX][i]   = new  TH1F(Form("%s_Vtx_Res_X%s",GetName(),str),               "",100,-1.,1.);
378     fhQA[kVtxResY][i]   = new  TH1F(Form("%s_Vtx_Res_Y%s",GetName(),str),               "",100,-1.,1.);
379     fhQA[kVtxResZ][i]   = new  TH1F(Form("%s_Vtx_Res_Z%s",GetName(),str),               "",100,-1.,1.);
380     fhQA[kVtxNCtrb][i]  = new  TH1F(Form("%s_Vtx_N_Ctrb%s",GetName(),str),              "",1000,0.,1000.);
381  
382     fhQA[kNTracks][i]   ->SetXTitle("Number of ESD tracks");
383     fhQA[kVtxPosX][i]   ->SetXTitle("Vertex Position X (cm)");
384     fhQA[kVtxPosY][i]   ->SetXTitle("Vertex Position Y (cm)");
385     fhQA[kVtxPosZ][i]   ->SetXTitle("Vertex Position Z (cm)");
386     fhQA[kVtxResX][i]   ->SetXTitle("Vertex Resolution X (cm)");
387     fhQA[kVtxResY][i]   ->SetXTitle("Vertex Resolution Y (cm)");
388     fhQA[kVtxResZ][i]   ->SetXTitle("Vertex Resolution Z (cm)");
389     fhQA[kVtxNCtrb][i]  ->SetXTitle("Number of contributors");
390   }
391
392   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
393
394 }
395
396 //_____________________________________________________________________________
397 void AliCFEventRecCuts::AddQAHistograms(TList *qaList) {
398   //
399   // saves the histograms in a TList
400   //
401
402   DefineHistograms();
403
404   for (Int_t j=0; j<kNStepQA; j++) {
405     for(Int_t i=0; i<kNCuts; i++)
406         qaList->Add(fhQA[i][j]);
407   }
408 }