Implementing better the kSigma2 approach - dEdx(measured)/dEdx(theory)
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysisBase.cxx
1 /**************************************************************************
2  * Author: Panos Christakoglou.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 /* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
15
16 //-----------------------------------------------------------------
17 //                 AliProtonAnalysisBase class
18 //   This is the class to deal with the proton analysis
19 //   Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21 #include <Riostream.h>
22 #include <TCanvas.h>
23 #include <TLatex.h>
24 #include <TF1.h>
25 #include <TList.h>
26 #include <TH1F.h>
27
28 #include <AliExternalTrackParam.h>
29 #include <AliESDEvent.h>
30 #include <AliPID.h>
31 #include <AliVertexerTracks.h>
32 #include <AliESDpid.h>
33 #include <AliTPCPIDResponse.h>
34 class AliLog;
35 class AliESDVertex;
36
37 #include "AliProtonAnalysisBase.h"
38
39 ClassImp(AliProtonAnalysisBase)
40
41 //____________________________________________________________________//
42 AliProtonAnalysisBase::AliProtonAnalysisBase() : 
43   TObject(),  fProtonAnalysisLevel("ESD"), fAnalysisMC(kFALSE),
44   fTriggerMode(kMB2), kUseOnlineTrigger(kFALSE), kUseOfflineTrigger(kFALSE), 
45   fPhysicsSelection(0),
46   fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
47   fAnalysisEtaMode(kFALSE),
48   fVxMax(100.), fVyMax(100.), fVzMax(100.), fMinNumOfContributors(0),
49   fNBinsX(0), fMinX(0), fMaxX(0),
50   fNBinsY(0), fMinY(0), fMaxY(0),
51   fMinTPCClusters(0), fMinITSClusters(0),
52   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
53   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
54   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
55   fMaxDCAXY(0), fMaxDCAXYTPC(0),
56   fMaxDCAZ(0), fMaxDCAZTPC(0),
57   fMaxDCA3D(0), fMaxDCA3DTPC(0),
58   fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
59   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
60   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
61   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
62   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
63   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
64   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
65   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
66   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
67   fMaxConstrainChi2Flag(kFALSE),
68   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
69   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
70   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
71   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
72   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
73   fMinTPCdEdxPointsFlag(kFALSE),
74   fFunctionProbabilityFlag(kFALSE), 
75   fNSigma(0),
76   fElectronFunction(0), fMuonFunction(0),
77   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
78   fDebugMode(kFALSE), fListVertexQA(new TList()) {
79   //Default constructor
80   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
81   for(Int_t i = 0; i < 24; i++) {
82     fdEdxMean[i] = 0.0;
83     fdEdxSigma[i] = 0.0;
84   }
85   fListVertexQA->SetName("fListVertexQA");
86   TH1F *gHistVx = new TH1F("gHistVx",
87                            "Vx distribution;V_{x} [cm];Entries",
88                            100,-5.,5.);
89   gHistVx->SetFillColor(kRed-2);
90   fListVertexQA->Add(gHistVx);
91   TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
92                                    "Vx distribution;V_{x} [cm];Entries",
93                                    100,-5.,5.);
94   fListVertexQA->Add(gHistVxAccepted);
95   TH1F *gHistVy = new TH1F("gHistVy",
96                            "Vy distribution;V_{y} [cm];Entries",
97                            100,-5.,5.);
98   gHistVy->SetFillColor(kRed-2);
99   fListVertexQA->Add(gHistVy);
100   TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
101                                    "Vy distribution;V_{y} [cm];Entries",
102                                    100,-5.,5.);
103   fListVertexQA->Add(gHistVyAccepted);
104   TH1F *gHistVz = new TH1F("gHistVz",
105                            "Vz distribution;V_{z} [cm];Entries",
106                            100,-25.,25.);
107   gHistVz->SetFillColor(kRed-2);
108   fListVertexQA->Add(gHistVz);
109   TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
110                                    "Vz distribution;V_{z} [cm];Entries",
111                                    100,-25.,25.);
112   fListVertexQA->Add(gHistVzAccepted);
113
114   TH1F *gHistNumberOfContributors = new TH1F("gHistNumberOfContributors",
115                                              "Number of contributors;N_{contr.};Entries",
116                                              100,0.,100.);
117   fListVertexQA->Add(gHistNumberOfContributors);
118
119
120 }
121
122 //____________________________________________________________________//
123 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
124   //Default destructor
125   if(fElectronFunction) delete fElectronFunction;
126   if(fMuonFunction) delete fMuonFunction;
127   if(fPionFunction) delete fPionFunction;
128   if(fKaonFunction) delete fKaonFunction;
129   if(fProtonFunction) delete fProtonFunction;
130   if(fListVertexQA) delete fListVertexQA;
131 }
132
133 //____________________________________________________________________//
134 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
135   //Return the a priori probs
136   Double_t partFrac=0;
137   if(fFunctionProbabilityFlag) {
138     if(i == 0) partFrac = fElectronFunction->Eval(p);
139     if(i == 1) partFrac = fMuonFunction->Eval(p);
140     if(i == 2) partFrac = fPionFunction->Eval(p);
141     if(i == 3) partFrac = fKaonFunction->Eval(p);
142     if(i == 4) partFrac = fProtonFunction->Eval(p);
143   }
144   else partFrac = fPartFrac[i];
145
146   return partFrac;
147 }
148
149 //____________________________________________________________________//
150 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
151   // Checks if the track is outside the analyzed y-Pt phase space
152   Double_t gP = 0.0, gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
153   Double_t eta = 0.0;
154
155   if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid)) {
156     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
157     if(!tpcTrack) {
158       gP = 0.0; gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
159     }
160     else {
161       gP = tpcTrack->P();
162       gPt = tpcTrack->Pt();
163       gPx = tpcTrack->Px();
164       gPy = tpcTrack->Py();
165       gPz = tpcTrack->Pz();
166       eta = tpcTrack->Eta();
167     }
168   }//standalone TPC or Hybrid TPC approaches
169   else {
170     gP = track->P();
171     gPt = track->Pt();
172     gPx = track->Px();
173     gPy = track->Py();
174     gPz = track->Pz();
175     eta = track->Eta();
176   }
177   
178   if((gPt < fMinY) || (gPt > fMaxY)) {
179       if(fDebugMode)
180         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
181       return kFALSE;
182   }
183   if((gP < fMinY) || (gP > fMaxY)) {
184       if(fDebugMode)
185         Printf("IsInPhaseSpace: Track rejected because it has a P value of %lf (accepted interval: %lf - %lf)",gP,fMinY,fMaxY);
186       return kFALSE;
187   }
188   if(fAnalysisEtaMode) {
189     if((eta < fMinX) || (eta > fMaxX)) {
190       if(fDebugMode)
191         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
192       return kFALSE;
193     }
194   }
195   else {
196     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
197       if(fDebugMode)
198         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
199       return kFALSE;
200     }
201   }
202
203   return kTRUE;
204 }
205
206 //____________________________________________________________________//
207 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDtrack* track) {
208   // Checks if the track is excluded from the cuts
209   Int_t  fIdxInt[200];
210   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
211   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
212
213   Float_t chi2PerClusterITS = -1;
214   if (nClustersITS!=0)
215     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
216   Float_t chi2PerClusterTPC = -1;
217   if (nClustersTPC!=0)
218     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
219
220   Double_t extCov[15];
221   track->GetExternalCovariance(extCov);
222
223   if(fPointOnITSLayer1Flag) {
224     if(!track->HasPointOnITSLayer(0)) {
225       if(fDebugMode)
226         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
227       return kFALSE;
228     }
229   }
230   if(fPointOnITSLayer2Flag) {
231     if(!track->HasPointOnITSLayer(1)) {
232       if(fDebugMode)
233         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
234       return kFALSE;
235     }
236   }
237   if(fPointOnITSLayer3Flag) {
238     if(!track->HasPointOnITSLayer(2)) {
239       if(fDebugMode)
240         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
241       return kFALSE;
242     }
243   }
244   if(fPointOnITSLayer4Flag) {
245     if(!track->HasPointOnITSLayer(3)) {
246       if(fDebugMode)
247         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
248       return kFALSE;
249     }
250   }
251   if(fPointOnITSLayer5Flag) {
252     if(!track->HasPointOnITSLayer(4)) {
253       if(fDebugMode)
254         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
255       return kFALSE;
256     }
257   }
258   if(fPointOnITSLayer6Flag) {
259     if(!track->HasPointOnITSLayer(5)) {
260       if(fDebugMode)
261         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
262       return kFALSE;
263     }
264   }
265   if(fMinITSClustersFlag) {
266     if(nClustersITS < fMinITSClusters) {
267       if(fDebugMode)
268         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
269       return kFALSE;
270     }
271   }
272   if(fMaxChi2PerITSClusterFlag) {
273     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
274       if(fDebugMode)
275         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
276       return kFALSE; 
277     }
278   }
279   if(fMinTPCClustersFlag) {
280     if(nClustersTPC < fMinTPCClusters) {
281       if(fDebugMode)
282         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
283       return kFALSE;
284     }
285   }
286   if(fMaxChi2PerTPCClusterFlag) {
287     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
288       if(fDebugMode)
289         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
290       return kFALSE; 
291     }
292   }
293   if(fMaxCov11Flag) {
294     if(extCov[0] > fMaxCov11) {
295       if(fDebugMode)
296         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
297       return kFALSE;
298     }
299   }
300   if(fMaxCov22Flag) {
301     if(extCov[2] > fMaxCov22) {
302       if(fDebugMode)
303         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
304       return kFALSE;
305     }
306   }
307   if(fMaxCov33Flag) {
308     if(extCov[5] > fMaxCov33) {
309       if(fDebugMode)
310         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
311       return kFALSE;
312     }
313   }
314   if(fMaxCov44Flag) {
315     if(extCov[9] > fMaxCov44) {
316       if(fDebugMode)
317         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
318       return kFALSE;
319     }
320   }
321   if(fMaxCov55Flag) {
322     if(extCov[14] > fMaxCov55) {
323       if(fDebugMode)
324         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
325       return kFALSE;
326     }
327   }
328   if(fMinTPCdEdxPointsFlag) {
329     if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
330       if(fDebugMode)
331         Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
332       return kFALSE;
333     }
334   }
335   if(fITSRefitFlag) {
336     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
337       if(fDebugMode)
338         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
339       return kFALSE;
340     }
341   }
342   if(fTPCRefitFlag) {
343     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
344       if(fDebugMode)
345         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
346       return kFALSE;
347     }
348   }
349   if(fESDpidFlag) {
350     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
351       if(fDebugMode)
352         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
353       return kFALSE;
354     }
355   }
356   if(fTPCpidFlag) {
357     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
358       if(fDebugMode)
359         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
360       return kFALSE;
361     }
362   }
363   if(fTOFpidFlag) {
364     if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
365       if(fDebugMode)
366         Printf("IsAccepted: Track rejected because it has no TOF pid flag");
367       return kFALSE;
368     }
369   }
370
371   return kTRUE;
372 }
373
374 //____________________________________________________________________//
375 Bool_t AliProtonAnalysisBase::IsPrimary(AliESDEvent *esd,
376                                         const AliESDVertex *vertex, 
377                                         AliESDtrack* track) {
378   // Checks if the track is a primary-like candidate
379   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
380   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
381   Double_t dca3D = 0.0;
382   
383   if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
384     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
385     if(!tpcTrack) {
386       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
387       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
388       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
389     }
390     else {
391       gPt = tpcTrack->Pt();
392       gPx = tpcTrack->Px();
393       gPy = tpcTrack->Py();
394       gPz = tpcTrack->Pz();
395       tpcTrack->PropagateToDCA(vertex,
396                                esd->GetMagneticField(),
397                                100.,dca,cov);
398     }
399   }//standalone TPC or hybrid TPC approaches
400   else {
401     gPt = track->Pt();
402     gPx = track->Px();
403     gPy = track->Py();
404     gPz = track->Pz();
405     track->PropagateToDCA(vertex,
406                           esd->GetMagneticField(),
407                           100.,dca,cov);
408   }
409   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
410                       TMath::Power(dca[1],2));
411      
412   if(fMaxSigmaToVertexFlag) {
413     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
414       if(fDebugMode)
415         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
416       return kFALSE;
417     }
418   }
419   if(fMaxSigmaToVertexTPCFlag) {
420     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
421       if(fDebugMode)
422         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
423       return kFALSE;
424     }
425   }
426   if(fMaxDCAXYFlag) { 
427     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
428       if(fDebugMode)
429         Printf("IsPrimary: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
430       return kFALSE;
431     }
432   }
433   if(fMaxDCAXYTPCFlag) { 
434     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
435       if(fDebugMode)
436         Printf("IsPrimary: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
437       return kFALSE;
438     }
439   }
440   if(fMaxDCAZFlag) { 
441     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
442       if(fDebugMode)
443         Printf("IsPrimary: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
444       return kFALSE;
445     }
446   }
447   if(fMaxDCAZTPCFlag) { 
448     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
449       if(fDebugMode)
450         Printf("IsPrimary: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
451       return kFALSE;
452     }
453   }
454   if(fMaxDCA3DFlag) { 
455     if(TMath::Abs(dca3D) > fMaxDCA3D) {
456       if(fDebugMode)
457         Printf("IsPrimary: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
458       return kFALSE;
459     }
460   }
461   if(fMaxDCA3DTPCFlag) { 
462     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
463       if(fDebugMode)
464         Printf("IsPrimary: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
465       return kFALSE;
466     }
467   }
468   if(fMaxConstrainChi2Flag) {
469     if(track->GetConstrainedChi2() > 0) 
470       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
471       if(fDebugMode)
472         Printf("IsPrimary: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
473       return kFALSE;
474       }
475   }
476
477   return kTRUE;
478 }
479
480 //____________________________________________________________________//
481 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
482   // Calculates the number of sigma to the vertex.
483   Float_t b[2];
484   Float_t bRes[2];
485   Float_t bCov[3];
486   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
487     esdTrack->GetImpactParametersTPC(b,bCov);
488   else
489     esdTrack->GetImpactParameters(b,bCov);
490   
491   if (bCov[0]<=0 || bCov[2]<=0) {
492     //AliDebug(1, "Estimated b resolution lower or equal zero!");
493     bCov[0]=0; bCov[2]=0;
494   }
495   bRes[0] = TMath::Sqrt(bCov[0]);
496   bRes[1] = TMath::Sqrt(bCov[2]);
497   
498   if (bRes[0] == 0 || bRes[1] ==0) return -1;
499   
500   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
501   
502   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
503   
504   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
505   
506   return d;
507 }
508
509 //____________________________________________________________________//
510 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx, 
511                                          Double_t gPy, 
512                                          Double_t gPz) const {
513   //returns the rapidity of the proton - to be removed
514   Double_t fMass = 9.38270000000000048e-01;
515   
516   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) + 
517                            TMath::Power(gPy,2) + 
518                            TMath::Power(gPz,2));
519   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
520   Double_t y = -999;
521   if(energy != gPz) 
522     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
523
524   return y;
525 }
526
527 //________________________________________________________________________
528 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd, 
529                                                      AnalysisMode mode,
530                                                      Double_t gVxMax,
531                                                      Double_t gVyMax,
532                                                      Double_t gVzMax) {
533   // Get the vertex from the ESD and returns it if the vertex is valid
534   // Second argument decides which vertex is used (this selects
535   // also the quality criteria that are applied)
536   const AliESDVertex* vertex = 0;
537   if(mode == kHybrid)
538     vertex = esd->GetPrimaryVertexSPD();
539   else if(mode == kTPC){
540     Double_t kBz = esd->GetMagneticField();
541     AliVertexerTracks vertexer(kBz);
542     vertexer.SetTPCMode();
543     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
544     esd->SetPrimaryVertexTPC(vTPC);
545     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
546       AliESDtrack *t = esd->GetTrack(i);
547       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
548     }
549     delete vTPC;
550     vertex = esd->GetPrimaryVertexTPC();
551   }
552   else if(mode == kGlobal)
553     vertex = esd->GetPrimaryVertex();
554   else
555     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
556   
557   if(!vertex) {
558     if(fDebugMode)
559       Printf("GetVertex: Event rejected because there is no valid vertex object");
560     return 0;
561   }
562   
563   // check Ncontributors
564   if(vertex->GetNContributors() <= 0) {
565     if(fDebugMode)
566       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
567     return 0;
568   }
569   
570   // check resolution
571   Double_t zRes = vertex->GetZRes();
572   if(zRes == 0) {
573     if(fDebugMode)
574       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
575     return 0;
576   }
577   ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
578   ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
579   ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
580
581   //check position
582   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
583     if(fDebugMode)
584       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
585     return 0;
586   }
587   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
588     if(fDebugMode)
589       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
590     return 0;
591   }
592   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
593     if(fDebugMode)
594       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
595     return 0;
596   }
597   ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
598   ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
599   ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
600   ((TH1F *)(fListVertexQA->At(6)))->Fill(vertex->GetNContributors());
601
602   //check number of contributors
603   if(fMinNumOfContributors > 0) {
604     if(fMinNumOfContributors > vertex->GetNContributors()) {
605       if(fDebugMode)
606         Printf("GetVertex: Event rejected because it has %d number of contributors (requested minimum: %d)",vertex->GetNContributors(),fMinNumOfContributors);
607       
608       return 0;
609     }
610   }
611   
612   return vertex;
613 }
614
615 //________________________________________________________________________
616 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd, 
617                                                TriggerMode trigger) {
618   // check if the event was triggered
619   ULong64_t triggerMask = esd->GetTriggerMask();
620   TString firedTriggerClass = esd->GetFiredTriggerClasses();
621
622   if(fAnalysisMC) {
623     // definitions from p-p.cfg
624     ULong64_t spdFO = (1 << 14);
625     ULong64_t v0left = (1 << 11);
626     ULong64_t v0right = (1 << 12);
627     
628     switch (trigger) {
629     case kMB1: {
630       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
631         return kTRUE;
632       break;
633     }
634     case kMB2: {
635       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
636         return kTRUE;
637       break;
638     }
639     case kSPDFASTOR: {
640       if (triggerMask & spdFO)
641         return kTRUE;
642       break;
643     }
644     }//switch
645   }
646   else {
647     if(kUseOnlineTrigger) {
648       if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
649         return kTRUE;
650     }
651     else if(!kUseOnlineTrigger)
652       return kTRUE;
653   }
654
655   return kFALSE;
656 }
657
658 //________________________________________________________________________
659 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
660   // return the list of cuts and their cut values for reference
661   TLatex l;
662   l.SetTextAlign(12);
663   l.SetTextSize(0.04);
664
665   TString listOfCuts;
666
667   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
668   c->SetFillColor(10); c->SetHighLightColor(41);
669   c->Divide(3,2);
670
671   c->cd(1)->SetFillColor(10);
672   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
673
674   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
675   l.DrawLatex(0.1,0.82,listOfCuts.Data());
676   listOfCuts = "Analysis mode: ";
677   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone"; 
678   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC"; 
679   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking"; 
680   l.DrawLatex(0.1,0.74,listOfCuts.Data());
681   listOfCuts = "Trigger mode: "; 
682   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1"; 
683   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2"; 
684   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)"; 
685   l.DrawLatex(0.1,0.66,listOfCuts.Data());
686   listOfCuts = "PID mode: "; 
687   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
688   if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})"; 
689   if(fProtonPIDMode == kSigma1) {
690     listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
691     listOfCuts += " #sigma";
692   }
693   if(fProtonPIDMode == kSigma2) {
694     listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
695     listOfCuts += " #sigma";
696   }
697   l.DrawLatex(0.1,0.58,listOfCuts.Data());
698   listOfCuts = "Accepted vertex diamond: "; 
699   l.DrawLatex(0.1,0.5,listOfCuts.Data());
700   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
701   l.DrawLatex(0.6,0.5,listOfCuts.Data());
702   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
703   l.DrawLatex(0.6,0.4,listOfCuts.Data());
704   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
705   l.DrawLatex(0.6,0.3,listOfCuts.Data());
706   listOfCuts = "Phase space: "; 
707   l.DrawLatex(0.1,0.2,listOfCuts.Data());
708   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";  
709   else listOfCuts = "|y| < ";
710   listOfCuts += TMath::Abs(fMinX); 
711   l.DrawLatex(0.35,0.2,listOfCuts.Data());
712   listOfCuts = "N_{bins} = ";
713   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
714   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
715   l.DrawLatex(0.35,0.15,listOfCuts.Data());
716   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
717   listOfCuts += fMaxY; listOfCuts += "GeV/c"; 
718   l.DrawLatex(0.35,0.1,listOfCuts.Data());
719   listOfCuts = "N_{bins} = ";
720   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
721   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
722   l.DrawLatex(0.35,0.05,listOfCuts.Data());
723
724   c->cd(2)->SetFillColor(2);
725   l.DrawLatex(0.3,0.9,"ITS related cuts");
726   listOfCuts = "Request a cluster on SPD1: "; 
727   listOfCuts += fPointOnITSLayer1Flag;
728   l.DrawLatex(0.1,0.8,listOfCuts.Data());
729   listOfCuts = "Request a cluster on SPD2: "; 
730   listOfCuts += fPointOnITSLayer2Flag;
731   l.DrawLatex(0.1,0.7,listOfCuts.Data());
732   listOfCuts = "Request a cluster on SDD1: "; 
733   listOfCuts += fPointOnITSLayer3Flag;
734   l.DrawLatex(0.1,0.6,listOfCuts.Data());
735   listOfCuts = "Request a cluster on SDD2: "; 
736   listOfCuts += fPointOnITSLayer4Flag;
737   l.DrawLatex(0.1,0.5,listOfCuts.Data());
738   listOfCuts = "Request a cluster on SSD1: "; 
739   listOfCuts += fPointOnITSLayer5Flag;
740   l.DrawLatex(0.1,0.4,listOfCuts.Data());
741   listOfCuts = "Request a cluster on SSD2: "; 
742   listOfCuts += fPointOnITSLayer6Flag; 
743   l.DrawLatex(0.1,0.3,listOfCuts.Data());  
744   listOfCuts = "Minimum number of ITS clusters: ";
745   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
746   else listOfCuts += "Not used";
747   l.DrawLatex(0.1,0.2,listOfCuts.Data());
748   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
749   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster; 
750   else listOfCuts += "Not used";
751   l.DrawLatex(0.1,0.1,listOfCuts.Data());
752
753   c->cd(3)->SetFillColor(3);
754   l.DrawLatex(0.3,0.9,"TPC related cuts");
755   listOfCuts = "Minimum number of TPC clusters: ";
756   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
757   else listOfCuts += "Not used";
758   l.DrawLatex(0.1,0.7,listOfCuts.Data());
759   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
760   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
761   else listOfCuts += "Not used";
762   l.DrawLatex(0.1,0.5,listOfCuts.Data());
763   listOfCuts = "Minimum number of TPC points for the dE/dx: ";
764   if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
765   else listOfCuts += "Not used";
766   l.DrawLatex(0.1,0.3,listOfCuts.Data());
767
768   c->cd(4)->SetFillColor(4);
769   l.DrawLatex(0.3,0.9,"Tracking related cuts");
770   listOfCuts = "Maximum cov11: ";
771   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
772   else listOfCuts += "Not used";
773   l.DrawLatex(0.1,0.75,listOfCuts.Data());
774   listOfCuts = "Maximum cov22: ";
775   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
776   else listOfCuts += "Not used";
777   l.DrawLatex(0.1,0.6,listOfCuts.Data());
778   listOfCuts = "Maximum cov33: ";
779   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
780   else listOfCuts += "Not used";
781   l.DrawLatex(0.1,0.45,listOfCuts.Data());
782   listOfCuts = "Maximum cov44: ";
783   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
784   else listOfCuts += "Not used";
785   l.DrawLatex(0.1,0.3,listOfCuts.Data());
786   listOfCuts = "Maximum cov55: ";
787   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
788   else listOfCuts += "Not used";
789   l.DrawLatex(0.1,0.15,listOfCuts.Data());
790
791   c->cd(5)->SetFillColor(5);
792   l.DrawLatex(0.3,0.9,"DCA related cuts");
793   listOfCuts = "Maximum sigma to vertex: ";
794   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
795   else listOfCuts += "Not used";
796   l.DrawLatex(0.1,0.8,listOfCuts.Data());
797   listOfCuts = "Maximum sigma to vertex (TPC): ";
798   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
799   else listOfCuts += "Not used";
800   l.DrawLatex(0.1,0.72,listOfCuts.Data());
801   listOfCuts = "Maximum DCA in xy: ";
802   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
803   else listOfCuts += "Not used";
804   l.DrawLatex(0.1,0.64,listOfCuts.Data());
805   listOfCuts = "Maximum DCA in z: ";
806   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
807   else listOfCuts += "Not used";
808   l.DrawLatex(0.1,0.56,listOfCuts.Data());
809   listOfCuts = "Maximum DCA in 3D: ";
810   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
811   else listOfCuts += "Not used";
812   l.DrawLatex(0.1,0.48,listOfCuts.Data());
813   listOfCuts = "Maximum DCA in xy (TPC): ";
814   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
815   else listOfCuts += "Not used";
816   l.DrawLatex(0.1,0.4,listOfCuts.Data());
817   listOfCuts = "Maximum DCA in z (TPC): ";
818   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
819   else listOfCuts += "Not used";
820   l.DrawLatex(0.1,0.32,listOfCuts.Data());
821   listOfCuts = "Maximum DCA in 3D (TPC): ";
822   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
823   else listOfCuts += "Not used";
824   l.DrawLatex(0.1,0.24,listOfCuts.Data());
825   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
826   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
827   else listOfCuts += "Not used";
828   l.DrawLatex(0.1,0.16,listOfCuts.Data());
829
830   c->cd(6)->SetFillColor(6);
831   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
832   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
833   l.DrawLatex(0.1,0.7,listOfCuts.Data());
834   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag; 
835   l.DrawLatex(0.1,0.5,listOfCuts.Data());
836   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
837   l.DrawLatex(0.1,0.3,listOfCuts.Data());
838   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
839   l.DrawLatex(0.1,0.1,listOfCuts.Data());
840
841   return c;
842 }
843
844 //________________________________________________________________________
845 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
846   //Function that checks if a track is a proton
847   Double_t probability[5];
848   Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
849   Long64_t fParticleType = 0;
850  
851   //Bayesian approach for the PID
852   if(fProtonPIDMode == kBayesian) {
853     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
854       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
855       if(tpcTrack) {
856         gPt = tpcTrack->Pt();
857         gP = tpcTrack->P();
858         track->GetTPCpid(probability);
859       }
860     }//TPC standalone or Hybrid TPC
861     else if(fProtonAnalysisMode == kGlobal) {
862       gPt = track->Pt();
863       gP = track->P();
864       track->GetESDpid(probability);
865     }//Global tracking    
866     
867     Double_t rcc = 0.0;
868     for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
869       rcc += probability[i]*GetParticleFraction(i,gP);
870     if(rcc != 0.0) {
871       Double_t w[5];
872       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
873         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
874       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
875     }
876     if(fParticleType == 4)
877       return kTRUE;
878   }
879   //Ratio of the measured over the theoretical dE/dx a la STAR
880   else if(fProtonPIDMode == kRatio) {
881     Printf("The kRatio mode is not implemented yet!!!");
882     return kFALSE;
883   }
884   //Definition of an N-sigma area around the dE/dx vs P band
885   else if(fProtonPIDMode == kSigma1) {
886     Double_t fAlephParameters[5];
887     if(fAnalysisMC) {
888       fAlephParameters[0] = 2.15898e+00/50.;
889       fAlephParameters[1] = 1.75295e+01;
890       fAlephParameters[2] = 3.40030e-09;
891       fAlephParameters[3] = 1.96178e+00;
892       fAlephParameters[4] = 3.91720e+00;
893     }
894     else {
895       fAlephParameters[0] = 0.0283086;
896       fAlephParameters[1] = 2.63394e+01;
897       fAlephParameters[2] = 5.04114e-11;
898       fAlephParameters[3] = 2.12543e+00;
899       fAlephParameters[4] = 4.88663e+00;
900     }
901
902     Double_t nsigma = 100.0;
903     AliESDpid *fESDpid = new AliESDpid(); 
904     fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
905     
906     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
907     if(tpcTrack)
908       nsigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton));
909   
910     if(nsigma <= fNSigma) 
911       return kTRUE;
912   }//kSigma1 PID method
913   //Another definition of an N-sigma area around the dE/dx vs P band
914   else if(fProtonPIDMode == kSigma2) {
915     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
916     if(tpcTrack) {
917       gPt = tpcTrack->Pt();
918       gP = tpcTrack->P();
919       gEta = tpcTrack->Eta();
920     }
921     Double_t fAlephParameters[5];
922     if(fAnalysisMC) {
923       fAlephParameters[0] = 2.15898e+00/50.;
924       fAlephParameters[1] = 1.75295e+01;
925       fAlephParameters[2] = 3.40030e-09;
926       fAlephParameters[3] = 1.96178e+00;
927       fAlephParameters[4] = 3.91720e+00;
928     }
929     else {
930       fAlephParameters[0] = 0.0283086;
931       fAlephParameters[1] = 2.63394e+01;
932       fAlephParameters[2] = 5.04114e-11;
933       fAlephParameters[3] = 2.12543e+00;
934       fAlephParameters[4] = 4.88663e+00;
935     }
936
937     AliESDpid *fESDpid = new AliESDpid(); 
938     AliTPCPIDResponse tpcResponse = fESDpid->GetTPCResponse(); 
939     tpcResponse.SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
940     Double_t normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse.GetExpectedSignal(gP,AliPID::kProton));
941
942     if(normalizeddEdx >= -0.15)
943       return kTRUE;
944   }
945
946   return kFALSE;
947 }
948
949 //________________________________________________________________________
950 void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
951   // This function is used in case the kSigma1 or kSigma2 PID mode is selected
952   // It takes as an argument the name of the ascii file (for the time being) 
953   // that is generated as a prior process.
954   // This ascii file has three columns: The min. P value (bins of 50MeV/c) 
955   // the mean and the sigma of the dE/dx distributions for protons coming 
956   // from a gaussian fit.
957   ifstream in;
958   in.open(filename);
959
960   Double_t gPtMin = 0.0;
961   Int_t iCounter = 0;
962   while(in.good()) {
963     in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
964     if(fDebugMode)
965       Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
966     iCounter += 1;
967   }
968 }
969
970 //________________________________________________________________________
971 Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
972   // This is the Bethe-Bloch function normalised to 1 at the minimum
973   // We renormalize it based on the MC information
974   // WARNING: To be revised soon!!!
975   // This is just a temporary fix
976   Double_t normalization = 49.2;
977   Double_t bg2=bg*bg;
978   Double_t beta2 = bg2/(1.+ bg2);
979   Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
980   //
981   const Float_t kmeanCorrection =0.1;
982   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
983   bb *= meanCorrection;
984
985   return normalization*bb; 
986 }
987
988