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