o small fix
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxBaseUtils.cxx
CommitLineData
6951a056 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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//
17// TRD dEdx base utils
18// xx
19// xx
20// xx
21// xx
22//
23// Xianguo Lu
24// lu@physi.uni-heidelberg.de
25// Xianguo.Lu@cern.ch
26//
27//
28
29#include "TF1.h"
30#include "TFile.h"
31#include "TH1D.h"
32#include "TH2D.h"
33#include "THnSparse.h"
34#include "TMath.h"
35#include "TMatrixD.h"
36#include "TMinuit.h"
37#include "TObjArray.h"
38#include "TRandom3.h"
39#include "TStopwatch.h"
40#include "TVectorD.h"
41
42#include "TTreeStream.h"
43
6951a056 44#include "AliESDEvent.h"
45#include "AliESDfriendTrack.h"
46#include "AliESDtrack.h"
6951a056 47#include "AliTRDtrackV1.h"
48
49#include "AliTRDdEdxBaseUtils.h"
50
6fa94550 51#define EPSILON 1e-8
6951a056 52
687aa844 53Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.55;
6951a056 54Double_t AliTRDdEdxBaseUtils::fgQ1Frac = 0.5;
55Double_t AliTRDdEdxBaseUtils::fgTimeBinCountCut = 0.0;
56Int_t AliTRDdEdxBaseUtils::fgCalibTPCnclsCut = 70;
57Bool_t AliTRDdEdxBaseUtils::fgExBOn = kTRUE;
58Bool_t AliTRDdEdxBaseUtils::fgPadGainOn = kTRUE;
687aa844 59Double_t AliTRDdEdxBaseUtils::fgQScale = 50;
6951a056 60
61//===================================================================================
62// Math and Histogram
63//===================================================================================
6856b6a8 64void AliTRDdEdxBaseUtils::BinLogX(TAxis *axis)
65{
66 //
67 // Method for the correct logarithmic binning of histograms
68 // copied and modified from AliTPCcalibBase
69
70 const Int_t bins = axis->GetNbins();
71
72 const Double_t from = axis->GetXmin();
73 const Double_t to = axis->GetXmax();
ca66da8f 74 if (from<EPSILON){
396eba11 75 //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from); exit(1);
76 return;
ca66da8f 77 }
78
6856b6a8 79 Double_t *new_bins = new Double_t[bins + 1];
80
81 new_bins[0] = from;
82 const Double_t factor = pow(to/from, 1./bins);
83
84 for (int i = 1; i <= bins; i++) {
85 new_bins[i] = factor * new_bins[i-1];
86 }
87 axis->Set(bins, new_bins);
88 delete [] new_bins;
89}
90
6951a056 91void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
92{
93 //
94 //counts of hh is sorted
95 //
96
97 for(Int_t ii=0; ii<ncut; ii++){
98 cuts[ii] = -999;
99 }
100
101 Int_t nsel = 0;
102 const Int_t nbin = hh->GetNbinsX();
103 Double_t datas[nbin];
104 for(Int_t ii=1; ii<=nbin; ii++){
105 const Double_t res = hh->GetBinContent(ii);
106 if(res<thres){
107 continue;
108 }
109
110 datas[nsel] = res;
111 nsel++;
112 }
113 if(!nsel)
114 return;
115
116 Int_t id[nsel];
117 TMath::Sort(nsel, datas, id, kFALSE);
118
119 for(Int_t ii=0; ii<ncut; ii++){
120 const Double_t icdf = cdfs[ii];
121 if(icdf<0 || icdf>1){
122 printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
123 }
124 cuts[ii] = datas[id[Int_t(icdf*nsel)]];
125 }
126}
127
128Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
129{
130 //
131 //calculate mean (with error) and rms from sum, w2s, nn
132 //if nn=0, mean, error, and rms all = 0
133 //
134
135 Double_t tmean = 0, trms = 0, terr = 0;
136
137 if(nn>EPSILON){
138 tmean = sum/nn;
139
140 const Double_t arg = w2s/nn-tmean*tmean;
141 if(TMath::Abs(arg)<EPSILON){
142 trms = 0;
143 }
144 else{
145 if( arg <0 ){
146 printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
147 }
148
149 trms = TMath::Sqrt(arg);
150 }
151
152 terr = trms/TMath::Sqrt(nn);
153 }
154
155 if(grms){
156 (*grms) = trms;
157 }
158
159 if(gerr){
160 (*gerr) = terr;
161 }
162
163 return tmean;
164}
165
166Double_t AliTRDdEdxBaseUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
167{
168 //
169 //calculate truncated mean
170 //return <x*w>_{low-high according to x}
171 //
172
173 /*
174 //test->
175 for(Int_t ii=0; ii<nx; ii++){
176 printf("test %d/%d %f\n", ii, nx, xdata[ii]);
177 }
178 //<--test
179 */
180
181 Int_t index[nx];
182 TMath::Sort(nx, xdata, index, kFALSE);
183
184 Int_t nused = 0;
185 Double_t sum = 0;
186 Double_t w2s = 0;
187 const Int_t istart = Int_t (nx*lowfrac);
188 const Int_t istop = Int_t (nx*highfrac);
189
190 //=,< correct, because when low=0, high=1 it is correct
191 for(Int_t ii=istart; ii<istop; ii++){
192 Double_t weight = 1;
193 if(wws){
194 weight = wws[index[ii]];
195 }
196 const Double_t sx = xdata[index[ii]]*weight;
197
198 sum += sx;
199 w2s += sx*sx;
200
201 nused++;
202 //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
203
204 }
205
206 return GetMeanRMS(nused, sum, w2s, grms, gerr);
207}
208
209Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
210{
211 //
212 //do truncation on histogram
213 //
214 //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
215
216 //with under-/over-flow
217 Double_t npreTrunc = 0;
218 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
219 const Double_t be = hh->GetBinError(itmp);
220 const Double_t bc = hh->GetBinContent(itmp);
221 if(be<EPSILON){
222 if(bc>EPSILON){
223 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
224 }
225 continue;
226 }
227 npreTrunc += bc*bc/be/be;
228 }
229
230 const Double_t nstart = npreTrunc*lowfrac;
231 const Double_t nstop = npreTrunc*highfrac;
232
233 //with Double_t this should also handle normalized hist
234 Double_t ntot = 0;
235 Double_t nused = 0;
236 Double_t sum = 0;
237 Double_t w2s = 0;
238 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
239 const Double_t be = hh->GetBinError(itmp);
240 const Double_t bc = hh->GetBinContent(itmp);
241 if(be<EPSILON){
242 if(bc>EPSILON){
243 printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
244 }
245 continue;
246 }
247 const Double_t weight = bc*bc/be/be;
248 ntot+=weight;
249 //<= correct, because when high=1, nstop has to be included
250 if(ntot>nstart && ntot<=nstop){
251
252 const Double_t val = hh->GetBinCenter(itmp);
253 sum += weight*val;
254 w2s += weight*val*val;
255
256 nused += weight;
257
258 //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
259 }
260 else if(ntot>nstop){
261 if(itmp>=hh->GetNbinsX()){
262 printf("AliTRDdEdxBaseUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);
263 nused = 0;
264 w2s = sum = 0;
265 }
266 break;
267 }
268 }
269
270 return GetMeanRMS(nused, sum, w2s, grms, gerr);
271}
272
273void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
274{
275 //
276 //fit slices of hh using truncation
277 //
278
279 const Int_t x0 = hh->GetXaxis()->GetFirst();
280 const Int_t x1 = hh->GetXaxis()->GetLast();
281 const Int_t y0 = hh->GetYaxis()->GetFirst();
282 const Int_t y1 = hh->GetYaxis()->GetLast();
283
284 const Int_t nx = hh->GetNbinsX();
285 const Int_t ny = hh->GetNbinsY();
286 const Double_t xmin = hh->GetXaxis()->GetXmin();
287 const Double_t xmax = hh->GetXaxis()->GetXmax();
288 const Double_t ymin = hh->GetYaxis()->GetXmin();
289 const Double_t ymax = hh->GetYaxis()->GetXmax();
290
291 hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
292 hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
293 hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
294 hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
295
687aa844 296 Double_t newbins[ny+1];
297 for(Int_t ii=1; ii<=ny+1; ii++){
298 const Double_t lowx= hh->GetYaxis()->GetBinLowEdge(ii);
299 newbins[ii-1]=lowx;
300 }
301
6951a056 302 for(Int_t ix=x0; ix<=x1; ix++){
303 //to speed up
304 const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
305 if(rawcount<EPSILON){
306 continue;
307 }
308
309 TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
687aa844 310 htmp->GetXaxis()->Set(ny, newbins);
6951a056 311 Double_t ntot = 0;
312 for(Int_t iy=y0; iy<=y1; iy++){
313 const Double_t be = hh->GetBinError(ix,iy);
314 const Double_t bc = hh->GetBinContent(ix, iy);
315
316 if(be<EPSILON){
317 if(bc>EPSILON){
318 printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
319 }
320 continue;
321 }
322
323 htmp->SetBinContent(iy, bc);
324 htmp->SetBinError(iy, be);
325
326 ntot += (bc/be)*(bc/be);
327
328 //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
329 }
330
331 hnor->SetBinContent(ix, ntot);
332 hnor->SetBinError( ix, 0);
333
334 if(ntot<thres || htmp->GetRMS()<EPSILON){
335 delete htmp;
336 continue;
337 }
338
339 //test htmp->Draw();
340 Double_t trms = -999, terr = -999;
341 const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
342
343 hmpv->SetBinContent(ix, tmean);
344 hmpv->SetBinError( ix, terr);
345
346 hwid->SetBinContent(ix, trms);
347 hwid->SetBinError( ix, 0);
348
349 hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
350 hres->SetBinError( ix, 0);
351
352 delete htmp;
353 }
354
355 TH1 *hhs[]={hnor, hmpv, hwid, hres};
356 const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
357 const Int_t nh = sizeof(hhs)/sizeof(TH1*);
358 for(Int_t ii=0; ii<nh; ii++){
359 hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
360 hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
361 hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
362 hhs[ii]->SetTitle(hh->GetTitle());
363 }
364}
365
366//===================================================================================
367// TRD Analysis Fast Tool
368//===================================================================================
369
370Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
371{
372 //
373 //number of trd tracklet in one esd event
374 //
375 const Int_t ntrk0 = esd->GetNumberOfTracks();
376 Int_t ntrdv1=0;
377 for(Int_t ii=0; ii<ntrk0; ii++){
378 ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
379 }
380 return ntrdv1;
381}
382
383AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
384{
385 //
386 //Get TRD friend track
387 //
388
389 AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
390 if(!friendtrk){
391 //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
392 return 0x0;
393 }
394
395 TObject *calibObject=0x0;
396 AliTRDtrackV1 * trdtrack=0x0;
397 for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
398 if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
399 break;
400 }
401
402 return trdtrack;
403}
404
405Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
406{
407 //
408 // to check if all tracklets are in the same stack, useful in cosmic
409 //
410
411 TVectorD secs(18), stks(5);
412
413 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
414 AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
415 if(!tracklet)
416 continue;
417
418 const Int_t det = tracklet->GetDetector();
419 const Int_t isector = AliTRDgeometry::GetSector(det);
420 const Int_t istack = AliTRDgeometry::GetStack(det);
421
422 secs[isector] = 1;
423 stks[istack] = 1;
424 }
425
426 if(secs.Sum()!=1 || stks.Sum()!=1){
427 return kFALSE;
428 }
429 else
430 return kTRUE;
431}
432
ca66da8f 433Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi)
434{
435 //
436 //redefine phi in pi/2 - pi/2 + 2pi
437 //
438
439 if(phi>0 && phi < TMath::Pi()/2){
440 phi += 2*TMath::Pi();
441 }
442
443 return phi;
444}
445
446Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet)
fe829aa0 447{
448 //
ca66da8f 449 //get dydx
fe829aa0 450 //
451 if(tracklet)
ca66da8f 452 return tracklet->GetYref(1);
fe829aa0 453 else
454 return -999;
455}
456
ca66da8f 457Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet)
458{
459 //
460 //get dzdx
461 //
462 if(tracklet)
463 return tracklet->GetZref(1);
464 else
465 return -999;
466}
467
468Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet)
469{
470 //
471 //return angular normalization factor = dldx
472 //
473 if(tracklet)
474 return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1));
475 else
476 return -999;
477}
478
479AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack)
480{
481 //
482 //as function name
483 //
484
485 AliTRDseedV1 *tracklet = 0x0;
486 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
487 tracklet=trdtrack->GetTracklet(ilayer);
488 if(!tracklet)
489 continue;
490
491 break;
492 }
493
494 return tracklet;
495}
496
fe829aa0 497AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
498{
499 //
500 //get last tracklet
501 //
502 AliTRDseedV1 *tracklet = 0x0;
503
504 for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
505 tracklet=trdtrack->GetTracklet(ilayer);
506 if(!tracklet)
507 continue;
508
509 break;
510 }
511
512 return tracklet;
513}
514
ca66da8f 515void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
6951a056 516{
517 //
518 //as function name
519 //
520 isec = istk = -999;
521 mom = -999;
522
ca66da8f 523 AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack);
524 if(tracklet){
6951a056 525 const Int_t det = tracklet->GetDetector();
526 isec = AliTRDgeometry::GetSector(det);
527 istk = AliTRDgeometry::GetStack(det);
ca66da8f 528
6951a056 529 mom = tracklet->GetMomentum();
6951a056 530 }
531
ca66da8f 532 return;
6951a056 533}
534
535//===================================================================================
536// Detector and Data Constant
537//===================================================================================
538Int_t AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
539{
540 //
541 //gtb = det*Ntb+itb
542 //
543 return gtb/AliTRDseedV1::kNtb;
544}
545
546Int_t AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
547{
548 //
549 //gtb = det*Ntb+itb
550 //
551 return gtb%AliTRDseedV1::kNtb;
552}
553
554Int_t AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
555{
556 //
557 //return sector
558 //
559 return AliTRDgeometry::GetSector(ToDetector(gtb));
560}
561
562Int_t AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
563{
564 //
565 //return stack
566 //
567 return AliTRDgeometry::GetStack(ToDetector(gtb));
568}
569
570Int_t AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
571{
572 //
573 //return layer
574 //
575 return AliTRDgeometry::GetLayer(ToDetector(gtb));
576}
577
ca66da8f 578void AliTRDdEdxBaseUtils::CheckRunB(const TString listrun1kg, const Int_t run, TString & type)
579{
580 //
581 //check run b field
582 //
583 if(listrun1kg.Contains(Form("%d",run))){
584 type+="1kG";
585 }
586 else {
587 type+="5kG";
588 }
589}
590
6951a056 591TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
592{
593 //
594 //return run type
595 //
596
597 TString type;
598 if(run>=121527 && run<= 126460)//LHC10d
ca66da8f 599 type="2010ppLHC10d";
6951a056 600 else if(run>=126461 && run<= 130930)//LHC10e
ca66da8f 601 type="2010ppLHC10e";
6951a056 602 else if(run>=136782 && run <= 139846)//LHC10h
ca66da8f 603 type="2010PbPbLHC10h";
6951a056 604 else if(run>= 143224 && run<= 143237)//2011Feb
ca66da8f 605 type="2011cosmicFeb";
6951a056 606 else if(run>= 150587 && run<= 154930){
ca66da8f 607 type="2011cosmicMayJun";
6951a056 608
ca66da8f 609 CheckRunB("154601 154602 154629 154634 154636 154639 154643", run, type);
610 }
611 else if(run>=173047 && run<=180164){
612 type="2012cosmic";
613
614 CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type);
6951a056 615 }
616 else{
617 type="unknown";
618 }
619
620 type.ToUpper();
621 return type;
622}
623
624void AliTRDdEdxBaseUtils::PrintControl()
625{
626 //
627 //print out control variable
628 //
687aa844 629 printf("\nAliTRDdEdxBaseUtils::PrintControl Q0Frac %.2f, Q1Frac %.2f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d, IsPadGainOn %d, QScale %.2f\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn(), IsPadGainOn(), QScale());
6951a056 630}
631
632//===================================================================================
633// dEdx Parameterization
634//===================================================================================
fe829aa0 635void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
636{
637 //
638 //fit dedx tr from mean
639 //
640
641 TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
642 Double_t par[8];
643 par[0]= 2.397001e-01;
644 par[1]= 1.334697e+00;
645 par[2]= 6.967470e+00;
646 par[3]= 9.055289e-02;
647 par[4]= 9.388760e+00;
648 par[5]= 9.452742e-04;
649 par[6]= -1.866091e+00;
650 par[7]= 1.403545e+00;
651
652 ff->SetParameters(par);
653 hh->Fit(ff,"N");
654
655 for(int ii=0; ii<8; ii++){
656 printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
657 }
658
659 TFile *fout=new TFile("fastfit.root","recreate");
660 hh->Write();
661 ff->Write();
662 fout->Save();
663 fout->Close();
664
665 delete fout;
666 delete ff;
667}
6951a056 668
669Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
670{
671 //
672 //truncated Mean Q_{xx} in TRD
673 //
674
675 Double_t par[8];
fe829aa0 676 Double_t scale = 1;
677
678 //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
679 scale = 0.9257;
680par[0]=8.316476e-01;
681par[1]=1.600907e+00;
682par[2]=7.728447e+00;
683par[3]=6.235622e-02;
684par[4]=1.136209e+01;
685par[5]=-1.495764e-06;
686par[6]=-2.536119e+00;
687par[7]=3.416031e+00;
688
689 return scale*MeandEdxTR(&bg, par);
6951a056 690}
691
692Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
693{
694 //
695 //truncated Mean Q_{xx} in TRD
696 //
697
698 Double_t par[8];
699
700 //03132012161259
701 //opt: PbPbQ0
702par[0]= 1.844912e-01;
703par[1]= 2.509702e+00;
704par[2]= 6.744031e+00;
705par[3]= 7.355123e-02;
706par[4]= 1.166023e+01;
707par[5]= 1.736186e-04;
708par[6]= -1.716063e+00;
709par[7]= 1.611366e+00;
710
711 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
712 return 0.460994*MeandEdxTR(&bg, par);
713}
714
715Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
716{
717 //
718 //truncated Mean 1/(1/Q)_{xx} in TRD
719 //
720
721 Double_t par[8];
722
723 //So 4. Mär 13:30:51 CET 2012
724 //opt: trdppQ1
725 par[0]= 2.434646e-01;
726 par[1]= 1.400211e+00;
727 par[2]= 6.937471e+00;
728 par[3]= 7.758118e-02;
729 par[4]= 1.097372e+01;
730 par[5]= 4.297518e-04;
731 par[6]= -1.806266e+00;
732 par[7]= 1.543811e+00;
733
734 //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
735
736 return 0.418629*MeandEdxTR(&bg, par);
737}
738
739Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
740{
741 //
742 //truncated Mean 1/(1/Q)_{xx} in TRD
743 //
744
745 Double_t par[8];
746
747 //So 4. Mär 13:30:52 CET 2012
748 //opt: trdPbPbQ1
749 par[0]= 2.193660e-01;
750 par[1]= 2.051864e+00;
751 par[2]= 6.825112e+00;
752 par[3]= 6.151693e-02;
753 par[4]= 1.390343e+01;
754 par[5]= 6.010032e-05;
755 par[6]= -1.676324e+00;
756 par[7]= 1.838873e+00;
757
758 //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
759
760 return 0.457988*MeandEdxTR(&bg, par);
761}
762
763Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
764{
765 //
766 //bethe bloch in TPC
767 //
768
769 Double_t par[5];
770 //Mi 15. Feb 14:48:05 CET 2012
771 //train_2012-02-13_1214.12001, tpcsig
772 par[0]= 4.401269;
773 par[1]= 9.725370;
774 par[2]= 0.000178;
775 par[3]= 1.904962;
776 par[4]= 1.426576;
777
778 return MeandEdx(&bg, par);
779}
780
ca66da8f 781Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
6951a056 782{
783 //
784 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
785 //npar = 8 = 3+5
786 //
787 Double_t ptr[4]={0};
788 for(int ii=0; ii<3; ii++){
789 ptr[ii+1]=pin[ii];
790 }
791 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
792}
793
ca66da8f 794Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
6951a056 795{
796 //
797 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
798 //npar = 4
799 //
800
801 const Double_t bg = xx[0];
802 const Double_t gamma = sqrt(1+bg*bg);
803
804 const Double_t p0 = TMath::Abs(par[1]);
805 const Double_t p1 = TMath::Abs(par[2]);
806 const Double_t p2 = TMath::Abs(par[3]);
807
808 const Double_t zz = TMath::Log(gamma);
809 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
810
811 return par[0]+tryield;
812}
813
ca66da8f 814Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
815{
816 return ALEPH(xx, par);
817}
818
819Double_t AliTRDdEdxBaseUtils::ALEPH(const Double_t * xx, const Double_t * par)
6951a056 820{
821 //
822 //ALEPH parametrization for dEdx
823 //npar = 5
824 //
825
826 const Double_t bg = xx[0];
827 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
828
829 const Double_t p0 = TMath::Abs(par[0]);
830 const Double_t p1 = TMath::Abs(par[1]);
ca66da8f 831
832 //after redefining p2 (<0) -> ln P2
833 //const Double_t p2 = par[2];
834
835 //restore back
6951a056 836 const Double_t p2 = TMath::Abs(par[2]);
ca66da8f 837
6951a056 838 const Double_t p3 = TMath::Abs(par[3]);
839 const Double_t p4 = TMath::Abs(par[4]);
840
841 const Double_t aa = TMath::Power(beta, p3);
ca66da8f 842
843 //numerically very bad when p2~1e-15, bg~1e3, p4~5.
6951a056 844 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
845
ca66da8f 846 //+1 to avoid the singularity when bg<1 and p4>>1
847 //very^inf important! with this hesse NOTPOS->OK and no nan or inf in migrad
848 //i.e. numerically very stable
849 //the reason is when bg<1 ln blows up very fast with p4>>1 and nan/inf appears in migrad search
850 //the fitting will adjust the parameters as if bg is not shifted, the fitted curve looks fine!!
851 //-- 2012 Aug 8
852 //----> fail for 10h, not used, restore back!
853 //const Double_t lbg = TMath::Log(bg);
854
855 //redefine p2(<0) -> ln P2
856 //corresponds to alephFit.C fAlephOPt = 11
857 //const Double_t bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
858
6951a056 859 //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
860
861 return (p1-aa-bb)*p0/aa;
862}
863
ca66da8f 864Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
6951a056 865{
866 //
867 //f(x)-> f(y) with y=log10(x)
868 //
869 const Double_t x2[]={TMath::Power(10, xx[0])};
870 return func(x2, par);
871}
872