Commented unnecessary include AliLog
[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-12
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", 0, 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     }
212     return kFALSE;
213   }
214   
215   delete fgHistArray; 
216   fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
217
218   return kTRUE;
219 }
220
221 void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale)
222 {
223   //
224   //fill calibration hist
225   //
226   if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
227
228   for(Int_t ii=0; ii<ncls; ii++){
229     Double_t dq = (*arrayQ)[ii]/scale;
230     const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
231     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
232
233     if(dq<qmin)
234       dq = qmin;
235     if(dq>qmax)
236       dq = qmax;
237
238     const Double_t xx = (*arrayX)[ii];
239     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
240     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
241
242     if(xx>=xmax || xx<xmin){
243       printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
244     }
245
246     const Double_t var[]={xx, dq};
247     hcalib->Fill(var);
248   }
249 }
250
251 void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale) 
252 {
253   //
254   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
255   //
256   if(!fgHistArray){
257     printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
258   }
259
260   THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
261
262   TVectorD arrayQ(200), arrayX(200);
263   const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
264   FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
265
266   static Int_t kprint = 100;
267   if(kprint<0){
268     printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
269     printf("\nkinvq= %d;\n", kinvq);
270     for(Int_t iq=0; iq<ncls; iq++){
271       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
272     }
273     printf("\n");
274   }
275   kprint++;
276 }
277
278 //============================================================
279 //
280 //============================================================
281
282 Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
283 {
284   //
285   //the scale used in calibration
286   //
287
288   if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
289     return -999;
290
291   if(tpcsig<EPSILON)
292     return -999;
293
294   return tpcsig/45;
295
296 }
297
298 void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
299 {
300   //
301   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
302   //
303   const Int_t ndet = 540;
304   TObjArray *obj=new TObjArray(ndet);
305   obj->SetOwner();
306   for(Int_t ii=0; ii<ndet; ii++){
307     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
308   }
309
310   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
311   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
312     const Double_t stat = hnor->GetBinContent(ibin+1);
313     if(stat<EPSILON){
314       continue;
315     }
316
317     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
318     const Int_t itb  = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
319     TVectorD *vec=(TVectorD *)obj->At(idet);
320     (*vec)[itb] = stat;
321   }
322
323   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
324   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
325     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
326     const TVectorD *vec=(TVectorD *)obj->At(idet);
327
328     Int_t nonzero = 0;
329     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
330       if((*vec)[ii]>EPSILON){
331         nonzero++;
332       }
333     }
334
335     Double_t mean = 0;
336     const Double_t lowfrac = 0.6;
337     //if there are too many 0's, reject this chamber by settig mean=rms=0
338     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
339       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
340       mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
341     }
342
343     hmean->SetBinContent(ibin+1, mean);
344   }
345
346   delete obj;
347 }
348
349   /*
350     http://root.cern.ch/root/html/TH1.html
351     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:
352     
353     h->SetDirectory(0);          for the current histogram h
354     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
355     
356     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. 
357   */
358
359 TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
360 {
361   //
362   //produce calibration objects
363   //
364
365   const TString hname = hh->GetName();
366   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
367
368   //----------------------------------------
369   //               Define nbin, tag, cobj0
370   //----------------------------------------
371   const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
372
373     
374   TString tag(hname);
375   tag.ReplaceAll("Hist","Obj");
376   
377   TObjArray * cobj0 = new TObjArray(1);
378   cobj0->SetOwner();
379   cobj0->SetName(tag);
380   cobj0->Add(new TVectorD(nbin));
381   
382   //----------------------------------------
383   //               Define lowFrac, highFrac
384   //----------------------------------------
385   Double_t lowFrac = -999, highFrac = -999;
386   if(!kinvq) {
387     lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
388   }
389   else{
390     lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
391   }
392   
393   //----------------------------------------
394   //               Get analysis result
395   //----------------------------------------
396   TH1::AddDirectory(kFALSE);//important!
397   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
398   TH2D *hpj = hh->Projection(1,0);
399   AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
400   GetPHCountMeanRMS(hnor, htrdphmean);
401   if(lout){
402     lout->Add(htrdphmean);
403   }
404   delete hpj;
405   
406   if(lout){
407     lout->Add(hnor);
408     lout->Add(hmpv);
409     lout->Add(hwid);
410     lout->Add(hres);
411   }
412   
413   //----------------------------------------
414   //               Define Counter
415   //----------------------------------------
416   TVectorD *countDet=0x0;
417   TObjArray *countSSL=0x0;
418
419   if(!kinvq){
420     countDet=new TVectorD(540);
421     countSSL=new TObjArray(90);//SectorStackLayer
422     countSSL->SetOwner();
423     for(Int_t ii=0; ii<90; ii++){
424       countSSL->Add(new TVectorD(6));
425     }
426   }
427
428   //----------------------------------------
429   //               Fill cobj0
430   //----------------------------------------
431
432   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
433   for(Int_t ibin=0; ibin<nbin; ibin++){
434     Double_t gnor = hnor->GetBinContent(ibin+1);
435     Double_t gmpv = hmpv->GetBinContent(ibin+1);
436     Double_t gwid = hwid->GetBinContent(ibin+1);
437     Double_t gres = hres->GetBinContent(ibin+1);
438
439     //--- set additional cut by kpass
440     Bool_t kpass = kTRUE;
441     Double_t gtrdphmean = -999;
442     if(htrdphmean){
443       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
444       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
445       if(gtrdphmean<EPSILON){
446         kpass = kFALSE;
447       }
448       if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
449         kpass = kFALSE;
450       }
451     }
452     
453     //--- set calibration constant p0
454     Double_t p0= 0;
455     
456     //reason for gmpv=0:
457     //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;
458     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
459     
460     if(gmpv>EPSILON && kpass){ 
461       if(tag.Contains("T0")){
462         p0 = gmpv;
463       }
464       else{
465         p0 = 1/gmpv;
466       }
467       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
468     }
469
470     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
471
472     //--- save optional record
473     if(p0>EPSILON && countDet && countSSL){
474       const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
475       (*countDet)[idet]=1;
476       
477       const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
478       const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
479       const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
480       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
481       (*vecsectorstack)[ilayer]=1;
482     }
483     
484     if(calibStream){
485       (*calibStream)<<tag<<
486         "run="<<run<<
487         "p0="<<p0<<
488         
489         "nor="<<gnor<<
490         "mpv="<<gmpv<<
491         "wid="<<gwid<<
492         "res="<<gres<<
493         "gtrdphmean="<<gtrdphmean<<
494         
495         "ibin="<<ibin<<
496         "\n";
497     }
498   }
499   
500   //----------------------------------------
501   //               Status Report
502   //----------------------------------------
503   if(countDet && countSSL){
504     TVectorD count2Dstack(90);
505     for(Int_t ii=0; ii<90; ii++){
506       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
507       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
508       if(nlayer==6){
509         count2Dstack[ii]=1;
510       }
511     }
512
513     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());
514   }
515
516   //----------------------------------------
517   //               Clean Up
518   //----------------------------------------
519   
520   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
521   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
522   for(Int_t ihh=0; ihh<nhh; ihh++){
523     if(!lout){
524       delete (*hhs[ihh]);
525     }
526   }
527   
528   delete countDet;
529   delete countSSL;
530
531   //----------------------------------------
532
533   return cobj0;
534 }
535