Isolation : Fiducial cut applied on detector origin of trigger and not calorimeter...
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxUtils.cxx
CommitLineData
fb3d369e 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(const Bool_t kinvq, const Double_t mag, const Int_t charge)
468{
469 //
470 //return calib obj
471 //
472 if(!fgObjPHQ){
473 printf("AliTRDdEdxUtils::GetObjPHQ error fgObjPHQ null!!\n"); exit(1);
474 }
475
476 return (TObjArray*) fgObjPHQ->At(GetPHQIterator(kinvq, mag, charge));
477}
478
479THnSparseD * AliTRDdEdxUtils::GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
480{
481 //
482 //return calib hist
483 //
484 return (THnSparseD*) fgHistPHQ->At(GetPHQIterator(kinvq, mag, charge));
485}
486
487TString AliTRDdEdxUtils::GetPHQName(const Bool_t kobj, const Int_t iter)
488{
489 //
490 //get name of calib obj/hist of PHQ
491 //
492 return Form("TRDCalib%sPHQ%d", kobj?"Obj":"Hist", iter);
493}
494
495void AliTRDdEdxUtils::DeleteCalibObj()
496{
497 //
498 //delete calib obj
499 //
500 delete fgObjGain;
501 delete fgObjT0;
502 delete fgObjVd;
503
504 fgObjGain = 0x0;
505 fgObjT0 = 0x0;
506 fgObjVd = 0x0;
507
508 if(fgObjPHQ){
509 fgObjPHQ->SetOwner();
510 delete fgObjPHQ;
511 fgObjPHQ = 0x0;
512 }
513}
514
515Bool_t AliTRDdEdxUtils::GenerateDefaultPHQOCDB(const TString path)
516{
517 //
518 //generate default OCDB object PHQ, do like
519 //AliTRDdEdxUtils::GenerateDefaultPHQOCDB("local://./")
520 //
521
522 TObjArray * arr8 = new TObjArray(8);
523 arr8->SetOwner();
524
525 for(Int_t ii=0; ii<8; ii++){
526 TObjArray * arr1 = new TObjArray(1);
527 arr1->SetOwner();
528 arr1->SetName(GetPHQName(1, ii));
529
530 const Int_t nbins = NTRDtimebin();
531 TVectorD * vec = new TVectorD(nbins);
532 for(Int_t jj=0; jj<nbins; jj++){
533 (*vec)[jj] = 1;
534 }
535 arr1->Add(vec);
536 arr8->Add(arr1);
537 }
538
539 AliCDBMetaData *metaData= new AliCDBMetaData();
540 metaData->SetObjectClassName("TObjArray");
541 metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
542
543 AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
544 AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
545 gStorage->Put(arr8, id1, metaData);
546
547 delete metaData;
548 delete arr8;
549
550 return kTRUE;
551}
552
553void AliTRDdEdxUtils::IniCalibObj()
554{
555 //
556 //set CalibObj from file, clone to static calib obj
557 //
558
559 DeleteCalibObj();
560
561 TFile *cfile=new TFile(fgCalibFile);
562 if(!cfile->IsOpen()){
563 printf("AliTRDdEdxUtils::IniCalibObj error fgCalibFile not open! %s\n", fgCalibFile.Data());exit(1);
564 }
565
566 printf("\nAliTRDdEdxUtils::IniCalibObj file: %s\n", fgCalibFile.Data());
567
568 //---
569 const TString objnames[] ={"TRDCalibObjGain", "TRDCalibObjT0", "TRDCalibObjVd"};
570 TObjArray ** gobjs[]={ &fgObjGain, &fgObjT0, &fgObjVd};
571
572 const Int_t nobj = sizeof(objnames)/sizeof(TString);
573 for(Int_t iobj=0; iobj<nobj; iobj++){
574 TObjArray *tmpo=0x0;
575 cfile->GetObject(objnames[iobj], tmpo);
576 if(!tmpo){
577 printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objnames[iobj].Data()); exit(1);
578 }
579
580 (*gobjs[iobj])=(TObjArray*)tmpo->Clone();
581 (*gobjs[iobj])->SetOwner();
582 }
583
584 fgObjPHQ = new TObjArray(8);
585 for(Int_t iter=0; iter<8; iter++){
586 const TString objn = GetPHQName(1, iter);
587 TObjArray *tmpo=0x0;
588 cfile->GetObject(objn, tmpo);
589 if(!tmpo){
590 printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objn.Data()); exit(1);
591 }
592
593 TObjArray *obji=(TObjArray*) tmpo->Clone();
594 obji->SetOwner();
595 fgObjPHQ->AddAt(obji, iter);
596 }
597
598 //---
599
600 cfile->Close();
601 delete cfile;
602}
603
604void AliTRDdEdxUtils::DeleteCalibHist()
605{
606 //
607 //delete calib hist
608 //
609 delete fgHistGain;
610 delete fgHistT0;
611 delete fgHistVd;
612
613 fgHistGain = 0x0;
614 fgHistT0 = 0x0;
615 fgHistVd = 0x0;
616
617 //fgHistPHQ owns the hists
618 fgHistPHQ->SetOwner();
619 fgHistPHQ->Clear();
620}
621
622void AliTRDdEdxUtils::IniCalibHist(TList *list, const Bool_t kPHQonly)
623{
624 //
625 //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
626 //
627
628 DeleteCalibHist();
629
630 Int_t nbin[2];
631 const Double_t xmin[2]={0, 0};
632 Double_t xmax[2];
633
634 nbin[0]=NTRDtimebin(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20;
635 for(Int_t iter=0; iter<8; iter++){
636 const TString hn = GetPHQName(0, iter);
637 THnSparseD *hi = new THnSparseD(hn, "", 2, nbin, xmin, xmax);
638 //fgHistPHQ owns the hists
639 fgHistPHQ->AddAt(hi, iter);
640 list->Add(hi);
641 }
642
643 if(kPHQonly)
644 return;
645
646 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20; fgHistGain = new THnSparseD("TRDCalibHistGain", "", 2, nbin, xmin, xmax);
647 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistT0 = new THnSparseD("TRDCalibHistT0", "", 2, nbin, xmin, xmax);
648 nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistVd = new THnSparseD("TRDCalibHistVd", "", 2, nbin, xmin, xmax);
649
650 list->Add(fgHistGain);
651 list->Add(fgHistT0);
652 list->Add(fgHistVd);
653}
654
655Bool_t AliTRDdEdxUtils::ReadCalibHist(const TString filename, const TString listname)
656{
657 //
658 //used in AliTRDPreprocessorOffline
659 //read in calib hist from file, only for PHQ
660 //
661 DeleteCalibHist();
662
663 //maybe already open by others... don't close
664 TFile fcalib(filename);
665
666 TObjArray * array = (TObjArray*)fcalib.Get(listname);
667
668 for(Int_t iter=0; iter<8; iter++){
669 const TString hn = GetPHQName(0, iter);
670 THnSparseD * tmph=0x0;
671 if(array){
672 tmph = (THnSparseD *) array->FindObject(hn);
673 }
674 else{
675 tmph = (THnSparseD *) fcalib.Get(hn);
676 }
677 if(!tmph){
678 printf("AliTRDdEdxUtils::ReadCalibHist warning calib hist not found! %s %s\n", filename.Data(), listname.Data());
679 fcalib.ls();
879b7730 680 if(array){
681 array->ls();
682 }
fb3d369e 683 return kFALSE;
684 }
685 THnSparseD *hi = (THnSparseD*)tmph->Clone();
686 fgHistPHQ->AddAt(hi, iter);
687 }
688
689 return kTRUE;
690}
691
692void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnSparseD * hcalib, const Double_t scale)
693{
694 //
695 //fill calibration hist
696 //
697 if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
698
699 for(Int_t ii=0; ii<ncls; ii++){
700 const Double_t dq = (*arrayQ)[ii];
701 const Double_t xx = (*arrayX)[ii];
702
703 const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
704 const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
705 const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
706
707 if(xx>=xmax || xx<xmin){
708 printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
709 }
710
711 const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
712 hcalib->Fill(var);
713 }
714}
715
716void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
717{
718 //
719 //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
720 //
721
722 THnSparseD * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge);
723
724 TVectorD arrayQ(200), arrayX(200);
725 const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
726 FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
727
728 static Int_t kprint = 100;
729 if(kprint<0){
730 printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n");
731 printf("\nkinvq= %d;\n", kinvq);
732 for(Int_t iq=0; iq<ncls; iq++){
733 printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
734 }
735 printf("\n");
736 }
737 kprint++;
738}
739
740Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
741{
742 //
743 //apply calibration on arrayQ
744 //
745 if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);}
746
747 TVectorD tmpq(arrayQ->GetNrows());
748 TVectorD tmpx(arrayX->GetNrows());
749 Int_t ncls = 0;
750
751 const TVectorD * gain = (TVectorD*) cobj->At(0);
752 for(Int_t ii=0; ii<nc0; ii++){
753 const Double_t dq = (*arrayQ)[ii];
754 const Int_t xx = (Int_t)(*arrayX)[ii];
755 const Double_t gg = (*gain)[xx];
756
757 if(gg<EPSILON){
758 continue;
759 }
760
761 tmpq[ncls] = dq*gg;
762 tmpx[ncls] = xx;
763 ncls++;
764 }
765
766 (*arrayQ)=tmpq;
767 (*arrayX)=tmpx;
768
769 return ncls;
770}
771
772void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
773{
774 //
775 //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
776 //
777 const Int_t ndet = 540;
778 TObjArray *obj=new TObjArray(ndet);
779 obj->SetOwner();
780 for(Int_t ii=0; ii<ndet; ii++){
781 obj->Add(new TVectorD(AliTRDseedV1::kNtb));
782 }
783
784 //ibin = binlowedge of bin(ibin+1) = the number fills this bin
785 for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
786 const Double_t stat = hnor->GetBinContent(ibin+1);
787 if(stat<EPSILON){
788 continue;
789 }
790
791 const Int_t idet = ToDetector(ibin);
792 const Int_t itb = ToTimeBin(ibin);
793 TVectorD *vec=(TVectorD *)obj->At(idet);
794 (*vec)[itb] = stat;
795 }
796
797 hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
798 for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
799 const Int_t idet = ToDetector(ibin);
800 const TVectorD *vec=(TVectorD *)obj->At(idet);
801
802 Int_t nonzero = 0;
803 for(Int_t ii=0; ii<vec->GetNrows(); ii++){
804 if((*vec)[ii]>EPSILON){
805 nonzero++;
806 }
807 }
808
809 Double_t mean = 0;
810 const Double_t lowfrac = 0.6;
811 //if there are too many 0's, reject this chamber by settig mean=rms=0
812 if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
813 //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
814 mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
815 }
816
817 hmean->SetBinContent(ibin+1, mean);
818 }
819
820 delete obj;
821}
822
823void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run)
824{
825 //
826 //produce calibration objects
827 //
828
829 TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd ");
830 for(Int_t iter=0; iter<8; iter++){
831 objnames+= GetPHQName(0, iter)+" ";
832 }
833
834 TList * lout = new TList;
835 lout->SetOwner();
836
837 TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
838
839 const Int_t nh=lin->GetEntries();
840 for(Int_t ii=0; ii<nh; ii++){
841 const THnSparseD *hh=(THnSparseD*)lin->At(ii);
842 const TString hname = hh->GetName();
843 if(!objnames.Contains(hname))
844 continue;
845
846 TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream);
847 lout->Add(cobj0);
848 }
849
850 //lout->ls();
851
852 //=============================================================
853 //=============================================================
854
855 TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
856 fout->cd();
857 const Int_t nout=lout->GetEntries();
858 for(Int_t ii=0; ii<nout; ii++){
859 const TString oname = lout->At(ii)->GetName();
860 if(oname.Contains("Obj")){
861 TObjArray * cobj = (TObjArray*) lout->At(ii);
862 cobj->Write(oname, TObjArray::kSingleKey);
863 }
864 }
865 fout->Save();
866 fout->Close();
867 delete fout;
868
869 fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
870 fout->cd();
871 lin->Write();
872 lout->Write();
873 fout->Save();
874 fout->Close();
875 delete fout;
876
877 delete calibStream;
878
879 /*
880 http://root.cern.ch/root/html/TH1.html
881 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:
882
883 h->SetDirectory(0); for the current histogram h
884 TH1::AddDirectory(kFALSE); sets a global switch disabling the reference
885
886 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.
887 */
888 delete lout;
889}
890
891TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnSparseD *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
892{
893 //
894 //produce calibration objects
895 //
896
897 const TString hname = hh->GetName();
898 const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
899
900 //----------------------------------------
901 // Define nbin, tag, cobj0
902 //----------------------------------------
903 Int_t nbin =-999;
904 if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){
905 nbin = NTRDchamber();
906 }
907 else if(hname.Contains("PHQ")){
908 nbin = NTRDtimebin();
909 }
910 else{
911 printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1);
912 }
913
914 TString tag(hname);
915 tag.ReplaceAll("Hist","Obj");
916
917 TObjArray * cobj0 = new TObjArray(1);
918 cobj0->SetOwner();
919 cobj0->SetName(tag);
920 cobj0->Add(new TVectorD(nbin));
921
922 //----------------------------------------
923 // Define lowFrac, highFrac
924 //----------------------------------------
925 Double_t lowFrac = -999, highFrac = -999;
926 if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){
927 lowFrac = 0.01; highFrac = Q0Frac();
928 }
929 else if(hname.Contains("PHQ") && kinvq){
930 lowFrac = Q1Frac(); highFrac = 0.99;
931 }
932 else{
933 lowFrac = 0.01;
934 highFrac = 0.99;
935 }
936
937 //----------------------------------------
938 // Get analysis result
939 //----------------------------------------
940 TH1::AddDirectory(kFALSE);//important!
941 TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
942 TH2D *hpj = hh->Projection(1,0);
943 FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
944 if(hname.Contains("PHQ")){
945 GetPHCountMeanRMS(hnor, htrdphmean);
946 if(lout){
947 lout->Add(htrdphmean);
948 }
949 }
950 delete hpj;
951
952 if(lout){
953 lout->Add(hnor);
954 lout->Add(hmpv);
955 lout->Add(hwid);
956 lout->Add(hres);
957 }
958
959 //----------------------------------------
960 // Define Counter
961 //----------------------------------------
962 TVectorD *countDet=0x0;
963 TObjArray *countSSL=0x0;
964
965 if(hname.Contains("PHQ") && !kinvq){
966 countDet=new TVectorD(540);
967 countSSL=new TObjArray(90);//SectorStackLayer
968 countSSL->SetOwner();
969 for(Int_t ii=0; ii<90; ii++){
970 countSSL->Add(new TVectorD(6));
971 }
972 }
973
974 //----------------------------------------
975 // Fill cobj0
976 //----------------------------------------
977
978 //ibin = binlowedge of bin(ibin+1) = the number fills this bin
979 for(Int_t ibin=0; ibin<nbin; ibin++){
980 Double_t gnor = hnor->GetBinContent(ibin+1);
981 Double_t gmpv = hmpv->GetBinContent(ibin+1);
982 Double_t gwid = hwid->GetBinContent(ibin+1);
983 Double_t gres = hres->GetBinContent(ibin+1);
984
985 //--- set additional cut by kpass
986 Bool_t kpass = kTRUE;
987 Double_t gtrdphmean = -999;
988 if(htrdphmean){
989 gtrdphmean = htrdphmean->GetBinContent(ibin+1);
990 //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
991 if(gtrdphmean<EPSILON){
992 kpass = kFALSE;
993 }
994 if(gnor<TimeBinCountCut()*gtrdphmean){
995 kpass = kFALSE;
996 }
997 }
998
999 //--- set calibration constant p0
1000 Double_t p0= 0;
1001
1002 //reason for gmpv=0:
1003 //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;
1004 //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
1005
1006 if(gmpv>EPSILON && kpass){
1007 if(tag.Contains("T0")){
1008 p0 = gmpv;
1009 }
1010 else{
1011 p0 = 1/gmpv;
1012 }
1013 //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
1014 }
1015
1016 (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
1017
1018 //--- save optional record
1019 if(p0>EPSILON && countDet && countSSL){
1020 const Int_t idet = ToDetector(ibin);
1021 (*countDet)[idet]=1;
1022
1023 const Int_t isector = ToSector(ibin);
1024 const Int_t istack = ToStack(ibin);
1025 const Int_t ilayer = ToLayer(ibin);
1026 TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
1027 (*vecsectorstack)[ilayer]=1;
1028 }
1029
1030 if(calibStream){
1031 (*calibStream)<<tag<<
1032 "run="<<run<<
1033 "p0="<<p0<<
1034
1035 "nor="<<gnor<<
1036 "mpv="<<gmpv<<
1037 "wid="<<gwid<<
1038 "res="<<gres<<
1039 "gtrdphmean="<<gtrdphmean<<
1040
1041 "ibin="<<ibin<<
1042 "\n";
1043 }
1044 }
1045
1046 //----------------------------------------
1047 // Status Report
1048 //----------------------------------------
1049 if(countDet && countSSL){
1050 TVectorD count2Dstack(90);
1051 for(Int_t ii=0; ii<90; ii++){
1052 TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
1053 const Int_t nlayer = (Int_t)vecsectorstack->Sum();
1054 if(nlayer==6){
1055 count2Dstack[ii]=1;
1056 }
1057 }
1058
1059 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());
1060 }
1061
1062 //----------------------------------------
1063 // Clean Up
1064 //----------------------------------------
1065
1066 TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
1067 const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
1068 for(Int_t ihh=0; ihh<nhh; ihh++){
1069 if(!lout){
1070 delete (*hhs[ihh]);
1071 }
1072 }
1073
1074 delete countDet;
1075 delete countSSL;
1076
1077 //----------------------------------------
1078
1079 return cobj0;
1080}
1081
1082//===================================================================================
1083// dEdx calculation
1084//===================================================================================
1085Double_t AliTRDdEdxUtils::GetSignal(const Int_t nch, const Int_t ncls, const Double_t qq)
1086{
1087 //
1088 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1089 //
1090 if(ncls>1e3){
1091 printf("AliTRDdEdxUtils::GetSignal error ncls>1e3! %d\n", ncls); exit(1);
1092 }
1093
1094 return nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq);
1095}
1096
1097Int_t AliTRDdEdxUtils::GetNch(const Double_t signal)
1098{
1099 //
1100 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1101 //
1102 return Int_t(signal/1e6);
1103
1104}
1105
1106Int_t AliTRDdEdxUtils::GetNcls(const Double_t signal)
1107{
1108 //
1109 //signal = Nch*1e6 + Ncls*1e3 + (Q>=1e3? 999 : Q)
1110 //
1111
1112 return Int_t ( (signal-GetNch(signal)*1e6)/1e3 );
1113}
1114
1115Double_t AliTRDdEdxUtils::GetQ(const Double_t signal)
1116{
1117 //
1118 //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1119 //
1120
1121 return signal-GetNch(signal)*1e6 - GetNcls(signal)*1e3;
1122}
1123
1124Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
1125{
1126 //
1127 //template for cookdedx
1128 //
1129 if(cobj){
1130 if(arrayQ && arrayX){
1131 ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
1132 }
1133 else{
1134 printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
1135 }
1136 }
1137
1138 Double_t lowFrac =-999, highFrac = -999;
1139 if(kinvq){
1140 lowFrac = Q1Frac(); highFrac = 0.99;
1141 }
1142 else{
1143 lowFrac = 0.01; highFrac = Q0Frac();
1144 }
1145
1146 Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
1147 if(kinvq){
1148 if(meanQ>EPSILON){
1149 meanQ = 1/meanQ;
1150 }
1151 }
1152
1153 return meanQ;
1154}
1155
1156Double_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)
1157{
1158 //
1159 //combine tpc and trd dedx
1160 //
1161
1162 for(Int_t iq=0; iq<tpcncls; iq++){
1163 (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
1164 if(tpcarrayX && trdarrayX && coarrayX){
1165 (*coarrayX)[iq]=(*tpcarrayX)[iq];
1166 }
1167 }
1168 for(Int_t iq=0; iq<trdncls; iq++){
1169 (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
1170 if(tpcarrayX && trdarrayX && coarrayX){
1171 (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
1172 }
1173 }
1174
1175 concls=trdncls+tpcncls;
1176
1177 const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
1178
1179 return coQ;
1180}
1181
1182
1183//===================================================================================
1184// dEdx Getter and Setter
1185//===================================================================================
1186Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
1187{
1188 //
1189 //return angular normalization factor
1190 //
1191
1192 return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
1193}
1194
1195Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
1196{
1197 //
1198 //get cluster charge
1199 //
1200 Double_t dq = 0;
1201 const AliTRDcluster *cl = 0x0;
1202
1203 cl = seed->GetClusters(itb); if(cl) dq+=cl->GetRawQ();
1204 cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl) dq+=cl->GetRawQ();
1205
1206 dq /= GetAngularCorrection(seed);
1207
1208 dq /= 45.;
1209
1210 if(kinvq){
1211 if(dq>EPSILON){
1212 dq = 1/dq;
1213 }
1214 }
1215
1216 return dq;
1217}
1218
1219Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
1220{
1221 //
1222 //return nclustter
1223 //(if kinvq, return 1/q array), size of array must be larger than 31*6
1224 //
1225 if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1226 printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
1227 }
879b7730 1228 if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
fb3d369e 1229 printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
1230 }
1231
1232 const Int_t mintb = 0;
1233 const Int_t maxtb = AliTRDseedV1::kNtb-1;
1234 if(timeBin0<mintb) timeBin0=mintb;
1235 if(timeBin1>maxtb) timeBin1=maxtb;
1236 if(tstep<=0) tstep=1;
1237
1238 //============
1239 Int_t tbN=0;
1240 Double_t tbQ[200];
1241 Int_t tbBin[200];
1242
1243 for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1244 const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1245 if(!seed)
1246 continue;
1247
1248 const Int_t det = seed->GetDetector();
1249
1250 for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
1251 const Double_t dq = GetClusterQ(kinvq, seed, itb);
1252 if(dq<EPSILON)
1253 continue;
1254
1255 const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
1256
1257 tbQ[tbN]=dq;
1258 tbBin[tbN]=gtb;
1259 tbN++;
1260 }
1261 }
1262
1263 Int_t ncls = 0;
1264 for(Int_t iq=0; iq<tbN; iq++){
1265 if(tbQ[iq]<EPSILON)
1266 continue;
1267
1268 (*arrayQ)[ncls] = tbQ[iq];
1269 (*arrayX)[ncls] = tbBin[iq];
1270
1271 ncls++;
1272 }
1273
879b7730 1274 static Int_t kprint = 100;
fb3d369e 1275 if(kprint<0){
1276 printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n");
1277 for(Int_t iq=0; iq<ncls; iq++){
1278 const Int_t ichamber = ToLayer((*arrayX)[iq]);
1279 const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1280 if(!seed){
1281 printf("error seed null!!\n"); exit(1);
1282 }
1283 const Double_t rawq = (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
1284 printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
1285 }
1286 printf("\n");
1287 }
879b7730 1288 kprint++;
fb3d369e 1289
1290 return ncls;
1291}
1292
1293Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
1294{
1295 //
1296 //arrayX det*Ntb+itb -> itb
1297 //
1298
1299 TVectorD countChamber(6);
1300 for(Int_t ii = 0; ii<ncls; ii++){
1301 const Int_t xx = (Int_t)(*arrayX)[ii];
1302 const Int_t idet = ToDetector(xx);
1303
1304 const Double_t ich = AliTRDgeometry::GetLayer(idet);
1305 const Double_t itb = ToTimeBin(xx);
1306 (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
1307
1308 countChamber[ich] = 1;
1309 }
1310
1311 const Double_t nch = countChamber.Sum();
1312 return (Int_t) nch;
1313}
1314
1315void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnSparseD * hgain, THnSparseD * ht0, THnSparseD * hvd)
1316{
1317 //
1318 //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution
1319 //
1320
1321 static Int_t kprint = 100;
1322
1323 fgNchamber = 0;
1324 for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1325 //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities
1326 fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0;
1327
1328 const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber);
1329 if (!seed)
1330 continue;
1331
1332 const Int_t idet = seed->GetDetector();
1333
1334 //-------------------------------------------------------------------------
1335
1336 Double_t qsum = 0, qtsum = 0, w2sum = 0;
1337 for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
1338 const Double_t dq = GetClusterQ(0, seed, itb);
1339 if(dq<EPSILON)
1340 continue;
1341
1342 qsum += dq;
1343 qtsum += dq*itb;
1344 w2sum += dq*itb*itb;
1345 }
1346 if(qsum<EPSILON)
1347 continue;
1348
1349 //-------------------------------------------------------------------------
1350
1351 Double_t tbm, tbr = 0;
1352 tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr);
1353
1354 qsum /= 1.25e3/45.;
1355
1356 if(hgain){
1357 const Double_t var[]={idet, qsum};
1358 hgain->Fill(var);
1359 }
1360 if(ht0){
1361 const Double_t var[]={idet, tbm};
1362 ht0->Fill(var);
1363 }
1364 if(hvd){
1365 const Double_t var[]={idet, tbr};
1366 hvd->Fill(var);
1367 }
1368
1369 Double_t gain = 1, t0 = 0, vd = 1;
1370 if(kcalib){
1371 if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);}
1372 if(! fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0 array null!!\n"); exit(1);}
1373 if(! fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd array null!!\n"); exit(1);}
1374
1375 const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet];
1376 const TVectorD * t0vec = (TVectorD*) fgObjT0->At(0); t0 = (* t0vec)[idet];
1377 const TVectorD * vdvec = (TVectorD*) fgObjVd->At(0); vd = (* vdvec)[idet];
1378 }
1379 if(kprint<0){
1380 printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1381 printf("idet = %d;\n", idet);
1382 printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd);
1383 printf("\n");
1384 }
1385
1386 qsum *= gain;
1387 tbm = (tbm-t0)*vd;
1388
1389 if(qsum<EPSILON)
1390 continue;
1391
1392 //-------------------------------------------------------------------------
1393
1394 //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0)
1395 fgChamberQ[ichamber] = qsum;
1396 fgChamberTmean[ichamber] = tbm;
1397 fgNchamber++;
1398 }
1399
1400 if(kprint<0){
1401 printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n");
1402
1403 printf("\nfgNchamber = %d\n", fgNchamber);
1404 for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1405 printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]);
1406 }
1407 }
1408
1409 fgTrackTmean = -999;
1410 if(fgNchamber){
1411 fgTrackTmean = 0;
1412 for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1413 fgTrackTmean += fgChamberTmean[ich];
1414 }
1415 fgTrackTmean /= fgNchamber;
1416 }
1417
1418 if(kprint<0){
1419 printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1420 printf("GetTrackTmean() %15f\n", GetTrackTmean());
1421 printf("\n");
1422 }
1423 kprint++;
1424
1425 return;
1426}
1427
1428
1429//===================================================================================
1430// dEdx Parameterization
1431//===================================================================================
1432
1433Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg)
1434{
1435 //
1436 //truncated Mean Q_{xx} in TRD
1437 //
1438
1439 Double_t par[8];
1440 //03132012161150
1441 //opt: ppQ0
1442par[0]= 2.397001e-01;
1443par[1]= 1.334697e+00;
1444par[2]= 6.967470e+00;
1445par[3]= 9.055289e-02;
1446par[4]= 9.388760e+00;
1447par[5]= 9.452742e-04;
1448par[6]= -1.866091e+00;
1449par[7]= 1.403545e+00;
1450
1451 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale 0.428543 at ltbg 0.650000
1452 return 0.428543*MeandEdxTR(&bg, par);
1453}
1454
1455Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg)
1456{
1457 //
1458 //truncated Mean Q_{xx} in TRD
1459 //
1460
1461 Double_t par[8];
1462
1463 //03132012161259
1464 //opt: PbPbQ0
1465par[0]= 1.844912e-01;
1466par[1]= 2.509702e+00;
1467par[2]= 6.744031e+00;
1468par[3]= 7.355123e-02;
1469par[4]= 1.166023e+01;
1470par[5]= 1.736186e-04;
1471par[6]= -1.716063e+00;
1472par[7]= 1.611366e+00;
1473
1474 ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale 0.460994 at ltbg 0.650000
1475 return 0.460994*MeandEdxTR(&bg, par);
1476}
1477
1478Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg)
1479{
1480 //
1481 //truncated Mean 1/(1/Q)_{xx} in TRD
1482 //
1483
1484 Double_t par[8];
1485
1486 //So 4. Mär 13:30:51 CET 2012
1487 //opt: trdppQ1
1488 par[0]= 2.434646e-01;
1489 par[1]= 1.400211e+00;
1490 par[2]= 6.937471e+00;
1491 par[3]= 7.758118e-02;
1492 par[4]= 1.097372e+01;
1493 par[5]= 4.297518e-04;
1494 par[6]= -1.806266e+00;
1495 par[7]= 1.543811e+00;
1496
1497 //hhtype2Q1b2c2 scale 0.418629 at ltbg 0.650000
1498
1499 return 0.418629*MeandEdxTR(&bg, par);
1500}
1501
1502Double_t AliTRDdEdxUtils::Q1MeanTRDPbPb(const Double_t bg)
1503{
1504 //
1505 //truncated Mean 1/(1/Q)_{xx} in TRD
1506 //
1507
1508 Double_t par[8];
1509
1510 //So 4. Mär 13:30:52 CET 2012
1511 //opt: trdPbPbQ1
1512 par[0]= 2.193660e-01;
1513 par[1]= 2.051864e+00;
1514 par[2]= 6.825112e+00;
1515 par[3]= 6.151693e-02;
1516 par[4]= 1.390343e+01;
1517 par[5]= 6.010032e-05;
1518 par[6]= -1.676324e+00;
1519 par[7]= 1.838873e+00;
1520
1521 //hhtype4Q1b2c2 scale 0.457988 at ltbg 0.650000
1522
1523 return 0.457988*MeandEdxTR(&bg, par);
1524}
1525
1526Double_t AliTRDdEdxUtils::QMeanTPC(const Double_t bg)
1527{
1528 //
1529 //bethe bloch in TPC
1530 //
1531
1532 Double_t par[5];
1533 //Mi 15. Feb 14:48:05 CET 2012
1534 //train_2012-02-13_1214.12001, tpcsig
1535 par[0]= 4.401269;
1536 par[1]= 9.725370;
1537 par[2]= 0.000178;
1538 par[3]= 1.904962;
1539 par[4]= 1.426576;
1540
1541 return MeandEdx(&bg, par);
1542}
1543
1544Double_t AliTRDdEdxUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
1545{
1546 //
1547 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1548 //npar = 8 = 3+5
1549 //
1550 Double_t ptr[4]={0};
1551 for(int ii=0; ii<3; ii++){
1552 ptr[ii+1]=pin[ii];
1553 }
1554 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
1555}
1556
1557Double_t AliTRDdEdxUtils::MeanTR(const Double_t * xx, const Double_t * par)
1558{
1559 //
1560 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1561 //npar = 4
1562 //
1563
1564 const Double_t bg = xx[0];
1565 const Double_t gamma = sqrt(1+bg*bg);
1566
1567 const Double_t p0 = TMath::Abs(par[1]);
1568 const Double_t p1 = TMath::Abs(par[2]);
1569 const Double_t p2 = TMath::Abs(par[3]);
1570
1571 const Double_t zz = TMath::Log(gamma);
1572 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
1573
1574 return par[0]+tryield;
1575}
1576
1577Double_t AliTRDdEdxUtils::MeandEdx(const Double_t * xx, const Double_t * par)
1578{
1579 //
1580 //ALEPH parametrization for dEdx
1581 //npar = 5
1582 //
1583
1584 const Double_t bg = xx[0];
1585 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
1586
1587 const Double_t p0 = TMath::Abs(par[0]);
1588 const Double_t p1 = TMath::Abs(par[1]);
1589 const Double_t p2 = TMath::Abs(par[2]);
1590 const Double_t p3 = TMath::Abs(par[3]);
1591 const Double_t p4 = TMath::Abs(par[4]);
1592
1593 const Double_t aa = TMath::Power(beta, p3);
1594 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
1595
1596 //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
1597
1598 return (p1-aa-bb)*p0/aa;
1599}
1600
1601Double_t AliTRDdEdxUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
1602{
1603 //
1604 //f(x)-> f(y) with y=log10(x)
1605 //
1606 const Double_t x2[]={TMath::Power(10, xx[0])};
1607 return func(x2, par);
1608}
1609
1610//===================================================================================
1611// Detector, Data and Control Constant
1612//===================================================================================
1613Int_t AliTRDdEdxUtils::ToDetector(const Int_t gtb)
1614{
1615 //
1616 //gtb = det*Ntb+itb
1617 //
1618 return gtb/AliTRDseedV1::kNtb;
1619}
1620
1621Int_t AliTRDdEdxUtils::ToTimeBin(const Int_t gtb)
1622{
1623 //
1624 //gtb = det*Ntb+itb
1625 //
1626 return gtb%AliTRDseedV1::kNtb;
1627}
1628
1629Int_t AliTRDdEdxUtils::ToSector(const Int_t gtb)
1630{
1631 //
1632 //return sector
1633 //
1634 return AliTRDgeometry::GetSector(ToDetector(gtb));
1635}
1636
1637Int_t AliTRDdEdxUtils::ToStack(const Int_t gtb)
1638{
1639 //
1640 //return stack
1641 //
1642 return AliTRDgeometry::GetStack(ToDetector(gtb));
1643}
1644
1645Int_t AliTRDdEdxUtils::ToLayer(const Int_t gtb)
1646{
1647 //
1648 //return layer
1649 //
1650 return AliTRDgeometry::GetLayer(ToDetector(gtb));
1651}
1652
1653TString AliTRDdEdxUtils::GetRunType(const Int_t run)
1654{
1655 //
1656 //return run type
1657 //
1658
1659 TString type;
1660 if(run>=121527 && run<= 126460)//LHC10d
1661 type="pp2010LHC10d";
1662 else if(run>=126461 && run<= 130930)//LHC10e
1663 type="pp2010LHC10e";
1664 else if(run>=136782 && run <= 139846)//LHC10h
1665 type="PbPb2010LHC10h";
1666 else if(run>= 143224 && run<= 143237)//2011Feb
1667 type="cosmic2011Feb";
1668 else if(run>= 150587 && run<= 154930){
1669 type="cosmic2011MayJun";
1670
1671 TString runstr(Form("%d",run));
1672 const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
1673 if(listrun1kg.Contains(runstr)){
1674 type+="1kG";
1675 }
1676 else{
1677 type+="5kG";
1678 }
1679 }
1680 else{
1681 type="unknown";
1682 }
1683
1684 type.ToUpper();
1685 return type;
1686}
1687
1688void AliTRDdEdxUtils::PrintControl()
1689{
1690 //
1691 //print out control variable
1692 //
1693 printf("\nAliTRDdEdxUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn());
1694}
1695//===================================================================================
1696//===================================================================================