Usage of multi-vertexer to tag pileup events in D meson analysis (Zaida)
[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   return;
1054 }
1055
1056 //--------------------------------------------------------------------------
1057 void AliRDHFCuts::PrintTrigger() const{
1058   // print the trigger selection 
1059
1060   printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
1061
1062   cout<<" Trigger selection pattern: ";
1063   if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
1064   if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
1065   if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
1066   if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
1067   if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
1068   if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
1069   if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
1070   if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
1071   if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
1072   if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
1073   cout << endl<< endl;
1074
1075 }
1076
1077 //---------------------------------------------------------------------------
1078 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
1079   //
1080   // get the cuts
1081   //
1082
1083   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
1084
1085
1086   Int_t iv,ib;
1087   if(!cutsRD) {
1088     //cout<<"Initialization..."<<endl;
1089     cutsRD=new Float_t*[fnVars];
1090     for(iv=0; iv<fnVars; iv++) {
1091       cutsRD[iv] = new Float_t[fnPtBins];
1092     }
1093   }
1094   
1095   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
1096     GetVarPtIndex(iGlobal,iv,ib);
1097     cutsRD[iv][ib] = fCutsRD[iGlobal];
1098   }
1099
1100   return;
1101 }
1102
1103 //---------------------------------------------------------------------------
1104 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
1105   //
1106   // give the global index from variable and pt bin
1107   //
1108   return iPtBin*fnVars+iVar;
1109 }
1110
1111 //---------------------------------------------------------------------------
1112 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
1113   //
1114   //give the index of the variable and of the pt bin from the global index
1115   //
1116   iPtBin=(Int_t)iGlob/fnVars;
1117   iVar=iGlob%fnVars;
1118
1119   return;
1120 }
1121
1122 //---------------------------------------------------------------------------
1123 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1124   //
1125   //give the pt bin where the pt lies.
1126   //
1127   Int_t ptbin=-1;
1128   if(pt<fPtBinLimits[0])return ptbin;
1129   for (Int_t i=0;i<fnPtBins;i++){
1130     if(pt<fPtBinLimits[i+1]) {
1131       ptbin=i;
1132       break;
1133     }
1134   }
1135   return ptbin;
1136 }
1137 //-------------------------------------------------------------------
1138 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1139   // 
1140   // Give the value of cut set for the variable iVar and the pt bin iPtBin
1141   //
1142   if(!fCutsRD){
1143     cout<<"Cuts not iniziaisez yet"<<endl;
1144     return 0;
1145   }
1146   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1147 }
1148 //-------------------------------------------------------------------
1149 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1150   //
1151   // Get centrality percentile
1152   //
1153
1154   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1155   if(mcArray) {fUseAOD049=kFALSE;}
1156
1157   AliAODHeader *header=aodEvent->GetHeader();
1158   AliCentrality *centrality=header->GetCentralityP();
1159   Float_t cent=-999.;
1160   Bool_t isSelRun=kFALSE;
1161   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1162   if(!centrality) return cent;
1163   else{
1164     if (estimator==kCentV0M){
1165       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1166       if(cent<0){
1167         Int_t quality = centrality->GetQuality();
1168         if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1169           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1170         }else{
1171           Int_t runnum=aodEvent->GetRunNumber();
1172           for(Int_t ir=0;ir<5;ir++){
1173             if(runnum==selRun[ir]){
1174               isSelRun=kTRUE;
1175               break;
1176             }
1177           }
1178           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1179         }
1180       }
1181
1182       //temporary fix for AOD049 outliers
1183       if(fUseAOD049&&cent>=0){
1184         Float_t v0=0;
1185         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1186         v0+=aodV0->GetMTotV0A();
1187         v0+=aodV0->GetMTotV0C();
1188         if(cent==0&&v0<19500)return -1;//filtering issue
1189         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1190         Float_t val= 1.30552 +  0.147931 * v0;
1191         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};
1192         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1193       }
1194     }
1195     else {
1196        if (estimator==kCentTRK) {
1197         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1198         if(cent<0){
1199           Int_t quality = centrality->GetQuality();
1200           if(quality<=1){
1201             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1202           }else{
1203             Int_t runnum=aodEvent->GetRunNumber();
1204             for(Int_t ir=0;ir<5;ir++){
1205               if(runnum==selRun[ir]){
1206                 isSelRun=kTRUE;
1207                 break;
1208               }
1209             }
1210             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1211           }
1212         }
1213        }
1214       else{
1215         if (estimator==kCentTKL){
1216           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1217           if(cent<0){
1218             Int_t quality = centrality->GetQuality();
1219             if(quality<=1){
1220               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1221             }else{
1222               Int_t runnum=aodEvent->GetRunNumber();
1223               for(Int_t ir=0;ir<5;ir++){
1224                 if(runnum==selRun[ir]){
1225                   isSelRun=kTRUE;
1226                   break;
1227             }
1228               }
1229               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1230             }   
1231           }
1232         }
1233         else{
1234           if (estimator==kCentCL1){
1235             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1236             if(cent<0){
1237               Int_t quality = centrality->GetQuality();
1238               if(quality<=1){
1239                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1240               }else{
1241                 Int_t runnum=aodEvent->GetRunNumber();
1242                 for(Int_t ir=0;ir<5;ir++){
1243                   if(runnum==selRun[ir]){
1244                     isSelRun=kTRUE;
1245                     break;
1246                   }
1247                 }
1248                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1249               }
1250             }
1251           }
1252         else{
1253           if (estimator==kCentZNA){
1254             cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1255             if(cent<0){
1256               Int_t quality = centrality->GetQuality();
1257               if(quality<=1){
1258                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1259               }else{
1260                 Int_t runnum=aodEvent->GetRunNumber();
1261                 for(Int_t ir=0;ir<5;ir++){
1262                   if(runnum==selRun[ir]){
1263                     isSelRun=kTRUE;
1264                     break;
1265                   }
1266                 }
1267                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1268               }
1269             }
1270           }
1271         else{
1272           if (estimator==kCentZPA){
1273             cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1274             if(cent<0){
1275               Int_t quality = centrality->GetQuality();
1276               if(quality<=1){
1277                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1278               }else{
1279                 Int_t runnum=aodEvent->GetRunNumber();
1280                 for(Int_t ir=0;ir<5;ir++){
1281                   if(runnum==selRun[ir]){
1282                     isSelRun=kTRUE;
1283                     break;
1284                   }
1285                 }
1286                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1287               }
1288             }
1289           }
1290         else{
1291           if (estimator==kCentV0A){
1292             cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1293             if(cent<0){
1294               Int_t quality = centrality->GetQuality();
1295               if(quality<=1){
1296                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1297               }else{
1298                 Int_t runnum=aodEvent->GetRunNumber();
1299                 for(Int_t ir=0;ir<5;ir++){
1300                   if(runnum==selRun[ir]){
1301                     isSelRun=kTRUE;
1302                     break;
1303                   }
1304                 }
1305                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1306               }
1307             }
1308           }
1309           else {
1310             AliWarning("Centrality estimator not valid");
1311             
1312           }
1313         }
1314     }
1315     }
1316     }
1317     }
1318     }
1319   }
1320   return cent;
1321 }
1322 //-------------------------------------------------------------------
1323 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1324   //
1325   // Compare two cuts objects
1326   //
1327
1328   Bool_t areEqual=kTRUE;
1329
1330   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1331
1332   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1333
1334   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1335
1336   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1337
1338   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1339
1340   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1341   if(fTrackCuts){
1342     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1343
1344     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1345
1346     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1347
1348     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;}
1349   }
1350
1351   if(fCutsRD) {
1352    for(Int_t iv=0;iv<fnVars;iv++) {
1353      for(Int_t ib=0;ib<fnPtBins;ib++) {
1354        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1355          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1356          areEqual=kFALSE;
1357        }
1358      }
1359    }
1360   }
1361
1362   return areEqual;
1363 }
1364 //---------------------------------------------------------------------------
1365 void AliRDHFCuts::MakeTable() const {
1366   //
1367   // print cuts values in table format
1368   // 
1369
1370         TString ptString = "pT range";
1371         if(fVarNames && fPtBinLimits && fCutsRD){
1372                 TString firstLine(Form("*       %-15s",ptString.Data()));
1373                 for (Int_t ivar=0; ivar<fnVars; ivar++){
1374                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
1375                         if (ivar == fnVars){
1376                                 firstLine+="*\n";
1377                         }
1378                 }
1379                 Printf("%s",firstLine.Data());
1380                 
1381                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1382                         TString line;
1383                         if (ipt==fnPtBins-1){
1384                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
1385                         }
1386                         else{
1387                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1388                         }
1389                         for (Int_t ivar=0; ivar<fnVars; ivar++){
1390                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1391                         }
1392                         Printf("%s",line.Data());
1393                 }
1394
1395         }
1396
1397
1398   return;
1399 }
1400 //--------------------------------------------------------------------------
1401 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1402                                         AliAODEvent *aod) const
1403 {
1404   //
1405   // Recalculate primary vertex without daughters
1406   //
1407
1408   if(!aod) {
1409     AliError("Can not remove daughters from vertex without AOD event");
1410     return 0;
1411   }   
1412
1413   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1414   if(!recvtx){
1415     AliDebug(2,"Removal of daughter tracks failed");
1416     return kFALSE;
1417   }
1418
1419
1420   //set recalculed primary vertex
1421   d->SetOwnPrimaryVtx(recvtx);
1422   delete recvtx;
1423
1424   return kTRUE;
1425 }
1426 //--------------------------------------------------------------------------
1427 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1428 {
1429   //
1430   // Recalculate primary vertex without daughters
1431   //
1432
1433   if(!aod) {
1434     AliError("Can not get MC vertex without AOD event");
1435     return kFALSE;
1436   }   
1437
1438   // load MC header
1439   AliAODMCHeader *mcHeader = 
1440     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1441   if(!mcHeader) {
1442     AliError("Can not get MC vertex without AODMCHeader event");
1443     return kFALSE;
1444   }
1445   Double_t pos[3];
1446   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1447   mcHeader->GetVertex(pos);
1448   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1449
1450   if(!recvtx){
1451     AliDebug(2,"Removal of daughter tracks failed");
1452     return kFALSE;
1453   }
1454
1455   //set recalculed primary vertex
1456   d->SetOwnPrimaryVtx(recvtx);
1457
1458   d->RecalculateImpPars(recvtx,aod);
1459
1460   delete recvtx;
1461
1462   return kTRUE;
1463 }
1464 //--------------------------------------------------------------------------
1465 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1466                                      AliAODEvent *aod,
1467                                      AliAODVertex *origownvtx) const
1468 {
1469   //
1470   // Clean-up own primary vertex if needed
1471   //
1472
1473   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1474     d->UnsetOwnPrimaryVtx();
1475     if(origownvtx) {
1476       d->SetOwnPrimaryVtx(origownvtx);
1477       delete origownvtx; origownvtx=NULL;
1478     }
1479     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1480   } else {
1481     if(origownvtx) {
1482       delete origownvtx; origownvtx=NULL;
1483     }
1484   }
1485   return;
1486 }
1487 //--------------------------------------------------------------------------
1488 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const 
1489 {
1490   //
1491   // Checks if this candidate is matched to MC signal
1492   // 
1493
1494   if(!aod) return kFALSE;
1495
1496   // get the MC array
1497   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1498
1499   if(!mcArray) return kFALSE;
1500
1501   // try to match  
1502   Int_t label = d->MatchToMC(pdg,mcArray);
1503   
1504   if(label>=0) {
1505     //printf("MATCH!\n");
1506     return kTRUE;
1507   }
1508
1509   return kFALSE;
1510 }
1511
1512
1513 //--------------------------------------------------------------------------
1514 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1515   // recompute event primary vertex from AOD tracks
1516
1517    AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1518    vertexer->SetITSMode();
1519    vertexer->SetMinClusters(3);
1520
1521    AliAODVertex* pvtx=event->GetPrimaryVertex(); 
1522    if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1523      Float_t diamondcovxy[3];
1524      event->GetDiamondCovXY(diamondcovxy);
1525      Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1526      Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1527      AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1528      vertexer->SetVtxStart(diamond);
1529      delete diamond; diamond=NULL;
1530    }
1531
1532    AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
1533    if(!vertexESD) return kFALSE;
1534    if(vertexESD->GetNContributors()<=0) { 
1535      //AliDebug(2,"vertexing failed"); 
1536      delete vertexESD; vertexESD=NULL;
1537      return kFALSE;
1538    }
1539    delete vertexer; vertexer=NULL;
1540
1541    // convert to AliAODVertex
1542    Double_t pos[3],cov[6],chi2perNDF;
1543    vertexESD->GetXYZ(pos); // position
1544    vertexESD->GetCovMatrix(cov); //covariance matrix
1545    chi2perNDF = vertexESD->GetChi2toNDF();
1546    delete vertexESD; vertexESD=NULL;
1547    
1548    pvtx->SetPosition(pos[0],pos[1],pos[2]);
1549    pvtx->SetChi2perNDF(chi2perNDF);
1550    pvtx->SetCovMatrix(cov);
1551
1552    return kTRUE;
1553 }