ee0e8ec6b641adcbee31f1803818b5579d4dcf4f
[u/mrichter/AliRoot.git] / PWGHF / 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 "AliAODMCHeader.h"
40 #include "AliAODMCParticle.h"
41 #include "AliVertexerTracks.h"
42 #include "AliRDHFCuts.h"
43 #include "AliAnalysisManager.h"
44 #include "AliInputEventHandler.h"
45 #include "AliPIDResponse.h"
46 #include "AliAnalysisUtils.h"
47 #include "TRandom.h"
48 #include <TF1.h>
49
50 using std::cout;
51 using std::endl;
52
53 ClassImp(AliRDHFCuts)
54
55 //--------------------------------------------------------------------------
56 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : 
57 AliAnalysisCuts(name,title),
58 fMinVtxType(3),
59 fMinVtxContr(1),
60 fMaxVtxRedChi2(1e6),
61 fMaxVtxZ(10.),
62 fMinSPDMultiplicity(0),
63 fTriggerMask(AliVEvent::kAnyINT),
64 fUseOnlyOneTrigger(kFALSE),
65 fTrackCuts(0),
66 fnPtBins(1),
67 fnPtBinLimits(1),
68 fPtBinLimits(0),
69 fnVars(1),
70 fVarNames(0),
71 fnVarsForOpt(0),
72 fVarsForOpt(0),
73 fGlobalIndex(1),
74 fCutsRD(0),
75 fIsUpperCut(0),
76 fUsePID(kFALSE),
77 fUseAOD049(kFALSE),
78 fPidHF(0),
79 fWhyRejection(0),
80 fEvRejectionBits(0),
81 fRemoveDaughtersFromPrimary(kFALSE),
82 fRecomputePrimVertex(kFALSE),
83 fUseMCVertex(kFALSE),
84 fUsePhysicsSelection(kTRUE),
85 fOptPileup(0),
86 fMinContrPileup(3),
87 fMinDzPileup(0.6),
88 fUseCentrality(0),
89 fMinCentrality(0.),
90 fMaxCentrality(100.),
91 fFixRefs(kFALSE),
92 fIsSelectedCuts(0),
93 fIsSelectedPID(0),
94 fMinPtCand(-1.),
95 fMaxPtCand(100000.),
96 fMaxRapidityCand(-999.),
97 fKeepSignalMC(kFALSE),
98 fIsCandTrackSPDFirst(kFALSE),
99 fMaxPtCandTrackSPDFirst(0.),
100 fApplySPDDeadPbPb2011(kFALSE),
101 fMaxDiffTRKV0Centr(-1.),
102 fRemoveTrackletOutliers(kFALSE),
103 fCutOnzVertexSPD(0),
104 fKinkReject(kFALSE),
105 fUseTrackSelectionWithFilterBits(kTRUE),
106 fUseCentrFlatteningInMC(kFALSE),
107 fHistCentrDistr(0x0),
108 fCutRatioClsOverCrossRowsTPC(0),
109 fCutRatioSignalNOverCrossRowsTPC(0),
110 fCutMinCrossedRowsTPCPtDep(""),
111 f1CutMinNCrossedRowsTPCPtDep(0x0)
112 {
113   //
114   // Default Constructor
115   //
116   fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
117 }
118 //--------------------------------------------------------------------------
119 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
120   AliAnalysisCuts(source),
121   fMinVtxType(source.fMinVtxType),
122   fMinVtxContr(source.fMinVtxContr),
123   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
124   fMaxVtxZ(source.fMaxVtxZ),
125   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
126   fTriggerMask(source.fTriggerMask),
127   fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
128   fTriggerClass(),
129   fTrackCuts(0),
130   fnPtBins(source.fnPtBins),
131   fnPtBinLimits(source.fnPtBinLimits),
132   fPtBinLimits(0),
133   fnVars(source.fnVars),
134   fVarNames(0),
135   fnVarsForOpt(source.fnVarsForOpt),
136   fVarsForOpt(0),
137   fGlobalIndex(source.fGlobalIndex),
138   fCutsRD(0),
139   fIsUpperCut(0),
140   fUsePID(source.fUsePID),
141   fUseAOD049(source.fUseAOD049),
142   fPidHF(0),
143   fWhyRejection(source.fWhyRejection),
144   fEvRejectionBits(source.fEvRejectionBits),
145   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
146   fRecomputePrimVertex(source.fRecomputePrimVertex),
147   fUseMCVertex(source.fUseMCVertex),
148   fUsePhysicsSelection(source.fUsePhysicsSelection),
149   fOptPileup(source.fOptPileup),
150   fMinContrPileup(source.fMinContrPileup),
151   fMinDzPileup(source.fMinDzPileup),
152   fUseCentrality(source.fUseCentrality),
153   fMinCentrality(source.fMinCentrality),
154   fMaxCentrality(source.fMaxCentrality),
155   fFixRefs(source.fFixRefs),
156   fIsSelectedCuts(source.fIsSelectedCuts),
157   fIsSelectedPID(source.fIsSelectedPID),
158   fMinPtCand(source.fMinPtCand),
159   fMaxPtCand(source.fMaxPtCand),
160   fMaxRapidityCand(source.fMaxRapidityCand),
161   fKeepSignalMC(source.fKeepSignalMC),
162   fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
163   fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
164   fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
165   fMaxDiffTRKV0Centr(source.fMaxDiffTRKV0Centr),
166   fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
167   fCutOnzVertexSPD(source.fCutOnzVertexSPD),
168   fKinkReject(source.fKinkReject),
169   fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
170   fUseCentrFlatteningInMC(source.fUseCentrFlatteningInMC),
171   fHistCentrDistr(0x0),
172   fCutRatioClsOverCrossRowsTPC(source.fCutRatioClsOverCrossRowsTPC),
173   fCutRatioSignalNOverCrossRowsTPC(source.fCutRatioSignalNOverCrossRowsTPC),
174   fCutMinCrossedRowsTPCPtDep(""),
175   f1CutMinNCrossedRowsTPCPtDep(0x0)
176 {
177   //
178   // Copy constructor
179   //
180   cout<<"Copy constructor"<<endl;
181   fTriggerClass[0] = source.fTriggerClass[0]; 
182   fTriggerClass[1] = source.fTriggerClass[1];
183   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
184   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
185   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
186   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
187   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
188   if(source.fPidHF) SetPidHF(source.fPidHF);
189   if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
190   if(source.fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
191   if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
192   PrintAll();
193
194 }
195 //--------------------------------------------------------------------------
196 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
197 {
198   //
199   // assignment operator
200   //
201   if(&source == this) return *this;
202
203   AliAnalysisCuts::operator=(source);
204
205   fMinVtxType=source.fMinVtxType;
206   fMinVtxContr=source.fMinVtxContr;
207   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
208   fMaxVtxZ=source.fMaxVtxZ;
209   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
210   fTriggerMask=source.fTriggerMask;
211   fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
212   fTriggerClass[0]=source.fTriggerClass[0];
213   fTriggerClass[1]=source.fTriggerClass[1];
214   fnPtBins=source.fnPtBins;
215   fnPtBinLimits=source.fnPtBinLimits;
216   fnVars=source.fnVars;
217   fGlobalIndex=source.fGlobalIndex;
218   fnVarsForOpt=source.fnVarsForOpt;
219   fUsePID=source.fUsePID;
220   fUseAOD049=source.fUseAOD049;
221   if(fPidHF) delete fPidHF;
222   fPidHF=new AliAODPidHF(*(source.GetPidHF()));
223   fWhyRejection=source.fWhyRejection;
224   fEvRejectionBits=source.fEvRejectionBits;
225   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
226   fRecomputePrimVertex=source.fRecomputePrimVertex;
227   fUseMCVertex=source.fUseMCVertex;
228   fUsePhysicsSelection=source.fUsePhysicsSelection;
229   fOptPileup=source.fOptPileup;
230   fMinContrPileup=source.fMinContrPileup;
231   fMinDzPileup=source.fMinDzPileup;
232   fUseCentrality=source.fUseCentrality;
233   fMinCentrality=source.fMinCentrality;
234   fMaxCentrality=source.fMaxCentrality;
235   fFixRefs=source.fFixRefs;
236   fIsSelectedCuts=source.fIsSelectedCuts;
237   fIsSelectedPID=source.fIsSelectedPID;
238   fMinPtCand=source.fMinPtCand;
239   fMaxPtCand=source.fMaxPtCand;
240   fMaxRapidityCand=source.fMaxRapidityCand;
241   fKeepSignalMC=source.fKeepSignalMC;
242   fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
243   fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
244   fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
245   fMaxDiffTRKV0Centr=source.fMaxDiffTRKV0Centr;
246   fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
247   fCutOnzVertexSPD=source.fCutOnzVertexSPD;
248   fKinkReject=source.fKinkReject;
249   fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
250   if(fHistCentrDistr) delete fHistCentrDistr;
251   fUseCentrFlatteningInMC=source.fUseCentrFlatteningInMC;
252   if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
253
254   if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
255   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
256   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
257   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
258   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
259   
260   if(fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
261   if(f1CutMinNCrossedRowsTPCPtDep) delete f1CutMinNCrossedRowsTPCPtDep;
262   if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
263   fCutRatioClsOverCrossRowsTPC=source.fCutRatioClsOverCrossRowsTPC;
264   fCutRatioSignalNOverCrossRowsTPC=source.fCutRatioSignalNOverCrossRowsTPC;
265   PrintAll();
266
267   return *this;
268 }
269 //--------------------------------------------------------------------------
270 AliRDHFCuts::~AliRDHFCuts() {
271   //  
272   // Default Destructor
273   //
274   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
275   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
276   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
277   if(fCutsRD) {
278     delete [] fCutsRD;
279     fCutsRD=0;
280   }
281   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
282   if(fPidHF){ 
283     delete fPidHF;
284     fPidHF=0;
285   }
286   if(fHistCentrDistr)delete fHistCentrDistr;
287
288   if(f1CutMinNCrossedRowsTPCPtDep) {
289     delete f1CutMinNCrossedRowsTPCPtDep;
290     f1CutMinNCrossedRowsTPCPtDep = 0;
291   }
292
293 }
294 //---------------------------------------------------------------------------
295 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
296   //
297   // Centrality selection
298   //
299   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){    
300     AliWarning("Centrality estimator not valid");    
301     return 3;  
302   }else{    
303     Float_t centvalue=GetCentrality((AliAODEvent*)event);          
304     if (centvalue<-998.){//-999 if no centralityP
305       return 0;    
306     }else{      
307       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
308         return 2;      
309       }
310       // selection to flatten the centrality distribution
311       if(fHistCentrDistr){
312         if(!IsEventSelectedForCentrFlattening(centvalue))return 4;     
313       }
314     } 
315   }  
316   return 0;
317 }
318
319
320 //-------------------------------------------------
321 void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
322   // set the histo for centrality flattening 
323   // the centrality is flatten in the range minCentr,maxCentr
324   // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference 
325   //                positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
326   //                negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat) 
327   // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom 
328   
329   if(maxCentr<minCentr){
330     AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
331   }
332   
333   if(fHistCentrDistr)delete fHistCentrDistr;
334   fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
335   fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
336   Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
337   Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
338   fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
339   Double_t ref=0.,bincont=0.,binrefwidth=1.;
340   Int_t binref=0;
341   if(TMath::Abs(centrRef)<0.0001){
342     binref=fHistCentrDistr->GetMinimumBin();
343     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
344     ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;   
345   }
346   else if(centrRef>0.){
347     binref=h->FindBin(centrRef);
348     if(binref<1||binref>h->GetNbinsX()){
349       AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
350     }
351     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
352     ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
353   }
354   else{
355     if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
356     binref=fHistCentrDistr->GetMaximumBin();
357     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
358     ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;   
359   }
360   
361   for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
362     if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
363       bincont=h->GetBinContent(j);
364       fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
365       fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
366     }
367     else{
368       h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
369     }
370   }
371
372   fHistCentrDistr->SetBinContent(0,switchTRand);
373   return;
374
375 }
376
377 //-------------------------------------------------
378 Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
379   //
380   //  Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
381   //  Can be faster if it was required that fHistCentrDistr covers
382   //  exactly the desired centrality range (e.g. part of the lines below should be done during the 
383   // setting of the histo) and TH1::SetMinimum called 
384   //
385
386   if(!fHistCentrDistr) return kTRUE;
387   // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
388   //   if(maxbin>fHistCentrDistr->GetNbinsX()){
389   //     AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
390   //   }
391   
392   Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
393   Double_t bincont=fHistCentrDistr->GetBinContent(bin);
394   Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
395   
396   if(fHistCentrDistr->GetBinContent(0)<-0.9999){
397     if(gRandom->Uniform(1.)<bincont)return kTRUE;
398     return kFALSE;
399   }
400
401   if(centDigits*100.<bincont)return kTRUE;
402   return kFALSE;   
403
404 }
405 //---------------------------------------------------------------------------
406 void AliRDHFCuts::SetupPID(AliVEvent *event) {
407   // Set the PID response object in the AliAODPidHF
408   // in case of old PID sets the TPC dE/dx BB parameterization
409
410   if(fPidHF){
411     if(fPidHF->GetPidResponse()==0x0){
412       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
413       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
414       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
415       fPidHF->SetPidResponse(pidResp);
416     }
417     if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
418     if(fPidHF->GetOldPid()) {
419
420       Bool_t isMC=kFALSE;
421       TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
422       if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
423
424       // pp, from LHC10d onwards
425       if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
426          event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
427       // pp, 2011 low energy run
428       if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
429         fPidHF->SetppLowEn2011(kTRUE);
430         fPidHF->SetOnePad(kFALSE);
431       }
432       // PbPb LHC10h
433       if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
434       // MC
435       if(isMC) fPidHF->SetMC(kTRUE);
436       if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
437         fPidHF->SetMClowenpp2011(kTRUE);
438       fPidHF->SetBetheBloch();
439     }else{
440       // check that AliPIDResponse object was properly set in case of using OADB
441       if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
442     }
443   }
444 }
445 //---------------------------------------------------------------------------
446 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
447   //
448   // Event selection
449   // 
450   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
451
452
453
454   fWhyRejection=0;
455   fEvRejectionBits=0;
456   Bool_t accept=kTRUE;
457
458   if(fRecomputePrimVertex){
459     Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
460     if(!vertOK){
461       fWhyRejection=6;
462       return kFALSE;
463     }
464   }
465
466   // check if it's MC
467   Bool_t isMC=kFALSE;
468   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
469   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
470
471
472   SetupPID(event);
473
474   // trigger class
475   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
476   // don't do for MC and for PbPb 2010 data
477   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
478     if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) && 
479        (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
480       fWhyRejection=5;
481       fEvRejectionBits+=1<<kNotSelTrigger;
482       accept=kFALSE;
483     }
484   }
485
486   // TEMPORARY FIX FOR GetEvent
487   Int_t nTracks=((AliAODEvent*)event)->GetNTracks();
488   for(Int_t itr=0; itr<nTracks; itr++){
489     AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
490     tr->SetAODEvent((AliAODEvent*)event);
491   }
492
493   // TEMPORARY FIX FOR REFERENCES
494   // Fix references to daughter tracks
495   //  if(fFixRefs) {
496   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
497   //    fixer->FixReferences((AliAODEvent*)event);
498   //    delete fixer;
499   //  }
500   //
501
502
503   // physics selection requirements
504   if(fUsePhysicsSelection){
505     Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
506     if(!isSelected) {
507       if(accept) fWhyRejection=7;
508       fEvRejectionBits+=1<<kPhysicsSelection;
509       accept=kFALSE;
510     }else{
511       if(fUseOnlyOneTrigger){
512         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
513           if(accept) fWhyRejection=7;
514           fEvRejectionBits+=1<<kPhysicsSelection;
515           accept=kFALSE;
516         }
517       }
518     }
519   }
520
521   // vertex requirements
522    
523   const AliVVertex *vertex = event->GetPrimaryVertex();
524
525   if(!vertex){
526     accept=kFALSE;
527     fEvRejectionBits+=1<<kNoVertex;
528   }else{
529     TString title=vertex->GetTitle();
530     if(title.Contains("Z") && fMinVtxType>1){
531       accept=kFALSE;
532       fEvRejectionBits+=1<<kNoVertex;
533     }
534     else if(title.Contains("3D") && fMinVtxType>2){
535       accept=kFALSE;
536       fEvRejectionBits+=1<<kNoVertex;
537     }
538     if(vertex->GetNContributors()<fMinVtxContr){
539       accept=kFALSE;
540       fEvRejectionBits+=1<<kTooFewVtxContrib;
541     }
542     if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
543       fEvRejectionBits+=1<<kZVtxOutFid;
544       if(accept) fWhyRejection=6;
545       accept=kFALSE;
546     } 
547   }
548
549   if(fCutOnzVertexSPD>0){
550     const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
551     if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
552       accept=kFALSE;
553       fEvRejectionBits+=1<<kBadSPDVertex;
554     }else{
555       if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
556         fEvRejectionBits+=1<<kZVtxSPDOutFid;
557         if(accept) fWhyRejection=6;
558         accept=kFALSE;
559       } 
560       if(fCutOnzVertexSPD==2 && vertex){
561         if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
562           fEvRejectionBits+=1<<kZVtxSPDOutFid;
563           if(accept) fWhyRejection=6;
564           accept=kFALSE;
565         } 
566       }
567     }
568   }
569
570   // pile-up rejection
571   if(fOptPileup==kRejectPileupEvent){
572     Int_t cutc=(Int_t)fMinContrPileup;
573     Double_t cutz=(Double_t)fMinDzPileup;
574     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
575       if(accept) fWhyRejection=1;
576       fEvRejectionBits+=1<<kPileup;
577       accept=kFALSE;
578     }
579   }
580   else if(fOptPileup==kRejectMVPileupEvent){
581     AliAnalysisUtils utils;
582     Bool_t isPUMV = utils.IsPileUpMV(event);
583     if(isPUMV) {
584       if(accept) fWhyRejection=1;
585       fEvRejectionBits+=1<<kPileup;
586       accept=kFALSE;
587     }
588   }
589
590   // centrality selection
591   if (fUseCentrality!=kCentOff) {  
592     Int_t rejection=IsEventSelectedInCentrality(event);    
593     Bool_t okCent=kFALSE;
594     if(rejection==0) okCent=kTRUE;
595     if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
596     if(!okCent){      
597       if(accept) fWhyRejection=rejection;      
598       if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
599       else fEvRejectionBits+=1<<kOutsideCentrality;
600       accept=kFALSE;
601     }
602    
603   }
604
605   // PbPb2011 outliers in tracklets vs. VZERO and centTRK vs. centV0
606   if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
607     if(fRemoveTrackletOutliers){
608       Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
609       Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
610       Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
611       if(ntracklets<1000. && v0cent<cutval){
612         if(accept) fWhyRejection=2;      
613         fEvRejectionBits+=1<<kOutsideCentrality;
614          accept=kFALSE;
615       }
616     }
617     if(fMaxDiffTRKV0Centr>0.){
618       Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
619       Double_t trkcent=GetCentrality((AliAODEvent*)event,kCentTRK);
620       if(TMath::Abs(trkcent-v0cent)>fMaxDiffTRKV0Centr){
621         if(accept) fWhyRejection=1;
622         fEvRejectionBits+=1<<kBadTrackV0Correl;
623         accept=kFALSE;  
624       }
625     }
626   }
627
628   return accept;
629 }
630 //---------------------------------------------------------------------------
631 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const{
632   //
633   // Daughter tracks selection
634   // 
635   if(!fTrackCuts) return kTRUE;
636  
637   Int_t ndaughters = d->GetNDaughters();
638   AliAODVertex *vAOD = d->GetPrimaryVtx();
639   Double_t pos[3],cov[6];
640   vAOD->GetXYZ(pos);
641   vAOD->GetCovarianceMatrix(cov);
642   const AliESDVertex vESD(pos,cov,100.,100);
643   
644   Bool_t retval=kTRUE;
645   
646   for(Int_t idg=0; idg<ndaughters; idg++) {
647     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
648     if(!dgTrack) {retval = kFALSE; continue;}
649     //printf("charge %d\n",dgTrack->Charge());
650     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
651
652     if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
653       { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
654
655     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
656   }
657   
658   return retval;
659 }
660 //---------------------------------------------------------------------------
661 Bool_t AliRDHFCuts::CheckPtDepCrossedRows(TString rows,Bool_t print) const {
662   //
663   // Check the correctness of the string syntax
664   //
665   Bool_t retval=kTRUE;
666
667   if(!rows.Contains("pt")) {
668     if(print) AliError("string must contain \"pt\"");
669     retval= kFALSE;
670   }
671   return retval;
672 }
673 //---------------------------------------------------------------------------
674 void AliRDHFCuts::SetMinCrossedRowsTPCPtDep(const char *rows){
675   //
676   //Create the TFormula from TString for TPC crossed rows pT dependent cut 
677   //
678
679
680   // setting data member that describes the TPC crossed rows pT dependent cut 
681   fCutMinCrossedRowsTPCPtDep = rows;
682
683   // creating TFormula from TString
684    if(f1CutMinNCrossedRowsTPCPtDep){
685      delete f1CutMinNCrossedRowsTPCPtDep;
686      // resetting TFormula
687      f1CutMinNCrossedRowsTPCPtDep = 0;
688    }
689    if(!CheckPtDepCrossedRows(rows,kTRUE))return;   
690    
691    TString tmp(rows);
692    tmp.ReplaceAll("pt","x");
693    f1CutMinNCrossedRowsTPCPtDep = new TFormula("f1CutMinNCrossedRowsTPCPtDep",tmp.Data());
694
695    
696 }
697 //---------------------------------------------------------------------------
698 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const{
699   //
700   // Convert to ESDtrack, relate to vertex and check cuts
701   //
702   if(!cuts) return kTRUE;
703
704   if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
705
706
707   // convert to ESD track here
708   AliESDtrack esdTrack(track);
709   // set the TPC cluster info
710   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
711   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
712   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
713   // needed to calculate the impact parameters
714   esdTrack.RelateToVertex(primary,0.,3.);
715
716   //applying ESDtrackCut
717   if(!cuts->IsSelected(&esdTrack)) return kFALSE; 
718
719   //appliyng kink rejection
720   if(fKinkReject){
721    AliAODVertex *maybeKink=track->GetProdVertex();
722    if(maybeKink->GetType()==AliAODVertex::kKink) return kFALSE;
723   }
724
725   //appliyng TPC crossed rows pT dependent cut
726   if(f1CutMinNCrossedRowsTPCPtDep){
727     Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
728     if(nCrossedRowsTPC<f1CutMinNCrossedRowsTPCPtDep->Eval(esdTrack.Pt())) return kFALSE;
729   }
730   
731   //appliyng NTPCcls/NTPCcrossedRows cut
732   if(fCutRatioClsOverCrossRowsTPC){
733     Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
734     Float_t nClustersTPC = esdTrack.GetTPCNcls();
735     if(nCrossedRowsTPC!=0){ 
736       Float_t ratio = nClustersTPC/nCrossedRowsTPC;
737       if(ratio<fCutRatioClsOverCrossRowsTPC) return kFALSE;
738     }
739     else return kFALSE;
740   }
741
742   //appliyng TPCsignalN/NTPCcrossedRows cut
743   if(fCutRatioSignalNOverCrossRowsTPC){
744     Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
745     Float_t nTPCsignal = esdTrack.GetTPCsignalN();
746     if(nCrossedRowsTPC!=0){
747       Float_t ratio = nTPCsignal/nCrossedRowsTPC;
748       if(ratio<fCutRatioSignalNOverCrossRowsTPC) return kFALSE;
749     }
750     else return kFALSE;
751   }
752
753  
754   if(fOptPileup==kRejectTracksFromPileupVertex){
755     // to be implemented
756     // we need either to have here the AOD Event, 
757     // or to have the pileup vertex object
758   }
759
760   if(fApplySPDDeadPbPb2011){
761     Bool_t deadSPDLay1PbPb2011[20][4]={
762       {kTRUE,kTRUE,kTRUE,kTRUE},
763       {kTRUE,kTRUE,kTRUE,kTRUE},
764       {kTRUE,kTRUE,kTRUE,kTRUE},
765       {kTRUE,kTRUE,kTRUE,kTRUE},
766       {kTRUE,kTRUE,kTRUE,kTRUE},
767       {kFALSE,kFALSE,kTRUE,kTRUE},
768       {kTRUE,kTRUE,kFALSE,kFALSE},
769       {kTRUE,kTRUE,kTRUE,kTRUE},
770       {kFALSE,kFALSE,kTRUE,kTRUE},
771       {kTRUE,kTRUE,kTRUE,kTRUE},
772       {kTRUE,kTRUE,kFALSE,kFALSE},
773       {kTRUE,kTRUE,kTRUE,kTRUE},
774       {kFALSE,kFALSE,kFALSE,kFALSE},
775       {kFALSE,kFALSE,kTRUE,kTRUE},
776       {kFALSE,kFALSE,kFALSE,kFALSE},
777       {kFALSE,kFALSE,kFALSE,kFALSE},
778       {kTRUE,kTRUE,kTRUE,kTRUE},
779       {kTRUE,kTRUE,kFALSE,kFALSE},
780       {kFALSE,kFALSE,kFALSE,kFALSE},
781       {kFALSE,kFALSE,kFALSE,kFALSE}
782     };
783     Bool_t deadSPDLay2PbPb2011[40][4]={
784       {kTRUE,kTRUE,kTRUE,kTRUE},
785       {kTRUE,kTRUE,kTRUE,kTRUE},
786       {kTRUE,kTRUE,kTRUE,kTRUE},
787       {kTRUE,kTRUE,kTRUE,kTRUE},
788       {kTRUE,kTRUE,kTRUE,kTRUE},
789       {kTRUE,kTRUE,kTRUE,kTRUE},
790       {kTRUE,kTRUE,kTRUE,kTRUE},
791       {kTRUE,kTRUE,kTRUE,kTRUE},
792       {kTRUE,kTRUE,kTRUE,kTRUE},
793       {kTRUE,kTRUE,kTRUE,kTRUE},
794       {kTRUE,kTRUE,kTRUE,kTRUE},
795       {kTRUE,kTRUE,kTRUE,kTRUE},
796       {kFALSE,kFALSE,kFALSE,kFALSE},
797       {kFALSE,kFALSE,kTRUE,kTRUE},
798       {kTRUE,kTRUE,kTRUE,kTRUE},
799       {kTRUE,kTRUE,kTRUE,kTRUE},
800       {kTRUE,kTRUE,kFALSE,kFALSE},
801       {kTRUE,kTRUE,kTRUE,kTRUE},
802       {kTRUE,kTRUE,kTRUE,kTRUE},
803       {kTRUE,kTRUE,kTRUE,kTRUE},
804       {kFALSE,kFALSE,kFALSE,kFALSE},
805       {kFALSE,kFALSE,kFALSE,kFALSE},
806       {kTRUE,kTRUE,kTRUE,kTRUE},
807       {kTRUE,kTRUE,kTRUE,kTRUE},
808       {kFALSE,kFALSE,kFALSE,kFALSE},
809       {kFALSE,kFALSE,kFALSE,kFALSE},
810       {kTRUE,kTRUE,kTRUE,kTRUE},
811       {kTRUE,kTRUE,kTRUE,kTRUE},
812       {kFALSE,kFALSE,kFALSE,kFALSE},
813       {kFALSE,kFALSE,kFALSE,kFALSE},
814       {kFALSE,kFALSE,kFALSE,kFALSE},
815       {kFALSE,kFALSE,kFALSE,kFALSE},
816       {kTRUE,kTRUE,kTRUE,kTRUE},
817       {kTRUE,kTRUE,kTRUE,kTRUE},
818       {kTRUE,kTRUE,kFALSE,kFALSE},
819       {kTRUE,kTRUE,kTRUE,kTRUE},
820       {kFALSE,kFALSE,kFALSE,kFALSE},
821       {kFALSE,kFALSE,kFALSE,kFALSE},
822       {kFALSE,kFALSE,kFALSE,kFALSE},
823       {kFALSE,kFALSE,kFALSE,kFALSE}     
824     };
825     Double_t xyz1[3],xyz2[3];
826     esdTrack.GetXYZAt(3.9,0.,xyz1);
827     esdTrack.GetXYZAt(7.6,0.,xyz2);
828     Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
829     if(phi1<0) phi1+=2*TMath::Pi();
830     Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
831     Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
832     if(phi2<0) phi2+=2*TMath::Pi();
833     Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
834     Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
835     Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
836     Bool_t lay1ok=kFALSE;
837     if(mod1>=0 && mod1<4 && lad1<20){
838       lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
839     }
840     Bool_t lay2ok=kFALSE;
841     if(mod2>=0 && mod2<4 && lad2<40){
842       lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
843     }
844     if(!lay1ok && !lay2ok) return kFALSE;
845   }
846
847   return kTRUE; 
848 }
849 //---------------------------------------------------------------------------
850 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
851   // Set the pt bins
852
853   if(fPtBinLimits) {
854     delete [] fPtBinLimits;
855     fPtBinLimits = NULL;
856     printf("Changing the pt bins\n");
857   }
858
859   if(nPtBinLimits != fnPtBins+1){
860     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
861     SetNPtBins(nPtBinLimits-1);
862   }
863
864   fnPtBinLimits = nPtBinLimits;
865   SetGlobalIndex();
866   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
867   fPtBinLimits = new Float_t[fnPtBinLimits];
868   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
869
870   return;
871 }
872 //---------------------------------------------------------------------------
873 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
874   // Set the variable names
875
876   if(fVarNames) {
877     delete [] fVarNames;
878     fVarNames = NULL;
879     //printf("Changing the variable names\n");
880   }
881   if(nVars!=fnVars){
882     printf("Wrong number of variables: it has to be %d\n",fnVars);
883     return;
884   }
885   //fnVars=nVars;
886   fVarNames = new TString[nVars];
887   fIsUpperCut = new Bool_t[nVars];
888   for(Int_t iv=0; iv<nVars; iv++) {
889     fVarNames[iv] = varNames[iv];
890     fIsUpperCut[iv] = isUpperCut[iv];
891   }
892
893   return;
894 }
895 //---------------------------------------------------------------------------
896 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
897   // Set the variables to be used for cuts optimization
898
899   if(fVarsForOpt) {
900     delete [] fVarsForOpt;
901     fVarsForOpt = NULL;
902     //printf("Changing the variables for cut optimization\n");
903   }
904   
905   if(nVars==0){//!=fnVars) {
906     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
907     return;
908   } 
909   
910   fnVarsForOpt = 0;
911   fVarsForOpt = new Bool_t[fnVars];
912   for(Int_t iv=0; iv<fnVars; iv++) {
913     fVarsForOpt[iv]=forOpt[iv];
914     if(fVarsForOpt[iv]) fnVarsForOpt++;
915   }
916
917   return;
918 }
919
920 //---------------------------------------------------------------------------
921 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
922   //
923   // set centrality estimator  
924   //
925   fUseCentrality=flag;
926   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
927  
928   return;
929 }
930
931
932 //---------------------------------------------------------------------------
933 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
934   //
935   // store the cuts
936   //
937   if(nVars!=fnVars) {
938     printf("Wrong number of variables: it has to be %d\n",fnVars);
939     AliFatal("exiting");
940   } 
941   if(nPtBins!=fnPtBins) {
942     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
943     AliFatal("exiting");
944   } 
945
946   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
947   
948
949   for(Int_t iv=0; iv<fnVars; iv++) {
950
951     for(Int_t ib=0; ib<fnPtBins; ib++) {
952
953       //check
954       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
955         cout<<"Overflow, exit..."<<endl;
956         return;
957       }
958
959       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
960
961     }
962   }
963   return;
964 }
965 //---------------------------------------------------------------------------
966 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
967   //
968   // store the cuts
969   //
970   if(glIndex != fGlobalIndex){
971     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
972     AliFatal("exiting");
973   }
974   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
975
976   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
977     fCutsRD[iGl] = cutsRDGlob[iGl];
978   }
979   return;
980 }
981 //---------------------------------------------------------------------------
982 void AliRDHFCuts::PrintAll() const {
983   //
984   // print all cuts values
985   // 
986
987   printf("Minimum vtx type %d\n",fMinVtxType);
988   printf("Minimum vtx contr %d\n",fMinVtxContr);
989   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
990   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
991   printf("Use PID %d  OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
992   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
993   printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
994   printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
995   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
996   if(fOptPileup==1) printf(" -- Reject pileup event");
997   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
998   if(fUseCentrality>0) {
999     TString estimator="";
1000     if(fUseCentrality==1) estimator = "V0";
1001     if(fUseCentrality==2) estimator = "Tracks";
1002     if(fUseCentrality==3) estimator = "Tracklets";
1003     if(fUseCentrality==4) estimator = "SPD clusters outer"; 
1004     if(fUseCentrality==5) estimator = "ZNA"; 
1005     if(fUseCentrality==6) estimator = "ZPA"; 
1006     if(fUseCentrality==7) estimator = "V0A"; 
1007     printf("Centrality class considered: %.1f-%.1f, estimated with %s\n",fMinCentrality,fMaxCentrality,estimator.Data());
1008   }
1009   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
1010
1011   if(fCutRatioClsOverCrossRowsTPC) printf("N TPC Clusters > %f N TPC Crossed Rows\n", fCutRatioClsOverCrossRowsTPC);
1012   if(fCutRatioSignalNOverCrossRowsTPC) printf("N TPC Points for dE/dx > %f N TPC Crossed Rows\n", fCutRatioSignalNOverCrossRowsTPC);
1013   if(f1CutMinNCrossedRowsTPCPtDep) printf("N TPC Crossed Rows pT-dependent cut: %s\n", fCutMinCrossedRowsTPCPtDep.Data());
1014
1015   if(fVarNames){
1016     cout<<"Array of variables"<<endl;
1017     for(Int_t iv=0;iv<fnVars;iv++){
1018       cout<<fVarNames[iv]<<"\t";
1019     }
1020     cout<<endl;
1021   }
1022   if(fVarsForOpt){
1023     cout<<"Array of optimization"<<endl;
1024     for(Int_t iv=0;iv<fnVars;iv++){
1025       cout<<fVarsForOpt[iv]<<"\t";
1026     }
1027     cout<<endl;
1028   }
1029   if(fIsUpperCut){
1030     cout<<"Array of upper/lower cut"<<endl;
1031    for(Int_t iv=0;iv<fnVars;iv++){
1032      cout<<fIsUpperCut[iv]<<"\t";
1033    }
1034    cout<<endl;
1035   }
1036   if(fPtBinLimits){
1037     cout<<"Array of ptbin limits"<<endl;
1038     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
1039       cout<<fPtBinLimits[ib]<<"\t";
1040     }
1041     cout<<endl;
1042   }
1043   if(fCutsRD){
1044     cout<<"Matrix of cuts"<<endl;
1045    for(Int_t iv=0;iv<fnVars;iv++){
1046      for(Int_t ib=0;ib<fnPtBins;ib++){
1047        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
1048      } 
1049      cout<<endl;
1050    }
1051    cout<<endl;
1052   }
1053   fPidHF->PrintAll();
1054   return;
1055 }
1056
1057 //--------------------------------------------------------------------------
1058 void AliRDHFCuts::PrintTrigger() const{
1059   // print the trigger selection 
1060
1061   printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
1062
1063   cout<<" Trigger selection pattern: ";
1064   if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
1065   if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
1066   if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
1067   if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
1068   if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
1069   if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
1070   if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
1071   if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
1072   if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
1073   if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
1074   cout << endl<< endl;
1075
1076 }
1077
1078 //---------------------------------------------------------------------------
1079 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
1080   //
1081   // get the cuts
1082   //
1083
1084   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
1085
1086
1087   Int_t iv,ib;
1088   if(!cutsRD) {
1089     //cout<<"Initialization..."<<endl;
1090     cutsRD=new Float_t*[fnVars];
1091     for(iv=0; iv<fnVars; iv++) {
1092       cutsRD[iv] = new Float_t[fnPtBins];
1093     }
1094   }
1095   
1096   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
1097     GetVarPtIndex(iGlobal,iv,ib);
1098     cutsRD[iv][ib] = fCutsRD[iGlobal];
1099   }
1100
1101   return;
1102 }
1103
1104 //---------------------------------------------------------------------------
1105 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
1106   //
1107   // give the global index from variable and pt bin
1108   //
1109   return iPtBin*fnVars+iVar;
1110 }
1111
1112 //---------------------------------------------------------------------------
1113 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
1114   //
1115   //give the index of the variable and of the pt bin from the global index
1116   //
1117   iPtBin=(Int_t)iGlob/fnVars;
1118   iVar=iGlob%fnVars;
1119
1120   return;
1121 }
1122
1123 //---------------------------------------------------------------------------
1124 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1125   //
1126   //give the pt bin where the pt lies.
1127   //
1128   Int_t ptbin=-1;
1129   if(pt<fPtBinLimits[0])return ptbin;
1130   for (Int_t i=0;i<fnPtBins;i++){
1131     if(pt<fPtBinLimits[i+1]) {
1132       ptbin=i;
1133       break;
1134     }
1135   }
1136   return ptbin;
1137 }
1138 //-------------------------------------------------------------------
1139 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1140   // 
1141   // Give the value of cut set for the variable iVar and the pt bin iPtBin
1142   //
1143   if(!fCutsRD){
1144     cout<<"Cuts not iniziaisez yet"<<endl;
1145     return 0;
1146   }
1147   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1148 }
1149 //-------------------------------------------------------------------
1150 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1151   //
1152   // Get centrality percentile
1153   //
1154
1155   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1156   if(mcArray) {fUseAOD049=kFALSE;}
1157
1158   AliAODHeader *header=aodEvent->GetHeader();
1159   AliCentrality *centrality=header->GetCentralityP();
1160   Float_t cent=-999.;
1161   Bool_t isSelRun=kFALSE;
1162   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1163   if(!centrality) return cent;
1164   else{
1165     if (estimator==kCentV0M){
1166       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1167       if(cent<0){
1168         Int_t quality = centrality->GetQuality();
1169         if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1170           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1171         }else{
1172           Int_t runnum=aodEvent->GetRunNumber();
1173           for(Int_t ir=0;ir<5;ir++){
1174             if(runnum==selRun[ir]){
1175               isSelRun=kTRUE;
1176               break;
1177             }
1178           }
1179           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1180         }
1181       }
1182
1183       //temporary fix for AOD049 outliers
1184       if(fUseAOD049&&cent>=0){
1185         Float_t v0=0;
1186         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1187         v0+=aodV0->GetMTotV0A();
1188         v0+=aodV0->GetMTotV0C();
1189         if(cent==0&&v0<19500)return -1;//filtering issue
1190         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1191         Float_t val= 1.30552 +  0.147931 * v0;
1192         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};
1193         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1194       }
1195     }
1196     else {
1197        if (estimator==kCentTRK) {
1198         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1199         if(cent<0){
1200           Int_t quality = centrality->GetQuality();
1201           if(quality<=1){
1202             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1203           }else{
1204             Int_t runnum=aodEvent->GetRunNumber();
1205             for(Int_t ir=0;ir<5;ir++){
1206               if(runnum==selRun[ir]){
1207                 isSelRun=kTRUE;
1208                 break;
1209               }
1210             }
1211             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1212           }
1213         }
1214        }
1215       else{
1216         if (estimator==kCentTKL){
1217           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1218           if(cent<0){
1219             Int_t quality = centrality->GetQuality();
1220             if(quality<=1){
1221               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1222             }else{
1223               Int_t runnum=aodEvent->GetRunNumber();
1224               for(Int_t ir=0;ir<5;ir++){
1225                 if(runnum==selRun[ir]){
1226                   isSelRun=kTRUE;
1227                   break;
1228             }
1229               }
1230               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1231             }   
1232           }
1233         }
1234         else{
1235           if (estimator==kCentCL1){
1236             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1237             if(cent<0){
1238               Int_t quality = centrality->GetQuality();
1239               if(quality<=1){
1240                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1241               }else{
1242                 Int_t runnum=aodEvent->GetRunNumber();
1243                 for(Int_t ir=0;ir<5;ir++){
1244                   if(runnum==selRun[ir]){
1245                     isSelRun=kTRUE;
1246                     break;
1247                   }
1248                 }
1249                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1250               }
1251             }
1252           }
1253         else{
1254           if (estimator==kCentZNA){
1255             cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1256             if(cent<0){
1257               Int_t quality = centrality->GetQuality();
1258               if(quality<=1){
1259                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1260               }else{
1261                 Int_t runnum=aodEvent->GetRunNumber();
1262                 for(Int_t ir=0;ir<5;ir++){
1263                   if(runnum==selRun[ir]){
1264                     isSelRun=kTRUE;
1265                     break;
1266                   }
1267                 }
1268                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1269               }
1270             }
1271           }
1272         else{
1273           if (estimator==kCentZPA){
1274             cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1275             if(cent<0){
1276               Int_t quality = centrality->GetQuality();
1277               if(quality<=1){
1278                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1279               }else{
1280                 Int_t runnum=aodEvent->GetRunNumber();
1281                 for(Int_t ir=0;ir<5;ir++){
1282                   if(runnum==selRun[ir]){
1283                     isSelRun=kTRUE;
1284                     break;
1285                   }
1286                 }
1287                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1288               }
1289             }
1290           }
1291         else{
1292           if (estimator==kCentV0A){
1293             cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1294             if(cent<0){
1295               Int_t quality = centrality->GetQuality();
1296               if(quality<=1){
1297                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1298               }else{
1299                 Int_t runnum=aodEvent->GetRunNumber();
1300                 for(Int_t ir=0;ir<5;ir++){
1301                   if(runnum==selRun[ir]){
1302                     isSelRun=kTRUE;
1303                     break;
1304                   }
1305                 }
1306                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1307               }
1308             }
1309           }
1310           else {
1311             AliWarning("Centrality estimator not valid");
1312             
1313           }
1314         }
1315     }
1316     }
1317     }
1318     }
1319     }
1320   }
1321   return cent;
1322 }
1323 //-------------------------------------------------------------------
1324 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1325   //
1326   // Compare two cuts objects
1327   //
1328
1329   Bool_t areEqual=kTRUE;
1330
1331   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1332
1333   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1334
1335   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1336
1337   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1338
1339   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1340
1341   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1342   if(fTrackCuts){
1343     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1344
1345     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1346
1347     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1348
1349     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;}
1350   }
1351
1352   if(fCutsRD) {
1353    for(Int_t iv=0;iv<fnVars;iv++) {
1354      for(Int_t ib=0;ib<fnPtBins;ib++) {
1355        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1356          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1357          areEqual=kFALSE;
1358        }
1359      }
1360    }
1361   }
1362
1363   return areEqual;
1364 }
1365 //---------------------------------------------------------------------------
1366 void AliRDHFCuts::MakeTable() const {
1367   //
1368   // print cuts values in table format
1369   // 
1370
1371         TString ptString = "pT range";
1372         if(fVarNames && fPtBinLimits && fCutsRD){
1373                 TString firstLine(Form("*       %-15s",ptString.Data()));
1374                 for (Int_t ivar=0; ivar<fnVars; ivar++){
1375                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
1376                         if (ivar == fnVars){
1377                                 firstLine+="*\n";
1378                         }
1379                 }
1380                 Printf("%s",firstLine.Data());
1381                 
1382                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1383                         TString line;
1384                         if (ipt==fnPtBins-1){
1385                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
1386                         }
1387                         else{
1388                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1389                         }
1390                         for (Int_t ivar=0; ivar<fnVars; ivar++){
1391                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1392                         }
1393                         Printf("%s",line.Data());
1394                 }
1395
1396         }
1397
1398
1399   return;
1400 }
1401 //--------------------------------------------------------------------------
1402 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1403                                         AliAODEvent *aod) const
1404 {
1405   //
1406   // Recalculate primary vertex without daughters
1407   //
1408
1409   if(!aod) {
1410     AliError("Can not remove daughters from vertex without AOD event");
1411     return 0;
1412   }   
1413
1414   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1415   if(!recvtx){
1416     AliDebug(2,"Removal of daughter tracks failed");
1417     return kFALSE;
1418   }
1419
1420
1421   //set recalculed primary vertex
1422   d->SetOwnPrimaryVtx(recvtx);
1423   delete recvtx;
1424
1425   return kTRUE;
1426 }
1427 //--------------------------------------------------------------------------
1428 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1429 {
1430   //
1431   // Recalculate primary vertex without daughters
1432   //
1433
1434   if(!aod) {
1435     AliError("Can not get MC vertex without AOD event");
1436     return kFALSE;
1437   }   
1438
1439   // load MC header
1440   AliAODMCHeader *mcHeader = 
1441     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1442   if(!mcHeader) {
1443     AliError("Can not get MC vertex without AODMCHeader event");
1444     return kFALSE;
1445   }
1446   Double_t pos[3];
1447   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1448   mcHeader->GetVertex(pos);
1449   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1450
1451   if(!recvtx){
1452     AliDebug(2,"Removal of daughter tracks failed");
1453     return kFALSE;
1454   }
1455
1456   //set recalculed primary vertex
1457   d->SetOwnPrimaryVtx(recvtx);
1458
1459   d->RecalculateImpPars(recvtx,aod);
1460
1461   delete recvtx;
1462
1463   return kTRUE;
1464 }
1465 //--------------------------------------------------------------------------
1466 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1467                                      AliAODEvent *aod,
1468                                      AliAODVertex *origownvtx) const
1469 {
1470   //
1471   // Clean-up own primary vertex if needed
1472   //
1473
1474   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1475     d->UnsetOwnPrimaryVtx();
1476     if(origownvtx) {
1477       d->SetOwnPrimaryVtx(origownvtx);
1478       delete origownvtx; origownvtx=NULL;
1479     }
1480     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1481   } else {
1482     if(origownvtx) {
1483       delete origownvtx; origownvtx=NULL;
1484     }
1485   }
1486   return;
1487 }
1488 //--------------------------------------------------------------------------
1489 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const 
1490 {
1491   //
1492   // Checks if this candidate is matched to MC signal
1493   // 
1494
1495   if(!aod) return kFALSE;
1496
1497   // get the MC array
1498   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1499
1500   if(!mcArray) return kFALSE;
1501
1502   // try to match  
1503   Int_t label = d->MatchToMC(pdg,mcArray);
1504   
1505   if(label>=0) {
1506     //printf("MATCH!\n");
1507     return kTRUE;
1508   }
1509
1510   return kFALSE;
1511 }
1512
1513
1514 //--------------------------------------------------------------------------
1515 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1516   // recompute event primary vertex from AOD tracks
1517
1518    AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1519    vertexer->SetITSMode();
1520    vertexer->SetMinClusters(3);
1521
1522    AliAODVertex* pvtx=event->GetPrimaryVertex(); 
1523    if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1524      Float_t diamondcovxy[3];
1525      event->GetDiamondCovXY(diamondcovxy);
1526      Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1527      Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1528      AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1529      vertexer->SetVtxStart(diamond);
1530      delete diamond; diamond=NULL;
1531    }
1532
1533    AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
1534    if(!vertexESD) return kFALSE;
1535    if(vertexESD->GetNContributors()<=0) { 
1536      //AliDebug(2,"vertexing failed"); 
1537      delete vertexESD; vertexESD=NULL;
1538      return kFALSE;
1539    }
1540    delete vertexer; vertexer=NULL;
1541
1542    // convert to AliAODVertex
1543    Double_t pos[3],cov[6],chi2perNDF;
1544    vertexESD->GetXYZ(pos); // position
1545    vertexESD->GetCovMatrix(cov); //covariance matrix
1546    chi2perNDF = vertexESD->GetChi2toNDF();
1547    delete vertexESD; vertexESD=NULL;
1548    
1549    pvtx->SetPosition(pos[0],pos[1],pos[2]);
1550    pvtx->SetChi2perNDF(chi2perNDF);
1551    pvtx->SetCovMatrix(cov);
1552
1553    return kTRUE;
1554 }