Updates for proper treatment of fQuality flag in centrality selection and for usage...
[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 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////
19 //
20 // Base class for cuts on AOD reconstructed heavy-flavour decay
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24 #include <Riostream.h>
25
26 #include "AliVEvent.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliVVertex.h"
30 #include "AliESDVertex.h"
31 #include "AliLog.h"
32 #include "AliAODVertex.h"
33 #include "AliESDtrack.h"
34 #include "AliAODTrack.h"
35 #include "AliESDtrackCuts.h"
36 #include "AliCentrality.h"
37 #include "AliAODRecoDecayHF.h"
38 #include "AliAnalysisVertexingHF.h"
39 #include "AliRDHFCuts.h"
40
41 ClassImp(AliRDHFCuts)
42
43 //--------------------------------------------------------------------------
44 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : 
45 AliAnalysisCuts(name,title),
46 fMinVtxType(3),
47 fMinVtxContr(1),
48 fMaxVtxRedChi2(1e6),
49 fMaxVtxZ(1e6),
50 fMinSPDMultiplicity(0),
51 fTriggerMask(0),
52 fTriggerClass("CINT1"),
53 fTrackCuts(0),
54 fnPtBins(1),
55 fnPtBinLimits(1),
56 fPtBinLimits(0),
57 fnVars(1),
58 fVarNames(0),
59 fnVarsForOpt(0),
60 fVarsForOpt(0),
61 fGlobalIndex(1),
62 fCutsRD(0),
63 fIsUpperCut(0),
64 fUsePID(kFALSE),
65 fPidHF(0),
66 fWhyRejection(0),
67 fRemoveDaughtersFromPrimary(kFALSE),
68 fOptPileup(0),
69 fMinContrPileup(3),
70 fMinDzPileup(0.6),
71 fUseCentrality(0),
72 fMinCentrality(0.),
73 fMaxCentrality(100.),
74 fFixRefs(kFALSE),
75 fIsSelectedCuts(0),
76 fIsSelectedPID(0),
77 fMinPtCand(-1.),
78 fMaxPtCand(100000.)
79 {
80   //
81   // Default Constructor
82   //
83 }
84 //--------------------------------------------------------------------------
85 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
86   AliAnalysisCuts(source),
87   fMinVtxType(source.fMinVtxType),
88   fMinVtxContr(source.fMinVtxContr),
89   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
90   fMaxVtxZ(source.fMaxVtxZ),
91   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
92   fTriggerMask(source.fTriggerMask),
93   fTriggerClass(source.fTriggerClass),
94   fTrackCuts(0),
95   fnPtBins(source.fnPtBins),
96   fnPtBinLimits(source.fnPtBinLimits),
97   fPtBinLimits(0),
98   fnVars(source.fnVars),
99   fVarNames(0),
100   fnVarsForOpt(source.fnVarsForOpt),
101   fVarsForOpt(0),
102   fGlobalIndex(source.fGlobalIndex),
103   fCutsRD(0),
104   fIsUpperCut(0),
105   fUsePID(source.fUsePID),
106   fPidHF(0),
107   fWhyRejection(source.fWhyRejection),
108   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
109   fOptPileup(source.fOptPileup),
110   fMinContrPileup(source.fMinContrPileup),
111   fMinDzPileup(source.fMinDzPileup),
112   fUseCentrality(source.fUseCentrality),
113   fMinCentrality(source.fMinCentrality),
114   fMaxCentrality(source.fMaxCentrality),
115   fFixRefs(source.fFixRefs),
116   fIsSelectedCuts(source.fIsSelectedCuts),
117   fIsSelectedPID(source.fIsSelectedPID),
118   fMinPtCand(source.fMinPtCand),
119   fMaxPtCand(source.fMaxPtCand)
120 {
121   //
122   // Copy constructor
123   //
124   cout<<"Copy constructor"<<endl;
125   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
126   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
127   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
128   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
129   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
130   if(source.fPidHF) SetPidHF(source.fPidHF);
131   PrintAll();
132
133 }
134 //--------------------------------------------------------------------------
135 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
136 {
137   //
138   // assignment operator
139   //
140   if(&source == this) return *this;
141
142   AliAnalysisCuts::operator=(source);
143
144   fMinVtxType=source.fMinVtxType;
145   fMinVtxContr=source.fMinVtxContr;
146   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
147   fMaxVtxZ=source.fMaxVtxZ;
148   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
149   fTriggerMask=source.fTriggerMask;
150   fTriggerClass=source.fTriggerClass;
151   fnPtBins=source.fnPtBins;
152   fnVars=source.fnVars;
153   fGlobalIndex=source.fGlobalIndex;
154   fnVarsForOpt=source.fnVarsForOpt;
155   fUsePID=source.fUsePID;
156   SetPidHF(source.GetPidHF());
157   fWhyRejection=source.fWhyRejection;
158   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
159   fOptPileup=source.fOptPileup;
160   fMinContrPileup=source.fMinContrPileup;
161   fMinDzPileup=source.fMinDzPileup;
162   fUseCentrality=source.fUseCentrality;
163   fMinCentrality=source.fMinCentrality;
164   fMaxCentrality=source.fMaxCentrality;
165   fFixRefs=source.fFixRefs;
166   fIsSelectedCuts=source.fIsSelectedCuts;
167   fIsSelectedPID=source.fIsSelectedPID;
168   fMinPtCand=source.fMinPtCand;
169   fMaxPtCand=source.fMaxPtCand;
170
171   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
172   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
173   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
174   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
175   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
176   PrintAll();
177
178   return *this;
179 }
180 //--------------------------------------------------------------------------
181 AliRDHFCuts::~AliRDHFCuts() {
182   //  
183   // Default Destructor
184   //
185   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
186   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
187   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
188   if(fCutsRD) {
189     delete [] fCutsRD;
190     fCutsRD=0;
191   }
192   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
193   if(fPidHF){ 
194     delete fPidHF;
195     fPidHF=0;
196   }
197 }
198 //---------------------------------------------------------------------------
199 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
200   //
201   // Centrality selection
202   //
203   
204   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){    
205     AliWarning("Centrality estimator not valid");    
206     return 3;  
207   }else{    
208     Float_t centvalue=GetCentrality((AliAODEvent*)event);          
209     if (centvalue<-998.){//-999 if no centralityP
210       return 0;    
211     }else{      
212       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
213         return 2;      
214       }    
215     } 
216   }  
217   return 0;
218 }
219 //---------------------------------------------------------------------------
220 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
221   //
222   // Event selection
223   // 
224   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
225
226   fWhyRejection=0;
227
228   // trigger class
229   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
230   // don't do for PbPb 2010 data
231   if(event->GetRunNumber()<136851 || event->GetRunNumber()>139517) {
232     if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {
233       fWhyRejection=5;
234       return kFALSE;
235     }
236   }
237
238   // TEMPORARY FIX FOR REFERENCES
239   // Fix references to daughter tracks
240   if(fFixRefs) {
241     AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
242     fixer->FixReferences((AliAODEvent*)event);
243     delete fixer;
244   }
245   //
246
247
248
249   // multiplicity cuts no implemented yet
250    
251   const AliVVertex *vertex = event->GetPrimaryVertex();
252
253   if(!vertex) return kFALSE;
254
255   TString title=vertex->GetTitle();
256   if(title.Contains("Z") && fMinVtxType>1) return kFALSE; 
257   if(title.Contains("3D") && fMinVtxType>2) return kFALSE; 
258
259   if(vertex->GetNContributors()<fMinVtxContr) return kFALSE; 
260
261   if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) return kFALSE;
262
263   // switch to settings for 1-pad cls in TPC
264   if(fPidHF) {
265     if(event->GetRunNumber()>121693 && event->GetRunNumber()<136851) 
266       fPidHF->SetOnePad(kTRUE);
267     if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) 
268       fPidHF->SetPbPb(kTRUE);
269   }
270
271   if(fOptPileup==kRejectPileupEvent){
272     Int_t cutc=(Int_t)fMinContrPileup;
273     Double_t cutz=(Double_t)fMinDzPileup;
274     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
275       fWhyRejection=1;
276       return kFALSE;
277     }
278   }
279
280   //centrality selection
281   if (fUseCentrality!=kCentOff) {  
282     Int_t rejection=IsEventSelectedInCentrality(event);    
283     if(rejection>1){      
284       fWhyRejection=rejection;      
285       return kFALSE;    
286     }
287   }
288
289   return kTRUE;
290 }
291 //---------------------------------------------------------------------------
292 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
293   //
294   // Daughter tracks selection
295   // 
296   if(!fTrackCuts) return kTRUE;
297  
298   Int_t ndaughters = d->GetNDaughters();
299   AliAODVertex *vAOD = d->GetPrimaryVtx();
300   Double_t pos[3],cov[6];
301   vAOD->GetXYZ(pos);
302   vAOD->GetCovarianceMatrix(cov);
303   const AliESDVertex vESD(pos,cov,100.,100);
304   
305   Bool_t retval=kTRUE;
306   
307   for(Int_t idg=0; idg<ndaughters; idg++) {
308     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
309     if(!dgTrack) {retval = kFALSE; continue;}
310     //printf("charge %d\n",dgTrack->Charge());
311     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
312
313     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
314   }
315   
316   return retval;
317 }
318 //---------------------------------------------------------------------------
319 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
320   //
321   // Convert to ESDtrack, relate to vertex and check cuts
322   //
323   if(!cuts) return kTRUE;
324
325   Bool_t retval=kTRUE;
326
327   // convert to ESD track here
328   AliESDtrack esdTrack(track);
329   // set the TPC cluster info
330   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
331   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
332   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
333   // needed to calculate the impact parameters
334   esdTrack.RelateToVertex(primary,0.,3.); 
335
336   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
337  
338   if(fOptPileup==kRejectTracksFromPileupVertex){
339     // to be implemented
340     // we need either to have here the AOD Event, 
341     // or to have the pileup vertex object
342   }
343   return retval; 
344 }
345 //---------------------------------------------------------------------------
346 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
347   // Set the pt bins
348
349   if(fPtBinLimits) {
350     delete [] fPtBinLimits;
351     fPtBinLimits = NULL;
352     printf("Changing the pt bins\n");
353   }
354
355   if(nPtBinLimits != fnPtBins+1){
356     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
357     SetNPtBins(nPtBinLimits-1);
358   }
359
360   fnPtBinLimits = nPtBinLimits;
361   SetGlobalIndex();
362   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
363   fPtBinLimits = new Float_t[fnPtBinLimits];
364   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
365
366   return;
367 }
368 //---------------------------------------------------------------------------
369 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
370   // Set the variable names
371
372   if(fVarNames) {
373     delete [] fVarNames;
374     fVarNames = NULL;
375     //printf("Changing the variable names\n");
376   }
377   if(nVars!=fnVars){
378     printf("Wrong number of variables: it has to be %d\n",fnVars);
379     return;
380   }
381   //fnVars=nVars;
382   fVarNames = new TString[nVars];
383   fIsUpperCut = new Bool_t[nVars];
384   for(Int_t iv=0; iv<nVars; iv++) {
385     fVarNames[iv] = varNames[iv];
386     fIsUpperCut[iv] = isUpperCut[iv];
387   }
388
389   return;
390 }
391 //---------------------------------------------------------------------------
392 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
393   // Set the variables to be used for cuts optimization
394
395   if(fVarsForOpt) {
396     delete [] fVarsForOpt;
397     fVarsForOpt = NULL;
398     //printf("Changing the variables for cut optimization\n");
399   }
400   
401   if(nVars==0){//!=fnVars) {
402     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
403     return;
404   } 
405   
406   fnVarsForOpt = 0;
407   fVarsForOpt = new Bool_t[fnVars];
408   for(Int_t iv=0; iv<fnVars; iv++) {
409     fVarsForOpt[iv]=forOpt[iv];
410     if(fVarsForOpt[iv]) fnVarsForOpt++;
411   }
412
413   return;
414 }
415
416 //---------------------------------------------------------------------------
417 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
418   //
419   // set centrality estimator  
420   //
421   fUseCentrality=flag;
422   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
423  
424   return;
425 }
426
427
428 //---------------------------------------------------------------------------
429 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
430   //
431   // store the cuts
432   //
433   if(nVars!=fnVars) {
434     printf("Wrong number of variables: it has to be %d\n",fnVars);
435     return;
436   } 
437   if(nPtBins!=fnPtBins) {
438     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
439     return;
440   } 
441
442   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
443   
444
445   for(Int_t iv=0; iv<fnVars; iv++) {
446
447     for(Int_t ib=0; ib<fnPtBins; ib++) {
448
449       //check
450       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
451         cout<<"Overflow, exit..."<<endl;
452         return;
453       }
454
455       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
456
457     }
458   }
459   return;
460 }
461 //---------------------------------------------------------------------------
462 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
463   //
464   // store the cuts
465   //
466   if(glIndex != fGlobalIndex){
467     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
468     return;
469   }
470   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
471
472   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
473     fCutsRD[iGl] = cutsRDGlob[iGl];
474   }
475   return;
476 }
477 //---------------------------------------------------------------------------
478 void AliRDHFCuts::PrintAll() const {
479   //
480   // print all cuts values
481   // 
482
483   printf("Minimum vtx type %d\n",fMinVtxType);
484   printf("Minimum vtx contr %d\n",fMinVtxContr);
485   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
486   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
487   printf("Use PID %d\n",(Int_t)fUsePID);
488   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
489   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
490   if(fOptPileup==1) printf(" -- Reject pileup event");
491   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
492   if(fUseCentrality>0) {
493     TString estimator="";
494     if(fUseCentrality==1) estimator = "V0";
495     if(fUseCentrality==2) estimator = "Tracks";
496     if(fUseCentrality==3) estimator = "Tracklets";
497     if(fUseCentrality==4) estimator = "SPD clusters outer"; 
498     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
499   }
500
501   if(fVarNames){
502     cout<<"Array of variables"<<endl;
503     for(Int_t iv=0;iv<fnVars;iv++){
504       cout<<fVarNames[iv]<<"\t";
505     }
506     cout<<endl;
507   }
508   if(fVarsForOpt){
509     cout<<"Array of optimization"<<endl;
510     for(Int_t iv=0;iv<fnVars;iv++){
511       cout<<fVarsForOpt[iv]<<"\t";
512     }
513     cout<<endl;
514   }
515   if(fIsUpperCut){
516     cout<<"Array of upper/lower cut"<<endl;
517    for(Int_t iv=0;iv<fnVars;iv++){
518      cout<<fIsUpperCut[iv]<<"\t";
519    }
520    cout<<endl;
521   }
522   if(fPtBinLimits){
523     cout<<"Array of ptbin limits"<<endl;
524     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
525       cout<<fPtBinLimits[ib]<<"\t";
526     }
527     cout<<endl;
528   }
529   if(fCutsRD){
530     cout<<"Matrix of cuts"<<endl;
531    for(Int_t iv=0;iv<fnVars;iv++){
532      for(Int_t ib=0;ib<fnPtBins;ib++){
533        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
534      } 
535      cout<<endl;
536    }
537    cout<<endl;
538   }
539   return;
540 }
541 //---------------------------------------------------------------------------
542 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
543   //
544   // get the cuts
545   //
546
547   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
548
549
550   Int_t iv,ib;
551   if(!cutsRD) {
552     //cout<<"Initialization..."<<endl;
553     cutsRD=new Float_t*[fnVars];
554     for(iv=0; iv<fnVars; iv++) {
555       cutsRD[iv] = new Float_t[fnPtBins];
556     }
557   }
558   
559   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
560     GetVarPtIndex(iGlobal,iv,ib);
561     cutsRD[iv][ib] = fCutsRD[iGlobal];
562   }
563
564   return;
565 }
566
567 //---------------------------------------------------------------------------
568 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
569   //
570   // give the global index from variable and pt bin
571   //
572   return iPtBin*fnVars+iVar;
573 }
574
575 //---------------------------------------------------------------------------
576 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
577   //
578   //give the index of the variable and of the pt bin from the global index
579   //
580   iPtBin=(Int_t)iGlob/fnVars;
581   iVar=iGlob%fnVars;
582
583   return;
584 }
585
586 //---------------------------------------------------------------------------
587 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
588   //
589   //give the pt bin where the pt lies.
590   //
591   Int_t ptbin=-1;
592   if(pt<fPtBinLimits[0])return ptbin;
593   for (Int_t i=0;i<fnPtBins;i++){
594     if(pt<fPtBinLimits[i+1]) {
595       ptbin=i;
596       break;
597     }
598   }
599   return ptbin;
600 }
601 //-------------------------------------------------------------------
602 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
603   // 
604   // Give the value of cut set for the variable iVar and the pt bin iPtBin
605   //
606   if(!fCutsRD){
607     cout<<"Cuts not iniziaisez yet"<<endl;
608     return 0;
609   }
610   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
611 }
612 //-------------------------------------------------------------------
613 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) const {
614   //
615   // Get centrality percentile
616   //
617   AliAODHeader *header=aodEvent->GetHeader();
618   AliCentrality *centrality=header->GetCentralityP();
619   Float_t cent=-999.;
620   Bool_t isSelRun=kFALSE;
621   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
622   if(!centrality) return cent;
623   else{
624     if (estimator==kCentV0M){
625       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
626       if(cent<0){
627         Int_t quality = centrality->GetQuality();
628         if(quality<=1){
629           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
630         }else{
631           Int_t runnum=aodEvent->GetRunNumber();
632           for(Int_t ir=0;ir<5;ir++){
633             if(runnum==selRun[ir]){
634               isSelRun=kTRUE;
635               break;
636             }
637           }
638           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
639         }
640       }
641     }
642     else {
643       if (estimator==kCentTRK) {
644         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
645         if(cent<0){
646           Int_t quality = centrality->GetQuality();
647           if(quality<=1){
648             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
649           }else{
650             Int_t runnum=aodEvent->GetRunNumber();
651             for(Int_t ir=0;ir<5;ir++){
652               if(runnum==selRun[ir]){
653                 isSelRun=kTRUE;
654                 break;
655               }
656             }
657             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
658           }
659         }
660       }
661       else{
662         if (estimator==kCentTKL){
663           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
664           if(cent<0){
665             Int_t quality = centrality->GetQuality();
666             if(quality<=1){
667               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
668             }else{
669               Int_t runnum=aodEvent->GetRunNumber();
670               for(Int_t ir=0;ir<5;ir++){
671                 if(runnum==selRun[ir]){
672                   isSelRun=kTRUE;
673                   break;
674             }
675               }
676               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
677             }   
678           }
679         }
680         else{
681           if (estimator==kCentCL1){
682             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
683             if(cent<0){
684               Int_t quality = centrality->GetQuality();
685               if(quality<=1){
686                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
687               }else{
688                 Int_t runnum=aodEvent->GetRunNumber();
689                 for(Int_t ir=0;ir<5;ir++){
690                   if(runnum==selRun[ir]){
691                     isSelRun=kTRUE;
692                     break;
693                   }
694                 }
695                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
696               }
697             }
698           }
699           else {
700             AliWarning("Centrality estimator not valid");
701             
702           }
703         }
704       }
705     } 
706   }
707   return cent;
708 }
709 //-------------------------------------------------------------------
710 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
711   //
712   // Compare two cuts objects
713   //
714
715   Bool_t areEqual=kTRUE;
716
717   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
718
719   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
720
721   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
722
723   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
724
725   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
726
727   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
728   if(fTrackCuts){
729     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
730
731     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
732
733     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
734
735     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;}
736   }
737
738   if(fCutsRD) {
739    for(Int_t iv=0;iv<fnVars;iv++) {
740      for(Int_t ib=0;ib<fnPtBins;ib++) {
741        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
742          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
743          areEqual=kFALSE;
744        }
745      }
746    }
747   }
748
749   return areEqual;
750 }
751 //---------------------------------------------------------------------------
752 void AliRDHFCuts::MakeTable() const {
753   //
754   // print cuts values in table format
755   // 
756
757         TString ptString = "pT range";
758         if(fVarNames && fPtBinLimits && fCutsRD){
759                 TString firstLine(Form("*       %-15s",ptString.Data()));
760                 for (Int_t ivar=0; ivar<fnVars; ivar++){
761                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
762                         if (ivar == fnVars){
763                                 firstLine+="*\n";
764                         }
765                 }
766                 Printf("%s",firstLine.Data());
767                 
768                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
769                         TString line;
770                         if (ipt==fnPtBins-1){
771                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
772                         }
773                         else{
774                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
775                         }
776                         for (Int_t ivar=0; ivar<fnVars; ivar++){
777                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
778                         }
779                         Printf("%s",line.Data());
780                 }
781
782         }
783
784
785   return;
786 }
787 //--------------------------------------------------------------------------
788 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod,
789                                       AliAODVertex *origownvtx,AliAODVertex *recvtx) const
790 {
791   //
792   // Recalculate primary vertex without daughters
793   //
794
795   if(!aod) {
796     AliError("Can not remove daughters from vertex without AOD event");
797     return kFALSE;
798   }   
799   if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
800   recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
801   if(!recvtx){
802     AliDebug(2,"Removal of daughter tracks failed");
803     if(origownvtx){
804       delete origownvtx;
805       origownvtx=NULL;
806     }
807     return kFALSE;
808   }
809   //set recalculed primary vertex
810   d->SetOwnPrimaryVtx(recvtx);
811   delete recvtx; recvtx=NULL;
812
813   return kTRUE;
814 }
815 //--------------------------------------------------------------------------
816 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,AliAODVertex *origownvtx) const
817 {
818   //
819   // Clean-up own primary vertex if needed
820   //
821
822   if(origownvtx) {
823     d->SetOwnPrimaryVtx(origownvtx);
824     delete origownvtx;
825     origownvtx=NULL;
826   } else if(fRemoveDaughtersFromPrimary) {
827     d->UnsetOwnPrimaryVtx();
828     AliDebug(3,"delete new vertex\n");
829   }
830
831   return;
832 }
833
834
835