Bug fix (Xianguo)
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxUtils.cxx
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
57 THnSparseD * AliTRDdEdxUtils::fgHistGain=0x0;
58 THnSparseD * AliTRDdEdxUtils::fgHistT0=0x0;
59 THnSparseD * AliTRDdEdxUtils::fgHistVd=0x0;
60 TObjArray * AliTRDdEdxUtils::fgHistPHQ=new TObjArray(8);
61
62 TString AliTRDdEdxUtils::fgCalibFile;
63 TObjArray * AliTRDdEdxUtils::fgObjGain = 0x0;
64 TObjArray * AliTRDdEdxUtils::fgObjT0 = 0x0;
65 TObjArray * AliTRDdEdxUtils::fgObjVd = 0x0;
66 TObjArray * AliTRDdEdxUtils::fgObjPHQ = 0x0;
67
68 Int_t AliTRDdEdxUtils::fgNchamber = -999;
69 Double_t AliTRDdEdxUtils::fgChamberQ[6];
70 Double_t AliTRDdEdxUtils::fgChamberTmean[6];
71
72 Double_t AliTRDdEdxUtils::fgTrackTmean = -999;
73
74 //===================================================================================
75 //                                   Math and Histogram
76 //===================================================================================
77 void 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
114 Double_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
152 Double_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
195 Double_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
259 void 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
349 Int_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
362 AliTRDtrackV1 * 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
384 Bool_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
412 Bool_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 //===================================================================================
443 Double_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
459 Int_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
467 TObjArray * 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
480 TObjArray * 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
492 THnSparseD * 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
500 TString 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
508 void 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
528 Bool_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
566 void 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
617 void 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
635 void 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
668 Bool_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
705 void 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
729 void 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
753 Int_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
785 void 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
836 void 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
904 TObjArray* 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 //===================================================================================
1098 Double_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
1110 Int_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
1119 Int_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
1128 Double_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
1137 Double_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
1169 Double_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 //===================================================================================
1199 Double_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
1208 Double_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
1232 Int_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
1306 Int_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
1328 void 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
1446 Double_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
1455 par[0]=   2.397001e-01;
1456 par[1]=   1.334697e+00;
1457 par[2]=   6.967470e+00;
1458 par[3]=   9.055289e-02;
1459 par[4]=   9.388760e+00;
1460 par[5]=   9.452742e-04;
1461 par[6]=  -1.866091e+00;
1462 par[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
1468 Double_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
1478 par[0]=   1.844912e-01;
1479 par[1]=   2.509702e+00;
1480 par[2]=   6.744031e+00;
1481 par[3]=   7.355123e-02;
1482 par[4]=   1.166023e+01;
1483 par[5]=   1.736186e-04;
1484 par[6]=  -1.716063e+00;
1485 par[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
1491 Double_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   //So 4. Mär 13:30:51 CET 2012
1500   //opt: trdppQ1
1501   par[0]=   2.434646e-01;
1502   par[1]=   1.400211e+00;
1503   par[2]=   6.937471e+00;
1504   par[3]=   7.758118e-02;
1505   par[4]=   1.097372e+01;
1506   par[5]=   4.297518e-04;
1507   par[6]=  -1.806266e+00;
1508   par[7]=   1.543811e+00;
1509
1510   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
1511
1512   return  0.418629*MeandEdxTR(&bg, par);
1513 }
1514
1515 Double_t AliTRDdEdxUtils::Q1MeanTRDPbPb(const Double_t bg)
1516 {
1517   //
1518   //truncated Mean 1/(1/Q)_{xx} in TRD
1519   //
1520  
1521   Double_t par[8];
1522
1523   //So 4. Mär 13:30:52 CET 2012
1524   //opt: trdPbPbQ1
1525   par[0]=   2.193660e-01;
1526   par[1]=   2.051864e+00;
1527   par[2]=   6.825112e+00;
1528   par[3]=   6.151693e-02;
1529   par[4]=   1.390343e+01;
1530   par[5]=   6.010032e-05;
1531   par[6]=  -1.676324e+00;
1532   par[7]=   1.838873e+00;
1533
1534   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
1535
1536   return  0.457988*MeandEdxTR(&bg, par);
1537 }
1538
1539 Double_t AliTRDdEdxUtils::QMeanTPC(const Double_t bg)
1540 {
1541   //
1542   //bethe bloch in TPC
1543   //
1544
1545   Double_t par[5];
1546   //Mi 15. Feb 14:48:05 CET 2012
1547   //train_2012-02-13_1214.12001, tpcsig
1548   par[0]=       4.401269;
1549   par[1]=       9.725370;
1550   par[2]=       0.000178;
1551   par[3]=       1.904962;
1552   par[4]=       1.426576;
1553
1554   return MeandEdx(&bg, par);
1555 }
1556
1557 Double_t AliTRDdEdxUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
1558 {
1559   //
1560   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1561   //npar = 8 = 3+5
1562   //
1563   Double_t ptr[4]={0};
1564   for(int ii=0; ii<3; ii++){
1565     ptr[ii+1]=pin[ii];
1566   }
1567   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
1568 }
1569
1570 Double_t AliTRDdEdxUtils::MeanTR(const Double_t * xx, const Double_t * par)
1571 {
1572   //
1573   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1574   //npar = 4
1575   //
1576
1577   const Double_t bg = xx[0];
1578   const Double_t gamma = sqrt(1+bg*bg);
1579
1580   const Double_t p0 = TMath::Abs(par[1]);
1581   const Double_t p1 = TMath::Abs(par[2]);
1582   const Double_t p2 = TMath::Abs(par[3]);
1583
1584   const Double_t zz = TMath::Log(gamma);
1585   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
1586
1587   return par[0]+tryield;
1588 }
1589
1590 Double_t AliTRDdEdxUtils::MeandEdx(const Double_t * xx, const Double_t * par)
1591 {
1592   //
1593   //ALEPH parametrization for dEdx
1594   //npar = 5
1595   //
1596
1597   const Double_t bg = xx[0];
1598   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
1599
1600   const Double_t p0 = TMath::Abs(par[0]);
1601   const Double_t p1 = TMath::Abs(par[1]);
1602   const Double_t p2 = TMath::Abs(par[2]);
1603   const Double_t p3 = TMath::Abs(par[3]);
1604   const Double_t p4 = TMath::Abs(par[4]);
1605
1606   const Double_t aa = TMath::Power(beta, p3);
1607   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
1608
1609   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
1610
1611   return (p1-aa-bb)*p0/aa;
1612 }
1613
1614 Double_t AliTRDdEdxUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
1615 {
1616   //
1617   //f(x)-> f(y) with y=log10(x)
1618   //
1619   const Double_t x2[]={TMath::Power(10, xx[0])};
1620   return func(x2, par);
1621 }
1622
1623 //===================================================================================
1624 //                                Detector, Data and Control Constant 
1625 //===================================================================================
1626 Int_t  AliTRDdEdxUtils::ToDetector(const Int_t gtb)
1627 {
1628   //
1629   //gtb = det*Ntb+itb
1630   //
1631   return gtb/AliTRDseedV1::kNtb;
1632 }
1633
1634 Int_t  AliTRDdEdxUtils::ToTimeBin(const Int_t gtb)
1635
1636   //
1637   //gtb = det*Ntb+itb
1638   //
1639   return gtb%AliTRDseedV1::kNtb;
1640 }
1641
1642 Int_t  AliTRDdEdxUtils::ToSector(const Int_t gtb)
1643 {
1644   //
1645   //return sector
1646   //
1647   return AliTRDgeometry::GetSector(ToDetector(gtb));
1648 }
1649
1650 Int_t  AliTRDdEdxUtils::ToStack(const Int_t gtb)
1651 {
1652   //
1653   //return stack
1654   //
1655   return AliTRDgeometry::GetStack(ToDetector(gtb));
1656 }
1657
1658 Int_t  AliTRDdEdxUtils::ToLayer(const Int_t gtb)
1659 {
1660   //
1661   //return layer
1662   //
1663   return AliTRDgeometry::GetLayer(ToDetector(gtb));
1664 }
1665
1666 TString AliTRDdEdxUtils::GetRunType(const Int_t run)
1667 {
1668   //
1669   //return run type
1670   //
1671
1672   TString type;
1673   if(run>=121527 && run<= 126460)//LHC10d
1674     type="pp2010LHC10d";
1675   else if(run>=126461 && run<= 130930)//LHC10e
1676     type="pp2010LHC10e";
1677   else if(run>=136782 && run <= 139846)//LHC10h
1678     type="PbPb2010LHC10h";
1679   else if(run>= 143224 && run<= 143237)//2011Feb
1680     type="cosmic2011Feb";
1681   else if(run>= 150587 && run<= 154930){
1682     type="cosmic2011MayJun";
1683
1684     TString runstr(Form("%d",run));
1685     const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
1686     if(listrun1kg.Contains(runstr)){
1687       type+="1kG";
1688     }
1689     else{
1690       type+="5kG";
1691     }      
1692   }
1693   else{
1694     type="unknown";
1695   }
1696
1697   type.ToUpper();
1698   return type;
1699 }
1700
1701 void AliTRDdEdxUtils::PrintControl()
1702 {
1703   //
1704   //print out control variable
1705   //
1706   printf("\nAliTRDdEdxUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn());
1707 }
1708 //===================================================================================
1709 //===================================================================================