Added method to compare two cuts objects
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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
16 /////////////////////////////////////////////////////////////
17 //
18 // Base class for cuts on AOD reconstructed heavy-flavour decay
19 //
20 // Author: A.Dainese, andrea.dainese@pd.infn.it
21 /////////////////////////////////////////////////////////////
22 #include <Riostream.h>
23
24 #include "AliVEvent.h"
25 #include "AliESDEvent.h"
26 #include "AliAODEvent.h"
27 #include "AliVVertex.h"
28 #include "AliESDVertex.h"
29 #include "AliAODVertex.h"
30 #include "AliESDtrack.h"
31 #include "AliAODTrack.h"
32 #include "AliESDtrackCuts.h"
33 #include "AliAODRecoDecayHF.h"
34 #include "AliRDHFCuts.h"
35
36 ClassImp(AliRDHFCuts)
37
38 //--------------------------------------------------------------------------
39 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : 
40 AliAnalysisCuts(name,title),
41 fMinVtxType(3),
42 fMinVtxContr(1),
43 fMaxVtxRedChi2(1e6),
44 fMinSPDMultiplicity(0),
45 fTriggerMask(0),
46 fTrackCuts(0),
47 fnPtBins(1),
48 fnPtBinLimits(1),
49 fPtBinLimits(0),
50 fnVars(1),
51 fVarNames(0),
52 fnVarsForOpt(0),
53 fVarsForOpt(0),
54 fGlobalIndex(1),
55 fCutsRD(0),
56 fIsUpperCut(0),
57 fUsePID(kFALSE),
58 fPidHF(0),
59 fWhyRejection(0),
60 fRemoveDaughtersFromPrimary(kFALSE)
61 {
62   //
63   // Default Constructor
64   //
65 }
66 //--------------------------------------------------------------------------
67 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
68   AliAnalysisCuts(source),
69   fMinVtxType(source.fMinVtxType),
70   fMinVtxContr(source.fMinVtxContr),
71   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
72   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
73   fTriggerMask(source.fTriggerMask),
74   fTrackCuts(0),
75   fnPtBins(source.fnPtBins),
76   fnPtBinLimits(source.fnPtBinLimits),
77   fPtBinLimits(0),
78   fnVars(source.fnVars),
79   fVarNames(0),
80   fnVarsForOpt(source.fnVarsForOpt),
81   fVarsForOpt(0),
82   fGlobalIndex(source.fGlobalIndex),
83   fCutsRD(0),
84   fIsUpperCut(0),
85   fUsePID(source.fUsePID),
86   fPidHF(0),
87   fWhyRejection(source.fWhyRejection),
88   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary)
89 {
90   //
91   // Copy constructor
92   //
93   cout<<"Copy constructor"<<endl;
94   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
95   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
96   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
97   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
98   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
99   if(source.fPidHF) SetPidHF(source.fPidHF);
100   PrintAll();
101
102 }
103 //--------------------------------------------------------------------------
104 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
105 {
106   //
107   // assignment operator
108   //
109   if(&source == this) return *this;
110
111   AliAnalysisCuts::operator=(source);
112
113   fMinVtxType=source.fMinVtxType;
114   fMinVtxContr=source.fMinVtxContr;
115   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
116   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
117   fTriggerMask=source.fTriggerMask;
118   fnPtBins=source.fnPtBins;
119   fnVars=source.fnVars;
120   fGlobalIndex=source.fGlobalIndex;
121   fnVarsForOpt=source.fnVarsForOpt;
122   fUsePID=source.fUsePID;
123   SetPidHF(source.GetPidHF());
124   fWhyRejection=source.fWhyRejection;
125   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
126
127   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
128   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
129   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
130   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
131   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
132   PrintAll();
133
134   return *this;
135 }
136 //--------------------------------------------------------------------------
137 AliRDHFCuts::~AliRDHFCuts() {
138   //  
139   // Default Destructor
140   //
141   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
142   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
143   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
144   if(fCutsRD) {
145     delete [] fCutsRD;
146     fCutsRD=0;
147   }
148   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
149   if(fPidHF){ 
150     delete fPidHF;
151     fPidHF=0;
152   }
153 }
154 //---------------------------------------------------------------------------
155 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) const {
156   //
157   // Event selection
158   // 
159   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
160
161   // multiplicity cuts no implemented yet
162
163   const AliVVertex *vertex = event->GetPrimaryVertex();
164
165   if(!vertex) return kFALSE;
166
167   TString title=vertex->GetTitle();
168   if(title.Contains("Z") && fMinVtxType>1) return kFALSE; 
169   if(title.Contains("3D") && fMinVtxType>2) return kFALSE; 
170
171   if(vertex->GetNContributors()<fMinVtxContr) return kFALSE; 
172
173   return kTRUE;
174 }
175 //---------------------------------------------------------------------------
176 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
177   //
178   // Daughter tracks selection
179   // 
180   if(!fTrackCuts) return kTRUE;
181
182   Int_t ndaughters = d->GetNDaughters();
183   AliAODVertex *vAOD = d->GetPrimaryVtx();
184   Double_t pos[3],cov[6];
185   vAOD->GetXYZ(pos);
186   vAOD->GetCovarianceMatrix(cov);
187   const AliESDVertex vESD(pos,cov,100.,100);
188
189   Bool_t retval=kTRUE;
190
191   for(Int_t idg=0; idg<ndaughters; idg++) {
192     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
193     if(!dgTrack) retval = kFALSE;
194     //printf("charge %d\n",dgTrack->Charge());
195     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
196
197     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
198   }
199
200   return retval;
201 }
202 //---------------------------------------------------------------------------
203 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
204   //
205   // Convert to ESDtrack, relate to vertex and check cuts
206   //
207   if(!cuts) return kTRUE;
208
209   Bool_t retval=kTRUE;
210
211   // convert to ESD track here
212   AliESDtrack esdTrack(track);
213   // needed to calculate the impact parameters
214   esdTrack.RelateToVertex(primary,0.,3.); 
215   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
216  
217   return retval; 
218 }
219 //---------------------------------------------------------------------------
220 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
221   // Set the pt bins
222
223   if(fPtBinLimits) {
224     delete [] fPtBinLimits;
225     fPtBinLimits = NULL;
226     printf("Changing the pt bins\n");
227   }
228
229   if(nPtBinLimits != fnPtBins+1){
230     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
231     SetNPtBins(nPtBinLimits-1);
232   }
233
234   fnPtBinLimits = nPtBinLimits;
235   SetGlobalIndex();
236   cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
237   fPtBinLimits = new Float_t[fnPtBinLimits];
238   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
239
240   return;
241 }
242 //---------------------------------------------------------------------------
243 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
244   // Set the variable names
245
246   if(fVarNames) {
247     delete [] fVarNames;
248     fVarNames = NULL;
249     printf("Changing the variable names\n");
250   }
251   if(nVars!=fnVars){
252     printf("Wrong number of variables: it has to be %d\n",fnVars);
253     return;
254   }
255   //fnVars=nVars;
256   fVarNames = new TString[nVars];
257   fIsUpperCut = new Bool_t[nVars];
258   for(Int_t iv=0; iv<nVars; iv++) {
259     fVarNames[iv] = varNames[iv];
260     fIsUpperCut[iv] = isUpperCut[iv];
261   }
262
263   return;
264 }
265 //---------------------------------------------------------------------------
266 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
267   // Set the variables to be used for cuts optimization
268
269   if(fVarsForOpt) {
270     delete [] fVarsForOpt;
271     fVarsForOpt = NULL;
272     printf("Changing the variables for cut optimization\n");
273   }
274   
275   if(nVars==0){//!=fnVars) {
276     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
277     return;
278   } 
279   
280   fnVarsForOpt = 0;
281   fVarsForOpt = new Bool_t[fnVars];
282   for(Int_t iv=0; iv<fnVars; iv++) {
283     fVarsForOpt[iv]=forOpt[iv];
284     if(fVarsForOpt[iv]) fnVarsForOpt++;
285   }
286
287   return;
288 }
289 //---------------------------------------------------------------------------
290 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
291   //
292   // store the cuts
293   //
294   if(nVars!=fnVars) {
295     printf("Wrong number of variables: it has to be %d\n",fnVars);
296     return;
297   } 
298   if(nPtBins!=fnPtBins) {
299     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
300     return;
301   } 
302
303   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
304   
305
306   for(Int_t iv=0; iv<fnVars; iv++) {
307
308     for(Int_t ib=0; ib<fnPtBins; ib++) {
309
310       //check
311       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
312         cout<<"Overflow, exit..."<<endl;
313         return;
314       }
315
316       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
317
318     }
319   }
320   return;
321 }
322 //---------------------------------------------------------------------------
323 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
324   //
325   // store the cuts
326   //
327   if(glIndex != fGlobalIndex){
328     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
329     return;
330   }
331   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
332
333   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
334     fCutsRD[iGl] = cutsRDGlob[iGl];
335   }
336   return;
337 }
338 //---------------------------------------------------------------------------
339 void AliRDHFCuts::PrintAll() const {
340   //
341   // print all cuts values
342   // 
343
344   printf("Minimum vtx type %d\n",fMinVtxType);
345   printf("Minimum vtx contr %d\n",fMinVtxContr);
346   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
347   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
348   printf("Use PID %d\n",(Int_t)fUsePID);
349   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
350
351   if(fVarNames){
352     cout<<"Array of variables"<<endl;
353     for(Int_t iv=0;iv<fnVars;iv++){
354       cout<<fVarNames[iv]<<"\t";
355     }
356     cout<<endl;
357   }
358   if(fVarsForOpt){
359     cout<<"Array of optimization"<<endl;
360     for(Int_t iv=0;iv<fnVars;iv++){
361       cout<<fVarsForOpt[iv]<<"\t";
362     }
363     cout<<endl;
364   }
365   if(fIsUpperCut){
366     cout<<"Array of upper/lower cut"<<endl;
367    for(Int_t iv=0;iv<fnVars;iv++){
368      cout<<fIsUpperCut[iv]<<"\t";
369    }
370    cout<<endl;
371   }
372   if(fPtBinLimits){
373     cout<<"Array of ptbin limits"<<endl;
374     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
375       cout<<fPtBinLimits[ib]<<"\t";
376     }
377     cout<<endl;
378   }
379   if(fCutsRD){
380     cout<<"Matrix of cuts"<<endl;
381    for(Int_t iv=0;iv<fnVars;iv++){
382      for(Int_t ib=0;ib<fnPtBins;ib++){
383        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
384      }
385      cout<<endl;
386    }
387    cout<<endl;
388   }
389   return;
390 }
391 //---------------------------------------------------------------------------
392 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
393   //
394   // get the cuts
395   //
396
397   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
398
399
400   Int_t iv,ib;
401   if(!cutsRD) {
402     //cout<<"Initialization..."<<endl;
403     cutsRD=new Float_t*[fnVars];
404     for(iv=0; iv<fnVars; iv++) {
405       cutsRD[iv] = new Float_t[fnPtBins];
406     }
407   }
408   
409   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
410     GetVarPtIndex(iGlobal,iv,ib);
411     cutsRD[iv][ib] = fCutsRD[iGlobal];
412   }
413
414   return;
415 }
416
417 //---------------------------------------------------------------------------
418 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
419   //
420   // give the global index from variable and pt bin
421   //
422   return iPtBin*fnVars+iVar;
423 }
424
425 //---------------------------------------------------------------------------
426 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
427   //
428   //give the index of the variable and of the pt bin from the global index
429   //
430   iPtBin=(Int_t)iGlob/fnVars;
431   iVar=iGlob%fnVars;
432
433   return;
434 }
435
436 //---------------------------------------------------------------------------
437 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
438   //
439   //give the pt bin where the pt lies.
440   //
441   Int_t ptbin=-1;
442   for (Int_t i=0;i<fnPtBins;i++){
443     if(pt<fPtBinLimits[i+1]) {
444       ptbin=i;
445       break;
446     }
447   }
448   return ptbin;
449 }
450 //-------------------------------------------------------------------
451 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
452   // 
453   // Give the value of cut set for the variable iVar and the pt bin iPtBin
454   //
455   if(!fCutsRD){
456     cout<<"Cuts not iniziaisez yet"<<endl;
457     return 0;
458   }
459   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
460 }
461 //-------------------------------------------------------------------
462 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
463   //
464   // Compare two cuts objects
465   //
466
467   Bool_t areEqual=kTRUE;
468
469   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
470
471   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
472
473   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
474
475   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
476
477   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
478
479   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
480
481   if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
482
483   if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
484
485   if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
486
487   if(fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)!=obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)) {printf("ClusterReq SPD %d  %d\n",fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD),obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)); areEqual=kFALSE;}
488
489   if(fCutsRD) {
490    for(Int_t iv=0;iv<fnVars;iv++) {
491      for(Int_t ib=0;ib<fnPtBins;ib++) {
492        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
493          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
494          areEqual=kFALSE;
495        }
496      }
497    }
498   }
499
500   return areEqual;
501 }