TOF + macros to each detector folder
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxCalibUtils.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 "AliTRDcalibDB.h"
52 #include "AliTRDCalROC.h"
53 #include "AliTRDtrackV1.h"
54
55 #include "AliTRDdEdxBaseUtils.h"
56 #include "AliTRDdEdxReconUtils.h"
57 #include "AliTRDdEdxCalibHistArray.h"
58 #include "AliTRDdEdxCalibUtils.h"
59
60 #define EPSILON 1e-8
61
62 AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
63 TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;
64
65 //============================================================
66 //                           CalibObj
67 //============================================================
68 TObjArray * AliTRDdEdxCalibUtils::GetObjArray()
69 {
70   //
71   //return fgObjArray, initialized if null
72   //
73
74   if(!fgObjArray){
75     fgObjArray = new TObjArray(8);
76   }
77
78   return fgObjArray;
79 }
80
81 TObjArray * AliTRDdEdxCalibUtils::GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge)
82 {
83   //
84   //return calib obj
85   //
86   if(!fgObjArray){
87     printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
88   }
89
90   return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
91 }
92
93 void AliTRDdEdxCalibUtils::DeleteObjArray()
94 {
95   //
96   //delete calib obj
97   //
98   if(fgObjArray){
99     fgObjArray->SetOwner();
100     delete fgObjArray;
101     fgObjArray = 0x0;
102   }
103 }
104
105 Bool_t AliTRDdEdxCalibUtils::GenerateOCDB(const Int_t run, const TString path)
106 {
107   //
108   //generate OCDB object PHQ, do like
109   //AliTRDdEdxCalibUtils::GenerateOCDB(run, "local://./")
110   //if fgObjArray==0x0, generate default one
111   //else generate according to fgObjArray
112   //
113
114   TObjArray * arr8 = 0x0;
115   if(fgObjArray){
116     arr8 = fgObjArray;
117   }
118   else{
119     arr8 = new TObjArray(8);
120     arr8->SetOwner();
121
122     for(Int_t ii=0; ii<8; ii++){
123       TObjArray * arr1 = new TObjArray(1);
124       arr1->SetOwner();
125       TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
126       objn.ReplaceAll("Hist","Obj");
127       arr1->SetName(objn);
128       
129       const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
130       TVectorD * vec = new TVectorD(nbins);
131       for(Int_t jj=0; jj<nbins; jj++){
132         (*vec)[jj] = 1;
133       }
134       arr1->Add(vec);
135       arr8->Add(arr1);
136     }
137   }
138
139   AliCDBMetaData *metaData= new AliCDBMetaData();
140   metaData->SetObjectClassName("TObjArray");
141   metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
142
143   AliCDBId id1("TRD/Calib/PHQ", run<0? 0 : run , run<0 ? 999999999 : run);
144   AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
145   gStorage->Put(arr8, id1, metaData);
146
147   delete metaData;
148   delete arr8;
149
150   return kTRUE;
151 }
152
153 //============================================================
154 //                           CalibHist
155 //============================================================
156 void AliTRDdEdxCalibUtils::DeleteHistArray()
157 {
158   //
159   //delete calib hist
160   //
161   delete fgHistArray;
162   fgHistArray = 0x0;
163 }
164
165 THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
166 {
167   //
168   //
169   //
170   if(iter<fgHistArray->GetSize())
171     return (THnBase*)fgHistArray->At(iter);
172   else
173     return 0x0;
174 }
175
176 void AliTRDdEdxCalibUtils::IniHistArray(TList *list, const Bool_t kNoInv)
177 {
178   //
179   //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
180   //
181
182   delete fgHistArray;
183   fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
184   list->Add(fgHistArray);
185 }
186
187 Bool_t AliTRDdEdxCalibUtils::ReadHistArray(const TString filename, const TString listname)
188 {
189   //
190   //used in AliTRDPreprocessorOffline
191   //read in calib hist from file, only for PHQ
192   //
193
194   //maybe already open by others... don't close
195   TFile fcalib(filename);
196   TObjArray * array = (TObjArray*)fcalib.Get(listname);
197   const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
198  
199   AliTRDdEdxCalibHistArray * tmph=0x0;
200   if(array){
201     tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
202   }
203   else{
204     tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
205   }
206   if(!tmph){
207     printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
208     fcalib.ls();
209     if(array){
210       array->ls();
211       delete array;
212     }
213     return kFALSE;
214   }
215   
216   delete fgHistArray; 
217   fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
218
219   return kTRUE;
220 }
221
222 void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale)
223 {
224   //
225   //fill calibration hist
226   //
227   if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
228
229   for(Int_t ii=0; ii<ncls; ii++){
230     Double_t dq = (*arrayQ)[ii]/scale;
231     const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
232     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
233
234     if(dq<qmin)
235       dq = qmin;
236     if(dq>qmax)
237       dq = qmax;
238
239     const Double_t xx = (*arrayX)[ii];
240     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
241     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
242
243     if(xx>=xmax || xx<xmin){
244       printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
245     }
246
247     const Double_t var[]={xx, dq};
248     hcalib->Fill(var);
249   }
250 }
251
252 void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale) 
253 {
254   //
255   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
256   //
257   if(!fgHistArray){
258     printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
259   }
260
261   THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
262
263   TVectorD arrayQ(200), arrayX(200);
264   const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
265   FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
266
267   static Int_t kprint = 100;
268   if(kprint<0){
269     printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
270     printf("\nkinvq= %d;\n", kinvq);
271     for(Int_t iq=0; iq<ncls; iq++){
272       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
273     }
274     printf("\n");
275   }
276   kprint++;
277 }
278
279 //============================================================
280 //
281 //============================================================
282
283 Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
284 {
285   //
286   //the scale used in calibration
287   //
288
289   if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
290     return -999;
291
292   if(tpcsig<EPSILON)
293     return -999;
294
295   return tpcsig/45;
296
297 }
298
299 void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
300 {
301   //
302   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
303   //
304   const Int_t ndet = 540;
305   TObjArray *obj=new TObjArray(ndet);
306   obj->SetOwner();
307   for(Int_t ii=0; ii<ndet; ii++){
308     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
309   }
310
311   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
312   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
313     const Double_t stat = hnor->GetBinContent(ibin+1);
314     if(stat<EPSILON){
315       continue;
316     }
317
318     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
319     const Int_t itb  = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
320     TVectorD *vec=(TVectorD *)obj->At(idet);
321     (*vec)[itb] = stat;
322   }
323
324   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
325   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
326     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
327     const TVectorD *vec=(TVectorD *)obj->At(idet);
328
329     Int_t nonzero = 0;
330     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
331       if((*vec)[ii]>EPSILON){
332         nonzero++;
333       }
334     }
335
336     Double_t mean = 0;
337     const Double_t lowfrac = 0.6;
338     //if there are too many 0's, reject this chamber by settig mean=rms=0
339     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
340       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
341       mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
342     }
343
344     hmean->SetBinContent(ibin+1, mean);
345   }
346
347   delete obj;
348 }
349
350   /*
351     http://root.cern.ch/root/html/TH1.html
352     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:
353     
354     h->SetDirectory(0);          for the current histogram h
355     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
356     
357     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. 
358   */
359
360 TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
361 {
362   //
363   //produce calibration objects
364   //
365
366   const TString hname = hh->GetName();
367   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
368
369   //----------------------------------------
370   //               Define nbin, tag, cobj0
371   //----------------------------------------
372   const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
373
374     
375   TString tag(hname);
376   tag.ReplaceAll("Hist","Obj");
377   
378   TObjArray * cobj0 = new TObjArray(1);
379   cobj0->SetOwner();
380   cobj0->SetName(tag);
381   cobj0->Add(new TVectorD(nbin));
382   
383   //----------------------------------------
384   //               Define lowFrac, highFrac
385   //----------------------------------------
386   Double_t lowFrac = -999, highFrac = -999;
387   if(!kinvq) {
388     lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
389   }
390   else{
391     lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
392   }
393   
394   //----------------------------------------
395   //               Get analysis result
396   //----------------------------------------
397   TH1::AddDirectory(kFALSE);//important!
398   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
399   TH2D *hpj = hh->Projection(1,0);
400   AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
401   GetPHCountMeanRMS(hnor, htrdphmean);
402   if(lout){
403     lout->Add(htrdphmean);
404   }
405   delete hpj;
406   
407   if(lout){
408     lout->Add(hnor);
409     lout->Add(hmpv);
410     lout->Add(hwid);
411     lout->Add(hres);
412   }
413   
414   //----------------------------------------
415   //               Define Counter
416   //----------------------------------------
417   TVectorD *countDet=0x0;
418   TObjArray *countSSL=0x0;
419
420   if(!kinvq){
421     countDet=new TVectorD(540);
422     countSSL=new TObjArray(90);//SectorStackLayer
423     countSSL->SetOwner();
424     for(Int_t ii=0; ii<90; ii++){
425       countSSL->Add(new TVectorD(6));
426     }
427   }
428
429   //----------------------------------------
430   //               Fill cobj0
431   //----------------------------------------
432
433   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
434   for(Int_t ibin=0; ibin<nbin; ibin++){
435     Double_t gnor = hnor->GetBinContent(ibin+1);
436     Double_t gmpv = hmpv->GetBinContent(ibin+1);
437     Double_t gwid = hwid->GetBinContent(ibin+1);
438     Double_t gres = hres->GetBinContent(ibin+1);
439
440     //--- set additional cut by kpass
441     Bool_t kpass = kTRUE;
442     Double_t gtrdphmean = -999;
443     if(htrdphmean){
444       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
445       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
446       if(gtrdphmean<EPSILON){
447         kpass = kFALSE;
448       }
449       if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
450         kpass = kFALSE;
451       }
452     }
453     
454     //--- set calibration constant p0
455     Double_t p0= 0;
456     
457     //reason for gmpv=0:
458     //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;
459     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
460     
461     if(gmpv>EPSILON && kpass){ 
462       if(tag.Contains("T0")){
463         p0 = gmpv;
464       }
465       else{
466         p0 = 1/gmpv;
467       }
468       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
469     }
470
471     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
472
473     //--- save optional record
474     if(p0>EPSILON && countDet && countSSL){
475       const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
476       (*countDet)[idet]=1;
477       
478       const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
479       const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
480       const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
481       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
482       (*vecsectorstack)[ilayer]=1;
483     }
484     
485     if(calibStream){
486       (*calibStream)<<tag<<
487         "run="<<run<<
488         "p0="<<p0<<
489         
490         "nor="<<gnor<<
491         "mpv="<<gmpv<<
492         "wid="<<gwid<<
493         "res="<<gres<<
494         "gtrdphmean="<<gtrdphmean<<
495         
496         "ibin="<<ibin<<
497         "\n";
498     }
499   }
500   
501   //----------------------------------------
502   //               Status Report
503   //----------------------------------------
504   if(countDet && countSSL){
505     TVectorD count2Dstack(90);
506     for(Int_t ii=0; ii<90; ii++){
507       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
508       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
509       if(nlayer==6){
510         count2Dstack[ii]=1;
511       }
512     }
513
514     printf("\nAliTRDdEdxCalibUtils::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());
515   }
516
517   //----------------------------------------
518   //               Clean Up
519   //----------------------------------------
520   
521   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
522   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
523   for(Int_t ihh=0; ihh<nhh; ihh++){
524     if(!lout){
525       delete (*hhs[ihh]);
526     }
527   }
528   
529   delete countDet;
530   delete countSSL;
531
532   //----------------------------------------
533
534   return cobj0;
535 }
536