]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDdEdxUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Bug fix (Xianguo)
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxUtils.cxx
... / ...
CommitLineData
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// class to calculate TRD dEdx
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
30#include "TF1.h"
31#include "TFile.h"
32#include "TH1D.h"
33#include "TH2D.h"
34#include "THnSparse.h"
35#include "TMath.h"
36#include "TMatrixD.h"
37#include "TMinuit.h"
38#include "TObjArray.h"
39#include "TRandom3.h"
40#include "TStopwatch.h"
41#include "TVectorD.h"
42
43#include "TTreeStream.h"
44
45#include "AliCDBId.h"
46#include "AliCDBMetaData.h"
47#include "AliCDBStorage.h"
48#include "AliESDEvent.h"
49#include "AliESDfriendTrack.h"
50#include "AliESDtrack.h"
51#include "AliTRDtrackV1.h"
52
53#include "AliTRDdEdxUtils.h"
54
55#define EPSILON 1e-12
56
57THnSparseD * AliTRDdEdxUtils::fgHistGain=0x0;
58THnSparseD * AliTRDdEdxUtils::fgHistT0=0x0;
59THnSparseD * AliTRDdEdxUtils::fgHistVd=0x0;
60TObjArray * AliTRDdEdxUtils::fgHistPHQ=new TObjArray(8);
61
62TString AliTRDdEdxUtils::fgCalibFile;
63TObjArray * AliTRDdEdxUtils::fgObjGain = 0x0;
64TObjArray * AliTRDdEdxUtils::fgObjT0 = 0x0;
65TObjArray * AliTRDdEdxUtils::fgObjVd = 0x0;
66TObjArray * AliTRDdEdxUtils::fgObjPHQ = 0x0;
67
68Int_t AliTRDdEdxUtils::fgNchamber = -999;
69Double_t AliTRDdEdxUtils::fgChamberQ[6];
70Double_t AliTRDdEdxUtils::fgChamberTmean[6];
71
72Double_t AliTRDdEdxUtils::fgTrackTmean = -999;
73
74//===================================================================================
75// Math and Histogram
76//===================================================================================
77void AliTRDdEdxUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
78{
79 //
80 //counts of hh is sorted
81 //
82
83 for(Int_t ii=0; ii<ncut; ii++){
84 cuts[ii] = -999;
85 }
86
87 Int_t nsel = 0;
88 const Int_t nbin = hh->GetNbinsX();
89 Double_t datas[nbin];
90 for(Int_t ii=1; ii<=nbin; ii++){
91 const Double_t res = hh->GetBinContent(ii);
92 if(res<thres){
93 continue;
94 }
95
96 datas[nsel] = res;
97 nsel++;
98 }
99 if(!nsel)
100 return;
101
102 Int_t id[nsel];
103 TMath::Sort(nsel, datas, id, kFALSE);
104
105 for(Int_t ii=0; ii<ncut; ii++){
106 const Double_t icdf = cdfs[ii];
107 if(icdf<0 || icdf>1){
108 printf("AliTRDdEdxUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
109 }
110 cuts[ii] = datas[id[Int_t(icdf*nsel)]];
111 }
112}
113
114Double_t AliTRDdEdxUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
115{
116 //
117 //calculate mean (with error) and rms from sum, w2s, nn
118 //if nn=0, mean, error, and rms all = 0
119 //
120
121 Double_t tmean = 0, trms = 0, terr = 0;
122
123 if(nn>EPSILON){
124 tmean = sum/nn;
125
126 const Double_t arg = w2s/nn-tmean*tmean;
127 if(TMath::Abs(arg)<EPSILON){
128 trms = 0;
129 }
130 else{
131 if( arg <0 ){
132 printf("AliTRDdEdxUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
133 }
134
135 trms = TMath::Sqrt(arg);
136 }
137
138 terr = trms/TMath::Sqrt(nn);
139 }
140
141 if(grms){
142 (*grms) = trms;
143 }
144
145 if(gerr){
146 (*gerr) = terr;
147 }
148
149 return tmean;
150}
151
152Double_t AliTRDdEdxUtils::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)
153{
154 //
155 //calculate truncated mean
156 //return <x*w>_{low-high according to x}
157 //
158
159 /*
160 //test->
161 for(Int_t ii=0; ii<nx; ii++){
162 printf("test %d/%d %f\n", ii, nx, xdata[ii]);
163 }
164 //<--test
165 */
166
167 Int_t index[nx];
168 TMath::Sort(nx, xdata, index, kFALSE);
169
170 Int_t nused = 0;
171 Double_t sum = 0;
172 Double_t w2s = 0;
173 const Int_t istart = Int_t (nx*lowfrac);
174 const Int_t istop = Int_t (nx*highfrac);
175
176 //=,< correct, because when low=0, high=1 it is correct
177 for(Int_t ii=istart; ii<istop; ii++){
178 Double_t weight = 1;
179 if(wws){
180 weight = wws[index[ii]];
181 }
182 const Double_t sx = xdata[index[ii]]*weight;
183
184 sum += sx;
185 w2s += sx*sx;
186
187 nused++;
188 //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
189
190 }
191
192 return GetMeanRMS(nused, sum, w2s, grms, gerr);
193}
194
195Double_t AliTRDdEdxUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
196{
197 //
198 //do truncation on histogram
199 //
200 //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
201
202 //with under-/over-flow
203 Double_t npreTrunc = 0;
204 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
205 const Double_t be = hh->GetBinError(itmp);
206 const Double_t bc = hh->GetBinContent(itmp);
207 if(be<EPSILON){
208 if(bc>EPSILON){
209 printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
210 }
211 continue;
212 }
213 npreTrunc += bc*bc/be/be;
214 }
215
216 const Double_t nstart = npreTrunc*lowfrac;
217 const Double_t nstop = npreTrunc*highfrac;
218
219 //with Double_t this should also handle normalized hist
220 Double_t ntot = 0;
221 Double_t nused = 0;
222 Double_t sum = 0;
223 Double_t w2s = 0;
224 for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
225 const Double_t be = hh->GetBinError(itmp);
226 const Double_t bc = hh->GetBinContent(itmp);
227 if(be<EPSILON){
228 if(bc>EPSILON){
229 printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
230 }
231 continue;
232 }
233 const Double_t weight = bc*bc/be/be;
234 ntot+=weight;
235 //<= correct, because when high=1, nstop has to be included
236 if(ntot>nstart && ntot<=nstop){
237
238 const Double_t val = hh->GetBinCenter(itmp);
239 sum += weight*val;
240 w2s += weight*val*val;
241
242 nused += weight;
243
244 //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
245 }
246 else if(ntot>nstop){
247 if(itmp>=hh->GetNbinsX()){
248 printf("AliTRDdEdxUtils::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);
249 nused = 0;
250 w2s = sum = 0;
251 }
252 break;
253 }
254 }
255
256 return GetMeanRMS(nused, sum, w2s, grms, gerr);
257}
258
259void AliTRDdEdxUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
260{
261 //
262 //fit slices of hh using truncation
263 //
264
265 const Int_t x0 = hh->GetXaxis()->GetFirst();
266 const Int_t x1 = hh->GetXaxis()->GetLast();
267 const Int_t y0 = hh->GetYaxis()->GetFirst();
268 const Int_t y1 = hh->GetYaxis()->GetLast();
269
270 const Int_t nx = hh->GetNbinsX();
271 const Int_t ny = hh->GetNbinsY();
272 const Double_t xmin = hh->GetXaxis()->GetXmin();
273 const Double_t xmax = hh->GetXaxis()->GetXmax();
274 const Double_t ymin = hh->GetYaxis()->GetXmin();
275 const Double_t ymax = hh->GetYaxis()->GetXmax();
276
277 hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
278 hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
279 hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
280 hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
281
282 for(Int_t ix=x0; ix<=x1; ix++){
283 //to speed up
284 const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
285 if(rawcount<EPSILON){
286 continue;
287 }
288
289 TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
290 Double_t ntot = 0;
291 for(Int_t iy=y0; iy<=y1; iy++){
292 const Double_t be = hh->GetBinError(ix,iy);
293 const Double_t bc = hh->GetBinContent(ix, iy);
294
295 if(be<EPSILON){
296 if(bc>EPSILON){
297 printf("AliTRDdEdxUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
298 }
299 continue;
300 }
301
302 htmp->SetBinContent(iy, bc);
303 htmp->SetBinError(iy, be);
304
305 ntot += (bc/be)*(bc/be);
306
307 //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
308 }
309
310 hnor->SetBinContent(ix, ntot);
311 hnor->SetBinError( ix, 0);
312
313 if(ntot<thres || htmp->GetRMS()<EPSILON){
314 delete htmp;
315 continue;
316 }
317
318 //test htmp->Draw();
319 Double_t trms = -999, terr = -999;
320 const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
321
322 hmpv->SetBinContent(ix, tmean);
323 hmpv->SetBinError( ix, terr);
324
325 hwid->SetBinContent(ix, trms);
326 hwid->SetBinError( ix, 0);
327
328 hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
329 hres->SetBinError( ix, 0);
330
331 delete htmp;
332 }
333
334 TH1 *hhs[]={hnor, hmpv, hwid, hres};
335 const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
336 const Int_t nh = sizeof(hhs)/sizeof(TH1*);
337 for(Int_t ii=0; ii<nh; ii++){
338 hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
339 hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
340 hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
341 hhs[ii]->SetTitle(hh->GetTitle());
342 }
343}
344
345//===================================================================================
346// TRD Analysis Fast Tool
347//===================================================================================
348
349Int_t AliTRDdEdxUtils::GetNtracklet(const AliESDEvent *esd)
350{
351 //
352 //number of trd tracklet in one esd event
353 //
354 const Int_t ntrk0 = esd->GetNumberOfTracks();
355 Int_t ntrdv1=0;
356 for(Int_t ii=0; ii<ntrk0; ii++){
357 ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
358 }
359 return ntrdv1;
360}
361
362AliTRDtrackV1 * AliTRDdEdxUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
363{
364 //
365 //Get TRD friend track
366 //
367
368 AliESDfriendTrack * friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
369 if(!friendtrk){
370 //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
371 return 0x0;
372 }
373
374 TObject *calibObject=0x0;
375 AliTRDtrackV1 * trdtrack=0x0;
376 for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
377 if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
378 break;
379 }
380
381 return trdtrack;
382}
383
384Bool_t AliTRDdEdxUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
385{
386 //
387 // to check if all tracklets are in the same stack, useful in cosmic
388 //
389
390 TVectorD secs(18), stks(5);
391
392 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
393 AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
394 if(!tracklet)
395 continue;
396
397 const Int_t det = tracklet->GetDetector();
398 const Int_t isector = AliTRDgeometry::GetSector(det);
399 const Int_t istack = AliTRDgeometry::GetStack(det);
400
401 secs[isector] = 1;
402 stks[istack] = 1;
403 }
404
405 if(secs.Sum()!=1 || stks.Sum()!=1){
406 return kFALSE;
407 }
408 else
409 return kTRUE;
410}
411
412Bool_t AliTRDdEdxUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
413{
414 //
415 //as function name
416 //
417 isec = istk = -999;
418 mom = -999;
419
420 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
421 AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
422 if(!tracklet)
423 continue;
424
425 const Int_t det = tracklet->GetDetector();
426 isec = AliTRDgeometry::GetSector(det);
427 istk = AliTRDgeometry::GetStack(det);
428
429 mom = tracklet->GetMomentum();
430
431 break;
432 }
433
434 if(isec<0)
435 return kFALSE;
436 else
437 return kTRUE;
438}
439
440//===================================================================================
441// Calibration
442//===================================================================================
443Double_t AliTRDdEdxUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
444{
445 //
446 //the scale used in calibration
447 //
448
449 if(tpcncls < CalibTPCnclsCut())
450 return -999;
451
452 if(tpcsig<EPSILON)
453 return -999;
454
455 return tpcsig/120;
456
457}
458
459Int_t AliTRDdEdxUtils::GetPHQIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge)
460{
461 //
462 //iterator for calib obj and hist
463 //
464 return kinvq*4 + (mag>0)*2 + (charge>0);
465}
466
467TObjArray * AliTRDdEdxUtils::GetObjPHQ()
468{
469 //
470 //return fgObjPHQ, initialized if null
471 //
472
473 if(!fgObjPHQ){
474 fgObjPHQ = new TObjArray(8);
475 }
476
477 return fgObjPHQ;
478}
479
480TObjArray * AliTRDdEdxUtils::GetObjPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
481{
482 //
483 //return calib obj
484 //
485 if(!fgObjPHQ){
486 printf("AliTRDdEdxUtils::GetObjPHQ(kinvq, mag, charge) error fgObjPHQ null!!\n"); exit(1);
487 }
488
489 return (TObjArray*) fgObjPHQ->At(GetPHQIterator(kinvq, mag, charge));
490}
491
492THnSparseD * AliTRDdEdxUtils::GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
493{
494 //
495 //return calib hist
496 //
497 return (THnSparseD*) fgHistPHQ->At(GetPHQIterator(kinvq, mag, charge));
498}
499
500TString AliTRDdEdxUtils::GetPHQName(const Bool_t kobj, const Int_t iter)
501{
502 //
503 //get name of calib obj/hist of PHQ
504 //
505 return Form("TRDCalib%sPHQ%d", kobj?"Obj":"Hist", iter);
506}
507
508void AliTRDdEdxUtils::DeleteCalibObj()
509{
510 //
511 //delete calib obj
512 //
513 delete fgObjGain;
514 delete fgObjT0;
515 delete fgObjVd;
516
517 fgObjGain = 0x0;
518 fgObjT0 = 0x0;
519 fgObjVd = 0x0;
520
521 if(fgObjPHQ){
522 fgObjPHQ->SetOwner();
523 delete fgObjPHQ;
524 fgObjPHQ = 0x0;
525 }
526}
527
528Bool_t AliTRDdEdxUtils::GenerateDefaultPHQOCDB(const TString path)
529{
530 //
531 //generate default OCDB object PHQ, do like
532 //AliTRDdEdxUtils::GenerateDefaultPHQOCDB("local://./")
533 //
534
535 TObjArray * arr8 = new TObjArray(8);
536 arr8->SetOwner();
537
538 for(Int_t ii=0; ii<8; ii++){
539 TObjArray * arr1 = new TObjArray(1);
540 arr1->SetOwner();
541 arr1->SetName(GetPHQName(1, ii));
542
543 const Int_t nbins = NTRDtimebin();
544 TVectorD * vec = new TVectorD(nbins);
545 for(Int_t jj=0; jj<nbins; jj++){
546 (*vec)[jj] = 1;
547 }
548 arr1->Add(vec);
549 arr8->Add(arr1);
550 }
551
552 AliCDBMetaData *metaData= new AliCDBMetaData();
553 metaData->SetObjectClassName("TObjArray");
554 metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
555
556 AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
557 AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
558 gStorage->Put(arr8, id1, metaData);
559
560 delete metaData;
561 delete arr8;
562
563 return kTRUE;
564}
565
566void AliTRDdEdxUtils::IniCalibObj()
567{
568 //
569 //set CalibObj from file, clone to static calib obj
570 //
571
572 DeleteCalibObj();
573
574 TFile *cfile=new TFile(fgCalibFile);
575 if(!cfile->IsOpen()){
576 printf("AliTRDdEdxUtils::IniCalibObj error fgCalibFile not open! %s\n", fgCalibFile.Data());exit(1);
577 }
578
579 printf("\nAliTRDdEdxUtils::IniCalibObj file: %s\n", fgCalibFile.Data());
580
581 //---
582 const TString objnames[] ={"TRDCalibObjGain", "TRDCalibObjT0", "TRDCalibObjVd"};
583 TObjArray ** gobjs[]={ &fgObjGain, &fgObjT0, &fgObjVd};
584
585 const Int_t nobj = sizeof(objnames)/sizeof(TString);
586 for(Int_t iobj=0; iobj<nobj; iobj++){
587 TObjArray *tmpo=0x0;
588 cfile->GetObject(objnames[iobj], tmpo);
589 if(!tmpo){
590 printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objnames[iobj].Data()); exit(1);
591 }
592
593 (*gobjs[iobj])=(TObjArray*)tmpo->Clone();
594 (*gobjs[iobj])->SetOwner();
595 }
596
597 fgObjPHQ = new TObjArray(8);
598 for(Int_t iter=0; iter<8; iter++){
599 const TString objn = GetPHQName(1, iter);
600 TObjArray *tmpo=0x0;
601 cfile->GetObject(objn, tmpo);
602 if(!tmpo){
603 printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objn.Data()); exit(1);
604 }
605
606 TObjArray *obji=(TObjArray*) tmpo->Clone();
607 obji->SetOwner();
608 fgObjPHQ->AddAt(obji, iter);
609 }
610
611 //---
612
613 cfile->Close();
614 delete cfile;
615}
616
617void AliTRDdEdxUtils::DeleteCalibHist()
618{
619 //
620 //delete calib hist
621 //
622 delete fgHistGain;
623 delete fgHistT0;
624 delete fgHistVd;
625
626 fgHistGain = 0x0;
627 fgHistT0 = 0x0;
628 fgHistVd = 0x0;
629
630 //fgHistPHQ owns the hists
631 fgHistPHQ->SetOwner();
632 fgHistPHQ->Clear();
633}
634
635void AliTRDdEdxUtils::IniCalibHist(TList *list, const Bool_t kPHQonly)
636{
637 //
638 //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
639 //
640
641 DeleteCalibHist();
642
643 Int_t nbin[2];
644 const Double_t xmin[2]={0, 0};
645 Double_t xmax[2];
646
647 nbin[0]=NTRDtimebin(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20;
648 for(Int_t iter=0; iter<8; iter++){
649 const TString hn = GetPHQName(0, iter);
650 THnSparseD *hi = new THnSparseD(hn, "", 2, nbin, xmin, xmax);
651 //fgHistPHQ owns the hists
652 fgHistPHQ->AddAt(hi, iter);
653 list->Add(hi);
654 }
655
656 if(kPHQonly)
657 return;
658
659 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20; fgHistGain = new THnSparseD("TRDCalibHistGain", "", 2, nbin, xmin, xmax);
660 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistT0 = new THnSparseD("TRDCalibHistT0", "", 2, nbin, xmin, xmax);
661 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistVd = new THnSparseD("TRDCalibHistVd", "", 2, nbin, xmin, xmax);
662
663 list->Add(fgHistGain);
664 list->Add(fgHistT0);
665 list->Add(fgHistVd);
666}
667
668Bool_t AliTRDdEdxUtils::ReadCalibHist(const TString filename, const TString listname)
669{
670 //
671 //used in AliTRDPreprocessorOffline
672 //read in calib hist from file, only for PHQ
673 //
674 DeleteCalibHist();
675
676 //maybe already open by others... don't close
677 TFile fcalib(filename);
678
679 TObjArray * array = (TObjArray*)fcalib.Get(listname);
680
681 for(Int_t iter=0; iter<8; iter++){
682 const TString hn = GetPHQName(0, iter);
683 THnSparseD * tmph=0x0;
684 if(array){
685 tmph = (THnSparseD *) array->FindObject(hn);
686 }
687 else{
688 tmph = (THnSparseD *) fcalib.Get(hn);
689 }
690 if(!tmph){
691 printf("AliTRDdEdxUtils::ReadCalibHist warning calib hist not found! %s %s\n", filename.Data(), listname.Data());
692 fcalib.ls();
693 if(array){
694 array->ls();
695 }
696 return kFALSE;
697 }
698 THnSparseD *hi = (THnSparseD*)tmph->Clone();
699 fgHistPHQ->AddAt(hi, iter);
700 }
701
702 return kTRUE;
703}
704
705void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnSparseD * hcalib, const Double_t scale)
706{
707 //
708 //fill calibration hist
709 //
710 if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
711
712 for(Int_t ii=0; ii<ncls; ii++){
713 const Double_t dq = (*arrayQ)[ii];
714 const Double_t xx = (*arrayX)[ii];
715
716 const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
717 const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
718 const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
719
720 if(xx>=xmax || xx<xmin){
721 printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
722 }
723
724 const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
725 hcalib->Fill(var);
726 }
727}
728
729void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
730{
731 //
732 //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
733 //
734
735 THnSparseD * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge);
736
737 TVectorD arrayQ(200), arrayX(200);
738 const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
739 FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
740
741 static Int_t kprint = 100;
742 if(kprint<0){
743 printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n");
744 printf("\nkinvq= %d;\n", kinvq);
745 for(Int_t iq=0; iq<ncls; iq++){
746 printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
747 }
748 printf("\n");
749 }
750 kprint++;
751}
752
753Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
754{
755 //
756 //apply calibration on arrayQ
757 //
758 if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);}
759
760 TVectorD tmpq(arrayQ->GetNrows());
761 TVectorD tmpx(arrayX->GetNrows());
762 Int_t ncls = 0;
763
764 const TVectorD * gain = (TVectorD*) cobj->At(0);
765 for(Int_t ii=0; ii<nc0; ii++){
766 const Double_t dq = (*arrayQ)[ii];
767 const Int_t xx = (Int_t)(*arrayX)[ii];
768 const Double_t gg = (*gain)[xx];
769
770 if(gg<EPSILON){
771 continue;
772 }
773
774 tmpq[ncls] = dq*gg;
775 tmpx[ncls] = xx;
776 ncls++;
777 }
778
779 (*arrayQ)=tmpq;
780 (*arrayX)=tmpx;
781
782 return ncls;
783}
784
785void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
786{
787 //
788 //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
789 //
790 const Int_t ndet = 540;
791 TObjArray *obj=new TObjArray(ndet);
792 obj->SetOwner();
793 for(Int_t ii=0; ii<ndet; ii++){
794 obj->Add(new TVectorD(AliTRDseedV1::kNtb));
795 }
796
797 //ibin = binlowedge of bin(ibin+1) = the number fills this bin
798 for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
799 const Double_t stat = hnor->GetBinContent(ibin+1);
800 if(stat<EPSILON){
801 continue;
802 }
803
804 const Int_t idet = ToDetector(ibin);
805 const Int_t itb = ToTimeBin(ibin);
806 TVectorD *vec=(TVectorD *)obj->At(idet);
807 (*vec)[itb] = stat;
808 }
809
810 hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
811 for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
812 const Int_t idet = ToDetector(ibin);
813 const TVectorD *vec=(TVectorD *)obj->At(idet);
814
815 Int_t nonzero = 0;
816 for(Int_t ii=0; ii<vec->GetNrows(); ii++){
817 if((*vec)[ii]>EPSILON){
818 nonzero++;
819 }
820 }
821
822 Double_t mean = 0;
823 const Double_t lowfrac = 0.6;
824 //if there are too many 0's, reject this chamber by settig mean=rms=0
825 if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
826 //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
827 mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
828 }
829
830 hmean->SetBinContent(ibin+1, mean);
831 }
832
833 delete obj;
834}
835
836void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run)
837{
838 //
839 //produce calibration objects
840 //
841
842 TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd ");
843 for(Int_t iter=0; iter<8; iter++){
844 objnames+= GetPHQName(0, iter)+" ";
845 }
846
847 TList * lout = new TList;
848 lout->SetOwner();
849
850 TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
851
852 const Int_t nh=lin->GetEntries();
853 for(Int_t ii=0; ii<nh; ii++){
854 const THnSparseD *hh=(THnSparseD*)lin->At(ii);
855 const TString hname = hh->GetName();
856 if(!objnames.Contains(hname))
857 continue;
858
859 TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream);
860 lout->Add(cobj0);
861 }
862
863 //lout->ls();
864
865 //=============================================================
866 //=============================================================
867
868 TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
869 fout->cd();
870 const Int_t nout=lout->GetEntries();
871 for(Int_t ii=0; ii<nout; ii++){
872 const TString oname = lout->At(ii)->GetName();
873 if(oname.Contains("Obj")){
874 TObjArray * cobj = (TObjArray*) lout->At(ii);
875 cobj->Write(oname, TObjArray::kSingleKey);
876 }
877 }
878 fout->Save();
879 fout->Close();
880 delete fout;
881
882 fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
883 fout->cd();
884 lin->Write();
885 lout->Write();
886 fout->Save();
887 fout->Close();
888 delete fout;
889
890 delete calibStream;
891
892 /*
893 http://root.cern.ch/root/html/TH1.html
894 When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
895
896 h->SetDirectory(0); for the current histogram h
897 TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
898
899 When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted.
900 */
901 delete lout;
902}
903
904TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnSparseD *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
905{
906 //
907 //produce calibration objects
908 //
909
910 const TString hname = hh->GetName();
911 const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
912
913 //----------------------------------------
914 // Define nbin, tag, cobj0
915 //----------------------------------------
916 Int_t nbin =-999;
917 if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){
918 nbin = NTRDchamber();
919 }
920 else if(hname.Contains("PHQ")){
921 nbin = NTRDtimebin();
922 }
923 else{
924 printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1);
925 }
926
927 TString tag(hname);
928 tag.ReplaceAll("Hist","Obj");
929
930 TObjArray * cobj0 = new TObjArray(1);
931 cobj0->SetOwner();
932 cobj0->SetName(tag);
933 cobj0->Add(new TVectorD(nbin));
934
935 //----------------------------------------
936 // Define lowFrac, highFrac
937 //----------------------------------------
938 Double_t lowFrac = -999, highFrac = -999;
939 if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){
940 lowFrac = 0.01; highFrac = Q0Frac();
941 }
942 else if(hname.Contains("PHQ") && kinvq){
943 lowFrac = Q1Frac(); highFrac = 0.99;
944 }
945 else{
946 lowFrac = 0.01;
947 highFrac = 0.99;
948 }
949
950 //----------------------------------------
951 // Get analysis result
952 //----------------------------------------
953 TH1::AddDirectory(kFALSE);//important!
954 TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
955 TH2D *hpj = hh->Projection(1,0);
956 FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
957 if(hname.Contains("PHQ")){
958 GetPHCountMeanRMS(hnor, htrdphmean);
959 if(lout){
960 lout->Add(htrdphmean);
961 }
962 }
963 delete hpj;
964
965 if(lout){
966 lout->Add(hnor);
967 lout->Add(hmpv);
968 lout->Add(hwid);
969 lout->Add(hres);
970 }
971
972 //----------------------------------------
973 // Define Counter
974 //----------------------------------------
975 TVectorD *countDet=0x0;
976 TObjArray *countSSL=0x0;
977
978 if(hname.Contains("PHQ") && !kinvq){
979 countDet=new TVectorD(540);
980 countSSL=new TObjArray(90);//SectorStackLayer
981 countSSL->SetOwner();
982 for(Int_t ii=0; ii<90; ii++){
983 countSSL->Add(new TVectorD(6));
984 }
985 }
986
987 //----------------------------------------
988 // Fill cobj0
989 //----------------------------------------
990
991 //ibin = binlowedge of bin(ibin+1) = the number fills this bin
992 for(Int_t ibin=0; ibin<nbin; ibin++){
993 Double_t gnor = hnor->GetBinContent(ibin+1);
994 Double_t gmpv = hmpv->GetBinContent(ibin+1);
995 Double_t gwid = hwid->GetBinContent(ibin+1);
996 Double_t gres = hres->GetBinContent(ibin+1);
997
998 //--- set additional cut by kpass
999 Bool_t kpass = kTRUE;
1000 Double_t gtrdphmean = -999;
1001 if(htrdphmean){
1002 gtrdphmean = htrdphmean->GetBinContent(ibin+1);
1003 //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
1004 if(gtrdphmean<EPSILON){
1005 kpass = kFALSE;
1006 }
1007 if(gnor<TimeBinCountCut()*gtrdphmean){
1008 kpass = kFALSE;
1009 }
1010 }
1011
1012 //--- set calibration constant p0
1013 Double_t p0= 0;
1014
1015 //reason for gmpv=0:
1016 //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
1017 //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
1018
1019 if(gmpv>EPSILON && kpass){
1020 if(tag.Contains("T0")){
1021 p0 = gmpv;
1022 }
1023 else{
1024 p0 = 1/gmpv;
1025 }
1026 //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
1027 }
1028
1029 (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
1030
1031 //--- save optional record
1032 if(p0>EPSILON && countDet && countSSL){
1033 const Int_t idet = ToDetector(ibin);
1034 (*countDet)[idet]=1;
1035
1036 const Int_t isector = ToSector(ibin);
1037 const Int_t istack = ToStack(ibin);
1038 const Int_t ilayer = ToLayer(ibin);
1039 TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
1040 (*vecsectorstack)[ilayer]=1;
1041 }
1042
1043 if(calibStream){
1044 (*calibStream)<<tag<<
1045 "run="<<run<<
1046 "p0="<<p0<<
1047
1048 "nor="<<gnor<<
1049 "mpv="<<gmpv<<
1050 "wid="<<gwid<<
1051 "res="<<gres<<
1052 "gtrdphmean="<<gtrdphmean<<
1053
1054 "ibin="<<ibin<<
1055 "\n";
1056 }
1057 }
1058
1059 //----------------------------------------
1060 // Status Report
1061 //----------------------------------------
1062 if(countDet && countSSL){
1063 TVectorD count2Dstack(90);
1064 for(Int_t ii=0; ii<90; ii++){
1065 TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
1066 const Int_t nlayer = (Int_t)vecsectorstack->Sum();
1067 if(nlayer==6){
1068 count2Dstack[ii]=1;
1069 }
1070 }
1071
1072 printf("\nAliTRDdEdxUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
1073 }
1074
1075 //----------------------------------------
1076 // Clean Up
1077 //----------------------------------------
1078
1079 TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
1080 const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
1081 for(Int_t ihh=0; ihh<nhh; ihh++){
1082 if(!lout){
1083 delete (*hhs[ihh]);
1084 }
1085 }
1086
1087 delete countDet;
1088 delete countSSL;
1089
1090 //----------------------------------------
1091
1092 return cobj0;
1093}
1094
1095//===================================================================================
1096// dEdx calculation
1097//===================================================================================
1098Double_t AliTRDdEdxUtils::GetSignal(const Int_t nch, const Int_t ncls, const Double_t qq)
1099{
1100 //
1101 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1102 //
1103 if(ncls>=1e3){
1104 printf("AliTRDdEdxUtils::GetSignal error ncls>1e3! %d\n", ncls); exit(1);
1105 }
1106
1107 return nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq);
1108}
1109
1110Int_t AliTRDdEdxUtils::GetNch(const Double_t signal)
1111{
1112 //
1113 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1114 //
1115 return Int_t(signal/1e6);
1116
1117}
1118
1119Int_t AliTRDdEdxUtils::GetNcls(const Double_t signal)
1120{
1121 //
1122 //signal = Nch*1e6 + Ncls*1e3 + (Q>=1e3? 999 : Q)
1123 //
1124
1125 return Int_t ( (signal-GetNch(signal)*1e6)/1e3 );
1126}
1127
1128Double_t AliTRDdEdxUtils::GetQ(const Double_t signal)
1129{
1130 //
1131 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1132 //
1133
1134 return signal-GetNch(signal)*1e6 - GetNcls(signal)*1e3;
1135}
1136
1137Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
1138{
1139 //
1140 //template for cookdedx
1141 //
1142 if(cobj){
1143 if(arrayQ && arrayX){
1144 ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
1145 }
1146 else{
1147 printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
1148 }
1149 }
1150
1151 Double_t lowFrac =-999, highFrac = -999;
1152 if(kinvq){
1153 lowFrac = Q1Frac(); highFrac = 0.99;
1154 }
1155 else{
1156 lowFrac = 0.01; highFrac = Q0Frac();
1157 }
1158
1159 Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
1160 if(kinvq){
1161 if(meanQ>EPSILON){
1162 meanQ = 1/meanQ;
1163 }
1164 }
1165
1166 return meanQ;
1167}
1168
1169Double_t AliTRDdEdxUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
1170{
1171 //
1172 //combine tpc and trd dedx
1173 //
1174
1175 for(Int_t iq=0; iq<tpcncls; iq++){
1176 (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
1177 if(tpcarrayX && trdarrayX && coarrayX){
1178 (*coarrayX)[iq]=(*tpcarrayX)[iq];
1179 }
1180 }
1181 for(Int_t iq=0; iq<trdncls; iq++){
1182 (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
1183 if(tpcarrayX && trdarrayX && coarrayX){
1184 (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
1185 }
1186 }
1187
1188 concls=trdncls+tpcncls;
1189
1190 const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
1191
1192 return coQ;
1193}
1194
1195
1196//===================================================================================
1197// dEdx Getter and Setter
1198//===================================================================================
1199Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
1200{
1201 //
1202 //return angular normalization factor
1203 //
1204
1205 return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
1206}
1207
1208Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
1209{
1210 //
1211 //get cluster charge
1212 //
1213 Double_t dq = 0;
1214 const AliTRDcluster *cl = 0x0;
1215
1216 cl = seed->GetClusters(itb); if(cl) dq+=cl->GetRawQ();
1217 cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl) dq+=cl->GetRawQ();
1218
1219 dq /= GetAngularCorrection(seed);
1220
1221 dq /= 45.;
1222
1223 if(kinvq){
1224 if(dq>EPSILON){
1225 dq = 1/dq;
1226 }
1227 }
1228
1229 return dq;
1230}
1231
1232Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
1233{
1234 //
1235 //return nclustter
1236 //(if kinvq, return 1/q array), size of array must be larger than 31*6
1237 //
1238 if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1239 printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
1240 }
1241 if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1242 printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
1243 }
1244
1245 const Int_t mintb = 0;
1246 const Int_t maxtb = AliTRDseedV1::kNtb-1;
1247 if(timeBin0<mintb) timeBin0=mintb;
1248 if(timeBin1>maxtb) timeBin1=maxtb;
1249 if(tstep<=0) tstep=1;
1250
1251 //============
1252 Int_t tbN=0;
1253 Double_t tbQ[200];
1254 Int_t tbBin[200];
1255
1256 for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1257 const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1258 if(!seed)
1259 continue;
1260
1261 const Int_t det = seed->GetDetector();
1262
1263 for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
1264 const Double_t dq = GetClusterQ(kinvq, seed, itb);
1265 if(dq<EPSILON)
1266 continue;
1267
1268 const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
1269
1270 tbQ[tbN]=dq;
1271 tbBin[tbN]=gtb;
1272 tbN++;
1273 }
1274 }
1275
1276 Int_t ncls = 0;
1277 for(Int_t iq=0; iq<tbN; iq++){
1278 if(tbQ[iq]<EPSILON)
1279 continue;
1280
1281 (*arrayQ)[ncls] = tbQ[iq];
1282 (*arrayX)[ncls] = tbBin[iq];
1283
1284 ncls++;
1285 }
1286
1287 static Int_t kprint = 100;
1288 if(kprint<0){
1289 printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n");
1290 for(Int_t iq=0; iq<ncls; iq++){
1291 const Int_t ichamber = ToLayer((*arrayX)[iq]);
1292 const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1293 if(!seed){
1294 printf("error seed null!!\n"); exit(1);
1295 }
1296 const Double_t rawq = (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
1297 printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
1298 }
1299 printf("\n");
1300 }
1301 kprint++;
1302
1303 return ncls;
1304}
1305
1306Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
1307{
1308 //
1309 //arrayX det*Ntb+itb -> itb
1310 //
1311
1312 TVectorD countChamber(6);
1313 for(Int_t ii = 0; ii<ncls; ii++){
1314 const Int_t xx = (Int_t)(*arrayX)[ii];
1315 const Int_t idet = ToDetector(xx);
1316
1317 const Double_t ich = AliTRDgeometry::GetLayer(idet);
1318 const Double_t itb = ToTimeBin(xx);
1319 (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
1320
1321 countChamber[ich] = 1;
1322 }
1323
1324 const Double_t nch = countChamber.Sum();
1325 return (Int_t) nch;
1326}
1327
1328void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnSparseD * hgain, THnSparseD * ht0, THnSparseD * hvd)
1329{
1330 //
1331 //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution
1332 //
1333
1334 static Int_t kprint = 100;
1335
1336 fgNchamber = 0;
1337 for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1338 //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities
1339 fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0;
1340
1341 const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber);
1342 if (!seed)
1343 continue;
1344
1345 const Int_t idet = seed->GetDetector();
1346
1347 //-------------------------------------------------------------------------
1348
1349 Double_t qsum = 0, qtsum = 0, w2sum = 0;
1350 for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
1351 const Double_t dq = GetClusterQ(0, seed, itb);
1352 if(dq<EPSILON)
1353 continue;
1354
1355 qsum += dq;
1356 qtsum += dq*itb;
1357 w2sum += dq*itb*itb;
1358 }
1359 if(qsum<EPSILON)
1360 continue;
1361
1362 //-------------------------------------------------------------------------
1363
1364 Double_t tbm, tbr = 0;
1365 tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr);
1366
1367 qsum /= 1.25e3/45.;
1368
1369 if(hgain){
1370 const Double_t var[]={idet, qsum};
1371 hgain->Fill(var);
1372 }
1373 if(ht0){
1374 const Double_t var[]={idet, tbm};
1375 ht0->Fill(var);
1376 }
1377 if(hvd){
1378 const Double_t var[]={idet, tbr};
1379 hvd->Fill(var);
1380 }
1381
1382 Double_t gain = 1, t0 = 0, vd = 1;
1383 if(kcalib){
1384 if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);}
1385 if(! fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0 array null!!\n"); exit(1);}
1386 if(! fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd array null!!\n"); exit(1);}
1387
1388 const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet];
1389 const TVectorD * t0vec = (TVectorD*) fgObjT0->At(0); t0 = (* t0vec)[idet];
1390 const TVectorD * vdvec = (TVectorD*) fgObjVd->At(0); vd = (* vdvec)[idet];
1391 }
1392 if(kprint<0){
1393 printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1394 printf("idet = %d;\n", idet);
1395 printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd);
1396 printf("\n");
1397 }
1398
1399 qsum *= gain;
1400 tbm = (tbm-t0)*vd;
1401
1402 if(qsum<EPSILON)
1403 continue;
1404
1405 //-------------------------------------------------------------------------
1406
1407 //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0)
1408 fgChamberQ[ichamber] = qsum;
1409 fgChamberTmean[ichamber] = tbm;
1410 fgNchamber++;
1411 }
1412
1413 if(kprint<0){
1414 printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n");
1415
1416 printf("\nfgNchamber = %d\n", fgNchamber);
1417 for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1418 printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]);
1419 }
1420 }
1421
1422 fgTrackTmean = -999;
1423 if(fgNchamber){
1424 fgTrackTmean = 0;
1425 for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1426 fgTrackTmean += fgChamberTmean[ich];
1427 }
1428 fgTrackTmean /= fgNchamber;
1429 }
1430
1431 if(kprint<0){
1432 printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1433 printf("GetTrackTmean() %15f\n", GetTrackTmean());
1434 printf("\n");
1435 }
1436 kprint++;
1437
1438 return;
1439}
1440
1441
1442//===================================================================================
1443// dEdx Parameterization
1444//===================================================================================
1445
1446Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg)
1447{
1448 //
1449 //truncated Mean Q_{xx} in TRD
1450 //
1451
1452 Double_t par[8];
1453 //03132012161150
1454 //opt: ppQ0
1455par[0]= 2.397001e-01;
1456par[1]= 1.334697e+00;
1457par[2]= 6.967470e+00;
1458par[3]= 9.055289e-02;
1459par[4]= 9.388760e+00;
1460par[5]= 9.452742e-04;
1461par[6]= -1.866091e+00;
1462par[7]= 1.403545e+00;
1463
1464 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000
1465 return 0.428543*MeandEdxTR(&bg, par);
1466}
1467
1468Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg)
1469{
1470 //
1471 //truncated Mean Q_{xx} in TRD
1472 //
1473
1474 Double_t par[8];
1475
1476 //03132012161259
1477 //opt: PbPbQ0
1478par[0]= 1.844912e-01;
1479par[1]= 2.509702e+00;
1480par[2]= 6.744031e+00;
1481par[3]= 7.355123e-02;
1482par[4]= 1.166023e+01;
1483par[5]= 1.736186e-04;
1484par[6]= -1.716063e+00;
1485par[7]= 1.611366e+00;
1486
1487 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
1488 return 0.460994*MeandEdxTR(&bg, par);
1489}
1490
1491Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg)
1492{
1493 //
1494 //truncated Mean 1/(1/Q)_{xx} in TRD
1495 //
1496
1497 Double_t par[8];
1498
1499