Coverity
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCuts.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 /* $Id$ */\r
17 \r
18 /////////////////////////////////////////////////////////////\r
19 //\r
20 // Base class for cuts on AOD reconstructed heavy-flavour decay\r
21 //\r
22 // Author: A.Dainese, andrea.dainese@pd.infn.it\r
23 /////////////////////////////////////////////////////////////\r
24 #include <Riostream.h>\r
25 \r
26 #include "AliVEvent.h"\r
27 #include "AliESDEvent.h"\r
28 #include "AliAODEvent.h"\r
29 #include "AliVVertex.h"\r
30 #include "AliESDVertex.h"\r
31 #include "AliLog.h"\r
32 #include "AliAODVertex.h"\r
33 #include "AliESDtrack.h"\r
34 #include "AliAODTrack.h"\r
35 #include "AliESDtrackCuts.h"\r
36 #include "AliCentrality.h"\r
37 #include "AliAODRecoDecayHF.h"\r
38 #include "AliAnalysisVertexingHF.h"\r
39 #include "AliAODMCHeader.h"\r
40 #include "AliAODMCParticle.h"\r
41 #include "AliRDHFCuts.h"\r
42 #include "AliAnalysisManager.h"\r
43 #include "AliInputEventHandler.h"\r
44 #include "AliPIDResponse.h"\r
45 \r
46 ClassImp(AliRDHFCuts)\r
47 \r
48 //--------------------------------------------------------------------------\r
49 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : \r
50 AliAnalysisCuts(name,title),\r
51 fMinVtxType(3),\r
52 fMinVtxContr(1),\r
53 fMaxVtxRedChi2(1e6),\r
54 fMaxVtxZ(10.),\r
55 fMinSPDMultiplicity(0),\r
56 fTriggerMask(0),\r
57 fTriggerClass("CINT1"),\r
58 fTrackCuts(0),\r
59 fnPtBins(1),\r
60 fnPtBinLimits(1),\r
61 fPtBinLimits(0),\r
62 fnVars(1),\r
63 fVarNames(0),\r
64 fnVarsForOpt(0),\r
65 fVarsForOpt(0),\r
66 fGlobalIndex(1),\r
67 fCutsRD(0),\r
68 fIsUpperCut(0),\r
69 fUsePID(kFALSE),\r
70 fUseAOD049(kFALSE),\r
71 fPidHF(0),\r
72 fWhyRejection(0),\r
73 fEvRejectionBits(0),\r
74 fRemoveDaughtersFromPrimary(kFALSE),\r
75 fUseMCVertex(kFALSE),\r
76 fUsePhysicsSelection(kTRUE),\r
77 fOptPileup(0),\r
78 fMinContrPileup(3),\r
79 fMinDzPileup(0.6),\r
80 fUseCentrality(0),\r
81 fMinCentrality(0.),\r
82 fMaxCentrality(100.),\r
83 fFixRefs(kFALSE),\r
84 fIsSelectedCuts(0),\r
85 fIsSelectedPID(0),\r
86 fMinPtCand(-1.),\r
87 fMaxPtCand(100000.),\r
88 fKeepSignalMC(kFALSE),\r
89 fIsCandTrackSPDFirst(kFALSE),\r
90 fMaxPtCandTrackSPDFirst(0.)\r
91 {\r
92   //\r
93   // Default Constructor\r
94   //\r
95 }\r
96 //--------------------------------------------------------------------------\r
97 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :\r
98   AliAnalysisCuts(source),\r
99   fMinVtxType(source.fMinVtxType),\r
100   fMinVtxContr(source.fMinVtxContr),\r
101   fMaxVtxRedChi2(source.fMaxVtxRedChi2),\r
102   fMaxVtxZ(source.fMaxVtxZ),\r
103   fMinSPDMultiplicity(source.fMinSPDMultiplicity),\r
104   fTriggerMask(source.fTriggerMask),\r
105   fTriggerClass(source.fTriggerClass),\r
106   fTrackCuts(0),\r
107   fnPtBins(source.fnPtBins),\r
108   fnPtBinLimits(source.fnPtBinLimits),\r
109   fPtBinLimits(0),\r
110   fnVars(source.fnVars),\r
111   fVarNames(0),\r
112   fnVarsForOpt(source.fnVarsForOpt),\r
113   fVarsForOpt(0),\r
114   fGlobalIndex(source.fGlobalIndex),\r
115   fCutsRD(0),\r
116   fIsUpperCut(0),\r
117   fUsePID(source.fUsePID),\r
118   fUseAOD049(source.fUseAOD049),\r
119   fPidHF(0),\r
120   fWhyRejection(source.fWhyRejection),\r
121   fEvRejectionBits(source.fEvRejectionBits),\r
122   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),\r
123   fUseMCVertex(source.fUseMCVertex),\r
124   fUsePhysicsSelection(source.fUsePhysicsSelection),\r
125   fOptPileup(source.fOptPileup),\r
126   fMinContrPileup(source.fMinContrPileup),\r
127   fMinDzPileup(source.fMinDzPileup),\r
128   fUseCentrality(source.fUseCentrality),\r
129   fMinCentrality(source.fMinCentrality),\r
130   fMaxCentrality(source.fMaxCentrality),\r
131   fFixRefs(source.fFixRefs),\r
132   fIsSelectedCuts(source.fIsSelectedCuts),\r
133   fIsSelectedPID(source.fIsSelectedPID),\r
134   fMinPtCand(source.fMinPtCand),\r
135   fMaxPtCand(source.fMaxPtCand),\r
136   fKeepSignalMC(source.fKeepSignalMC),\r
137   fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),\r
138   fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst)\r
139 {\r
140   //\r
141   // Copy constructor\r
142   //\r
143   cout<<"Copy constructor"<<endl;\r
144   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());\r
145   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);\r
146   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);\r
147   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);\r
148   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);\r
149   if(source.fPidHF) SetPidHF(source.fPidHF);\r
150   PrintAll();\r
151 \r
152 }\r
153 //--------------------------------------------------------------------------\r
154 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)\r
155 {\r
156   //\r
157   // assignment operator\r
158   //\r
159   if(&source == this) return *this;\r
160 \r
161   AliAnalysisCuts::operator=(source);\r
162 \r
163   fMinVtxType=source.fMinVtxType;\r
164   fMinVtxContr=source.fMinVtxContr;\r
165   fMaxVtxRedChi2=source.fMaxVtxRedChi2;\r
166   fMaxVtxZ=source.fMaxVtxZ;\r
167   fMinSPDMultiplicity=source.fMinSPDMultiplicity;\r
168   fTriggerMask=source.fTriggerMask;\r
169   fTriggerClass=source.fTriggerClass;\r
170   fnPtBins=source.fnPtBins;\r
171   fnVars=source.fnVars;\r
172   fGlobalIndex=source.fGlobalIndex;\r
173   fnVarsForOpt=source.fnVarsForOpt;\r
174   fUsePID=source.fUsePID;\r
175   fUseAOD049=source.fUseAOD049;\r
176   SetPidHF(source.GetPidHF());\r
177   fWhyRejection=source.fWhyRejection;\r
178   fEvRejectionBits=source.fEvRejectionBits;\r
179   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;\r
180   fUseMCVertex=source.fUseMCVertex;\r
181   fUsePhysicsSelection=source.fUsePhysicsSelection;\r
182   fOptPileup=source.fOptPileup;\r
183   fMinContrPileup=source.fMinContrPileup;\r
184   fMinDzPileup=source.fMinDzPileup;\r
185   fUseCentrality=source.fUseCentrality;\r
186   fMinCentrality=source.fMinCentrality;\r
187   fMaxCentrality=source.fMaxCentrality;\r
188   fFixRefs=source.fFixRefs;\r
189   fIsSelectedCuts=source.fIsSelectedCuts;\r
190   fIsSelectedPID=source.fIsSelectedPID;\r
191   fMinPtCand=source.fMinPtCand;\r
192   fMaxPtCand=source.fMaxPtCand;\r
193   fKeepSignalMC=source.fKeepSignalMC;\r
194   fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;\r
195   fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;\r
196 \r
197   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());\r
198   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);\r
199   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);\r
200   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);\r
201   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);\r
202   PrintAll();\r
203 \r
204   return *this;\r
205 }\r
206 //--------------------------------------------------------------------------\r
207 AliRDHFCuts::~AliRDHFCuts() {\r
208   //  \r
209   // Default Destructor\r
210   //\r
211   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}\r
212   if(fVarNames) {delete [] fVarNames; fVarNames=0;}\r
213   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}\r
214   if(fCutsRD) {\r
215     delete [] fCutsRD;\r
216     fCutsRD=0;\r
217   }\r
218   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}\r
219   if(fPidHF){ \r
220     delete fPidHF;\r
221     fPidHF=0;\r
222   }\r
223 }\r
224 //---------------------------------------------------------------------------\r
225 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {\r
226   //\r
227   // Centrality selection\r
228   //\r
229   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){    \r
230     AliWarning("Centrality estimator not valid");    \r
231     return 3;  \r
232   }else{    \r
233     Float_t centvalue=GetCentrality((AliAODEvent*)event);          \r
234     if (centvalue<-998.){//-999 if no centralityP\r
235       return 0;    \r
236     }else{      \r
237       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){\r
238         return 2;      \r
239       }    \r
240     } \r
241   }  \r
242   return 0;\r
243 }\r
244 //---------------------------------------------------------------------------\r
245 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {\r
246   //\r
247   // Event selection\r
248   // \r
249   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;\r
250 \r
251   fWhyRejection=0;\r
252   fEvRejectionBits=0;\r
253   Bool_t accept=kTRUE;\r
254 \r
255   // check if it's MC\r
256   Bool_t isMC=kFALSE;\r
257   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
258   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}\r
259 \r
260   // settings for the TPC dE/dx BB parameterization\r
261   if(fPidHF && fPidHF->GetOldPid()) {\r
262     // pp, from LHC10d onwards\r
263     if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||\r
264        event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);\r
265     // pp, 2011 low energy run\r
266     if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){\r
267       fPidHF->SetppLowEn2011(kTRUE);\r
268       fPidHF->SetOnePad(kFALSE);\r
269     }\r
270     // PbPb LHC10h\r
271     if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);\r
272     // MC\r
273     if(isMC) fPidHF->SetMC(kTRUE);\r
274     if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))\r
275       fPidHF->SetMClowenpp2011(kTRUE);\r
276   } \r
277   else if(fPidHF && !fPidHF->GetOldPid()) {\r
278    if(fPidHF->GetPidResponse()==0x0){\r
279     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
280     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();\r
281     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();\r
282     fPidHF->SetPidResponse(pidResp);\r
283    }\r
284   }\r
285 \r
286 \r
287   // trigger class\r
288   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();\r
289   // don't do for MC and for PbPb 2010 data\r
290   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {\r
291     if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {\r
292       fWhyRejection=5;\r
293       fEvRejectionBits+=1<<kNotSelTrigger;\r
294       accept=kFALSE;\r
295     }\r
296   }\r
297 \r
298   // TEMPORARY FIX FOR REFERENCES\r
299   // Fix references to daughter tracks\r
300   //  if(fFixRefs) {\r
301   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();\r
302   //    fixer->FixReferences((AliAODEvent*)event);\r
303   //    delete fixer;\r
304   //  }\r
305   //\r
306 \r
307 \r
308   // physics selection requirements\r
309   if(fUsePhysicsSelection){\r
310     Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
311     if(!isSelected) {\r
312       if(accept) fWhyRejection=7;\r
313       fEvRejectionBits+=1<<kPhysicsSelection;\r
314       accept=kFALSE;\r
315     }\r
316   }\r
317 \r
318   // vertex requirements\r
319    \r
320   const AliVVertex *vertex = event->GetPrimaryVertex();\r
321 \r
322   if(!vertex){\r
323     accept=kFALSE;\r
324     fEvRejectionBits+=1<<kNoVertex;\r
325   }else{\r
326     TString title=vertex->GetTitle();\r
327     if(title.Contains("Z") && fMinVtxType>1){\r
328       accept=kFALSE;\r
329       fEvRejectionBits+=1<<kNoVertex;\r
330     }\r
331     else if(title.Contains("3D") && fMinVtxType>2){\r
332       accept=kFALSE;\r
333       fEvRejectionBits+=1<<kNoVertex;\r
334     }\r
335     if(vertex->GetNContributors()<fMinVtxContr){\r
336       accept=kFALSE;\r
337       fEvRejectionBits+=1<<kTooFewVtxContrib;\r
338     }\r
339     if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {\r
340       fEvRejectionBits+=1<<kZVtxOutFid;\r
341       if(accept) fWhyRejection=6;\r
342       accept=kFALSE;\r
343     } \r
344   }\r
345 \r
346 \r
347   // pile-up rejection\r
348   if(fOptPileup==kRejectPileupEvent){\r
349     Int_t cutc=(Int_t)fMinContrPileup;\r
350     Double_t cutz=(Double_t)fMinDzPileup;\r
351     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {\r
352       if(accept) fWhyRejection=1;\r
353       fEvRejectionBits+=1<<kPileupSPD;\r
354       accept=kFALSE;\r
355     }\r
356   }\r
357 \r
358   // centrality selection\r
359   if (fUseCentrality!=kCentOff) {  \r
360     Int_t rejection=IsEventSelectedInCentrality(event);    \r
361     if(rejection>1){      \r
362       if(accept) fWhyRejection=rejection;      \r
363       fEvRejectionBits+=1<<kOutsideCentrality;\r
364       accept=kFALSE;\r
365     }\r
366   }\r
367 \r
368   return accept;\r
369 }\r
370 //---------------------------------------------------------------------------\r
371 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {\r
372   //\r
373   // Daughter tracks selection\r
374   // \r
375   if(!fTrackCuts) return kTRUE;\r
376  \r
377   Int_t ndaughters = d->GetNDaughters();\r
378   AliAODVertex *vAOD = d->GetPrimaryVtx();\r
379   Double_t pos[3],cov[6];\r
380   vAOD->GetXYZ(pos);\r
381   vAOD->GetCovarianceMatrix(cov);\r
382   const AliESDVertex vESD(pos,cov,100.,100);\r
383   \r
384   Bool_t retval=kTRUE;\r
385   \r
386   for(Int_t idg=0; idg<ndaughters; idg++) {\r
387     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);\r
388     if(!dgTrack) {retval = kFALSE; continue;}\r
389     //printf("charge %d\n",dgTrack->Charge());\r
390     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0\r
391 \r
392     if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)\r
393       { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }\r
394 \r
395     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;\r
396   }\r
397   \r
398   return retval;\r
399 }\r
400 //---------------------------------------------------------------------------\r
401 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {\r
402   //\r
403   // Convert to ESDtrack, relate to vertex and check cuts\r
404   //\r
405   if(!cuts) return kTRUE;\r
406 \r
407   Bool_t retval=kTRUE;\r
408 \r
409   // convert to ESD track here\r
410   AliESDtrack esdTrack(track);\r
411   // set the TPC cluster info\r
412   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());\r
413   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());\r
414   esdTrack.SetTPCPointsF(track->GetTPCNclsF());\r
415   // needed to calculate the impact parameters\r
416   esdTrack.RelateToVertex(primary,0.,3.); \r
417 \r
418   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;\r
419  \r
420   if(fOptPileup==kRejectTracksFromPileupVertex){\r
421     // to be implemented\r
422     // we need either to have here the AOD Event, \r
423     // or to have the pileup vertex object\r
424   }\r
425   return retval; \r
426 }\r
427 //---------------------------------------------------------------------------\r
428 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {\r
429   // Set the pt bins\r
430 \r
431   if(fPtBinLimits) {\r
432     delete [] fPtBinLimits;\r
433     fPtBinLimits = NULL;\r
434     printf("Changing the pt bins\n");\r
435   }\r
436 \r
437   if(nPtBinLimits != fnPtBins+1){\r
438     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;\r
439     SetNPtBins(nPtBinLimits-1);\r
440   }\r
441 \r
442   fnPtBinLimits = nPtBinLimits;\r
443   SetGlobalIndex();\r
444   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;\r
445   fPtBinLimits = new Float_t[fnPtBinLimits];\r
446   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];\r
447 \r
448   return;\r
449 }\r
450 //---------------------------------------------------------------------------\r
451 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){\r
452   // Set the variable names\r
453 \r
454   if(fVarNames) {\r
455     delete [] fVarNames;\r
456     fVarNames = NULL;\r
457     //printf("Changing the variable names\n");\r
458   }\r
459   if(nVars!=fnVars){\r
460     printf("Wrong number of variables: it has to be %d\n",fnVars);\r
461     return;\r
462   }\r
463   //fnVars=nVars;\r
464   fVarNames = new TString[nVars];\r
465   fIsUpperCut = new Bool_t[nVars];\r
466   for(Int_t iv=0; iv<nVars; iv++) {\r
467     fVarNames[iv] = varNames[iv];\r
468     fIsUpperCut[iv] = isUpperCut[iv];\r
469   }\r
470 \r
471   return;\r
472 }\r
473 //---------------------------------------------------------------------------\r
474 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {\r
475   // Set the variables to be used for cuts optimization\r
476 \r
477   if(fVarsForOpt) {\r
478     delete [] fVarsForOpt;\r
479     fVarsForOpt = NULL;\r
480     //printf("Changing the variables for cut optimization\n");\r
481   }\r
482   \r
483   if(nVars==0){//!=fnVars) {\r
484     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);\r
485     return;\r
486   } \r
487   \r
488   fnVarsForOpt = 0;\r
489   fVarsForOpt = new Bool_t[fnVars];\r
490   for(Int_t iv=0; iv<fnVars; iv++) {\r
491     fVarsForOpt[iv]=forOpt[iv];\r
492     if(fVarsForOpt[iv]) fnVarsForOpt++;\r
493   }\r
494 \r
495   return;\r
496 }\r
497 \r
498 //---------------------------------------------------------------------------\r
499 void AliRDHFCuts::SetUseCentrality(Int_t flag) {\r
500   //\r
501   // set centrality estimator  \r
502   //\r
503   fUseCentrality=flag;\r
504   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");\r
505  \r
506   return;\r
507 }\r
508 \r
509 \r
510 //---------------------------------------------------------------------------\r
511 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {\r
512   //\r
513   // store the cuts\r
514   //\r
515   if(nVars!=fnVars) {\r
516     printf("Wrong number of variables: it has to be %d\n",fnVars);\r
517     AliFatal("exiting");\r
518   } \r
519   if(nPtBins!=fnPtBins) {\r
520     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);\r
521     AliFatal("exiting");\r
522   } \r
523 \r
524   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];\r
525   \r
526 \r
527   for(Int_t iv=0; iv<fnVars; iv++) {\r
528 \r
529     for(Int_t ib=0; ib<fnPtBins; ib++) {\r
530 \r
531       //check\r
532       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {\r
533         cout<<"Overflow, exit..."<<endl;\r
534         return;\r
535       }\r
536 \r
537       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];\r
538 \r
539     }\r
540   }\r
541   return;\r
542 }\r
543 //---------------------------------------------------------------------------\r
544 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){\r
545   //\r
546   // store the cuts\r
547   //\r
548   if(glIndex != fGlobalIndex){\r
549     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;\r
550     AliFatal("exiting");\r
551   }\r
552   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];\r
553 \r
554   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){\r
555     fCutsRD[iGl] = cutsRDGlob[iGl];\r
556   }\r
557   return;\r
558 }\r
559 //---------------------------------------------------------------------------\r
560 void AliRDHFCuts::PrintAll() const {\r
561   //\r
562   // print all cuts values\r
563   // \r
564 \r
565   printf("Minimum vtx type %d\n",fMinVtxType);\r
566   printf("Minimum vtx contr %d\n",fMinVtxContr);\r
567   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);\r
568   printf("Min SPD mult %d\n",fMinSPDMultiplicity);\r
569   printf("Use PID %d\n",(Int_t)fUsePID);\r
570   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);\r
571   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");\r
572   if(fOptPileup==1) printf(" -- Reject pileup event");\r
573   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");\r
574   if(fUseCentrality>0) {\r
575     TString estimator="";\r
576     if(fUseCentrality==1) estimator = "V0";\r
577     if(fUseCentrality==2) estimator = "Tracks";\r
578     if(fUseCentrality==3) estimator = "Tracklets";\r
579     if(fUseCentrality==4) estimator = "SPD clusters outer"; \r
580     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());\r
581   }\r
582   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);\r
583 \r
584   if(fVarNames){\r
585     cout<<"Array of variables"<<endl;\r
586     for(Int_t iv=0;iv<fnVars;iv++){\r
587       cout<<fVarNames[iv]<<"\t";\r
588     }\r
589     cout<<endl;\r
590   }\r
591   if(fVarsForOpt){\r
592     cout<<"Array of optimization"<<endl;\r
593     for(Int_t iv=0;iv<fnVars;iv++){\r
594       cout<<fVarsForOpt[iv]<<"\t";\r
595     }\r
596     cout<<endl;\r
597   }\r
598   if(fIsUpperCut){\r
599     cout<<"Array of upper/lower cut"<<endl;\r
600    for(Int_t iv=0;iv<fnVars;iv++){\r
601      cout<<fIsUpperCut[iv]<<"\t";\r
602    }\r
603    cout<<endl;\r
604   }\r
605   if(fPtBinLimits){\r
606     cout<<"Array of ptbin limits"<<endl;\r
607     for(Int_t ib=0;ib<fnPtBinLimits;ib++){\r
608       cout<<fPtBinLimits[ib]<<"\t";\r
609     }\r
610     cout<<endl;\r
611   }\r
612   if(fCutsRD){\r
613     cout<<"Matrix of cuts"<<endl;\r
614    for(Int_t iv=0;iv<fnVars;iv++){\r
615      for(Int_t ib=0;ib<fnPtBins;ib++){\r
616        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";\r
617      } \r
618      cout<<endl;\r
619    }\r
620    cout<<endl;\r
621   }\r
622   return;\r
623 }\r
624 //---------------------------------------------------------------------------\r
625 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{\r
626   //\r
627   // get the cuts\r
628   //\r
629 \r
630   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;\r
631 \r
632 \r
633   Int_t iv,ib;\r
634   if(!cutsRD) {\r
635     //cout<<"Initialization..."<<endl;\r
636     cutsRD=new Float_t*[fnVars];\r
637     for(iv=0; iv<fnVars; iv++) {\r
638       cutsRD[iv] = new Float_t[fnPtBins];\r
639     }\r
640   }\r
641   \r
642   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {\r
643     GetVarPtIndex(iGlobal,iv,ib);\r
644     cutsRD[iv][ib] = fCutsRD[iGlobal];\r
645   }\r
646 \r
647   return;\r
648 }\r
649 \r
650 //---------------------------------------------------------------------------\r
651 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{\r
652   //\r
653   // give the global index from variable and pt bin\r
654   //\r
655   return iPtBin*fnVars+iVar;\r
656 }\r
657 \r
658 //---------------------------------------------------------------------------\r
659 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {\r
660   //\r
661   //give the index of the variable and of the pt bin from the global index\r
662   //\r
663   iPtBin=(Int_t)iGlob/fnVars;\r
664   iVar=iGlob%fnVars;\r
665 \r
666   return;\r
667 }\r
668 \r
669 //---------------------------------------------------------------------------\r
670 Int_t AliRDHFCuts::PtBin(Double_t pt) const {\r
671   //\r
672   //give the pt bin where the pt lies.\r
673   //\r
674   Int_t ptbin=-1;\r
675   if(pt<fPtBinLimits[0])return ptbin;\r
676   for (Int_t i=0;i<fnPtBins;i++){\r
677     if(pt<fPtBinLimits[i+1]) {\r
678       ptbin=i;\r
679       break;\r
680     }\r
681   }\r
682   return ptbin;\r
683 }\r
684 //-------------------------------------------------------------------\r
685 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {\r
686   // \r
687   // Give the value of cut set for the variable iVar and the pt bin iPtBin\r
688   //\r
689   if(!fCutsRD){\r
690     cout<<"Cuts not iniziaisez yet"<<endl;\r
691     return 0;\r
692   }\r
693   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];\r
694 }\r
695 //-------------------------------------------------------------------\r
696 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {\r
697   //\r
698   // Get centrality percentile\r
699   //\r
700 \r
701   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
702   if(mcArray) {fUseAOD049=kFALSE;}\r
703 \r
704   AliAODHeader *header=aodEvent->GetHeader();\r
705   AliCentrality *centrality=header->GetCentralityP();\r
706   Float_t cent=-999.;\r
707   Bool_t isSelRun=kFALSE;\r
708   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};\r
709   if(!centrality) return cent;\r
710   else{\r
711     if (estimator==kCentV0M){\r
712       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));\r
713       if(cent<0){\r
714         Int_t quality = centrality->GetQuality();\r
715         if(quality<=1){\r
716           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
717         }else{\r
718           Int_t runnum=aodEvent->GetRunNumber();\r
719           for(Int_t ir=0;ir<5;ir++){\r
720             if(runnum==selRun[ir]){\r
721               isSelRun=kTRUE;\r
722               break;\r
723             }\r
724           }\r
725           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
726         }\r
727       }\r
728 \r
729       //temporary fix for AOD049 outliers\r
730       if(fUseAOD049&&cent>=0){\r
731         Float_t v0=0;\r
732         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();\r
733         v0+=aodV0->GetMTotV0A();\r
734         v0+=aodV0->GetMTotV0C();\r
735         if(cent==0&&v0<19500)return -1;//filtering issue\r
736         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());\r
737         Float_t val= 1.30552 +  0.147931 * v0;\r
738         Float_t tklSigma[101]={176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86, 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654, 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334, 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224, 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255, 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398, 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235, 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504, 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};\r
739         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier\r
740       }\r
741     }\r
742     else {\r
743       if (estimator==kCentTRK) {\r
744         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));\r
745         if(cent<0){\r
746           Int_t quality = centrality->GetQuality();\r
747           if(quality<=1){\r
748             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
749           }else{\r
750             Int_t runnum=aodEvent->GetRunNumber();\r
751             for(Int_t ir=0;ir<5;ir++){\r
752               if(runnum==selRun[ir]){\r
753                 isSelRun=kTRUE;\r
754                 break;\r
755               }\r
756             }\r
757             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
758           }\r
759         }\r
760       }\r
761       else{\r
762         if (estimator==kCentTKL){\r
763           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));\r
764           if(cent<0){\r
765             Int_t quality = centrality->GetQuality();\r
766             if(quality<=1){\r
767               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
768             }else{\r
769               Int_t runnum=aodEvent->GetRunNumber();\r
770               for(Int_t ir=0;ir<5;ir++){\r
771                 if(runnum==selRun[ir]){\r
772                   isSelRun=kTRUE;\r
773                   break;\r
774             }\r
775               }\r
776               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
777             }   \r
778           }\r
779         }\r
780         else{\r
781           if (estimator==kCentCL1){\r
782             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));\r
783             if(cent<0){\r
784               Int_t quality = centrality->GetQuality();\r
785               if(quality<=1){\r
786                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
787               }else{\r
788                 Int_t runnum=aodEvent->GetRunNumber();\r
789                 for(Int_t ir=0;ir<5;ir++){\r
790                   if(runnum==selRun[ir]){\r
791                     isSelRun=kTRUE;\r
792                     break;\r
793                   }\r
794                 }\r
795                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
796               }\r
797             }\r
798           }\r
799           else {\r
800             AliWarning("Centrality estimator not valid");\r
801             \r
802           }\r
803         }\r
804       }\r
805     } \r
806   }\r
807   return cent;\r
808 }\r
809 //-------------------------------------------------------------------\r
810 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {\r
811   //\r
812   // Compare two cuts objects\r
813   //\r
814 \r
815   Bool_t areEqual=kTRUE;\r
816 \r
817   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}\r
818 \r
819   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}\r
820 \r
821   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}\r
822 \r
823   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}\r
824 \r
825   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}\r
826 \r
827   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}\r
828   if(fTrackCuts){\r
829     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}\r
830 \r
831     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}\r
832 \r
833     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}\r
834 \r
835     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;}\r
836   }\r
837 \r
838   if(fCutsRD) {\r
839    for(Int_t iv=0;iv<fnVars;iv++) {\r
840      for(Int_t ib=0;ib<fnPtBins;ib++) {\r
841        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {\r
842          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";\r
843          areEqual=kFALSE;\r
844        }\r
845      }\r
846    }\r
847   }\r
848 \r
849   return areEqual;\r
850 }\r
851 //---------------------------------------------------------------------------\r
852 void AliRDHFCuts::MakeTable() const {\r
853   //\r
854   // print cuts values in table format\r
855   // \r
856 \r
857         TString ptString = "pT range";\r
858         if(fVarNames && fPtBinLimits && fCutsRD){\r
859                 TString firstLine(Form("*       %-15s",ptString.Data()));\r
860                 for (Int_t ivar=0; ivar<fnVars; ivar++){\r
861                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());\r
862                         if (ivar == fnVars){\r
863                                 firstLine+="*\n";\r
864                         }\r
865                 }\r
866                 Printf("%s",firstLine.Data());\r
867                 \r
868                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){\r
869                         TString line;\r
870                         if (ipt==fnPtBins-1){\r
871                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);\r
872                         }\r
873                         else{\r
874                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);\r
875                         }\r
876                         for (Int_t ivar=0; ivar<fnVars; ivar++){\r
877                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);\r
878                         }\r
879                         Printf("%s",line.Data());\r
880                 }\r
881 \r
882         }\r
883 \r
884 \r
885   return;\r
886 }\r
887 //--------------------------------------------------------------------------\r
888 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
889                                         AliAODEvent *aod) const\r
890 {\r
891   //\r
892   // Recalculate primary vertex without daughters\r
893   //\r
894 \r
895   if(!aod) {\r
896     AliError("Can not remove daughters from vertex without AOD event");\r
897     return 0;\r
898   }   \r
899 \r
900   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);\r
901   if(!recvtx){\r
902     AliDebug(2,"Removal of daughter tracks failed");\r
903     return kFALSE;\r
904   }\r
905 \r
906 \r
907   //set recalculed primary vertex\r
908   d->SetOwnPrimaryVtx(recvtx);\r
909   delete recvtx;\r
910 \r
911   return kTRUE;\r
912 }\r
913 //--------------------------------------------------------------------------\r
914 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const\r
915 {\r
916   //\r
917   // Recalculate primary vertex without daughters\r
918   //\r
919 \r
920   if(!aod) {\r
921     AliError("Can not get MC vertex without AOD event");\r
922     return kFALSE;\r
923   }   \r
924 \r
925   // load MC header\r
926   AliAODMCHeader *mcHeader = \r
927     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
928   if(!mcHeader) {\r
929     AliError("Can not get MC vertex without AODMCHeader event");\r
930     return kFALSE;\r
931   }\r
932   Double_t pos[3];\r
933   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};\r
934   mcHeader->GetVertex(pos);\r
935   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);\r
936 \r
937   if(!recvtx){\r
938     AliDebug(2,"Removal of daughter tracks failed");\r
939     return kFALSE;\r
940   }\r
941 \r
942   //set recalculed primary vertex\r
943   d->SetOwnPrimaryVtx(recvtx);\r
944 \r
945   d->RecalculateImpPars(recvtx,aod);\r
946 \r
947   delete recvtx;\r
948 \r
949   return kTRUE;\r
950 }\r
951 //--------------------------------------------------------------------------\r
952 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
953                                      AliAODEvent *aod,\r
954                                      AliAODVertex *origownvtx) const\r
955 {\r
956   //\r
957   // Clean-up own primary vertex if needed\r
958   //\r
959 \r
960   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {\r
961     d->UnsetOwnPrimaryVtx();\r
962     if(origownvtx) {\r
963       d->SetOwnPrimaryVtx(origownvtx);\r
964       delete origownvtx; origownvtx=NULL;\r
965     }\r
966     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);\r
967   } else {\r
968     if(origownvtx) {\r
969       delete origownvtx; origownvtx=NULL;\r
970     }\r
971   }\r
972   return;\r
973 }\r
974 //--------------------------------------------------------------------------\r
975 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const \r
976 {\r
977   //\r
978   // Checks if this candidate is matched to MC signal\r
979   // \r
980 \r
981   if(!aod) return kFALSE;\r
982 \r
983   // get the MC array\r
984   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
985 \r
986   if(!mcArray) return kFALSE;\r
987 \r
988   // try to match  \r
989   Int_t label = d->MatchToMC(pdg,mcArray);\r
990   \r
991   if(label>=0) {\r
992     //printf("MATCH!\n");\r
993     return kTRUE;\r
994   }\r
995 \r
996   return kFALSE;\r
997 }\r
998 \r
999 \r