Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxCalibUtils.cxx
CommitLineData
6951a056 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
6fa94550 60#define EPSILON 1e-8
6951a056 61
62AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
63TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;
64
65//============================================================
66// CalibObj
67//============================================================
68TObjArray * 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
81TObjArray * 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
93void AliTRDdEdxCalibUtils::DeleteObjArray()
94{
95 //
96 //delete calib obj
97 //
98 if(fgObjArray){
99 fgObjArray->SetOwner();
100 delete fgObjArray;
101 fgObjArray = 0x0;
102 }
103}
104
687aa844 105Bool_t AliTRDdEdxCalibUtils::GenerateOCDB(const Int_t run, const TString path)
6951a056 106{
107 //
687aa844 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
6951a056 112 //
113
687aa844 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);
6951a056 136 }
6951a056 137 }
138
139 AliCDBMetaData *metaData= new AliCDBMetaData();
140 metaData->SetObjectClassName("TObjArray");
141 metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
142
fe829aa0 143 AliCDBId id1("TRD/Calib/PHQ", run<0? 0 : run , run<0 ? 999999999 : run);
6951a056 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//============================================================
156void AliTRDdEdxCalibUtils::DeleteHistArray()
157{
158 //
159 //delete calib hist
160 //
161 delete fgHistArray;
162 fgHistArray = 0x0;
163}
164
165THnBase * 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
176void 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
187Bool_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();
177f532f 211 delete array;
6951a056 212 }
213 return kFALSE;
214 }
215
216 delete fgHistArray;
217 fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
218
219 return kTRUE;
220}
221
222void 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++){
6856b6a8 230 Double_t dq = (*arrayQ)[ii]/scale;
231 const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
6951a056 232 const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
6856b6a8 233
234 if(dq<qmin)
235 dq = qmin;
236 if(dq>qmax)
237 dq = qmax;
238
239 const Double_t xx = (*arrayX)[ii];
6951a056 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
6856b6a8 247 const Double_t var[]={xx, dq};
6951a056 248 hcalib->Fill(var);
249 }
250}
251
252void 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
283Double_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
687aa844 295 return tpcsig/45;
6951a056 296
297}
298
299void 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
6951a056 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 */
6951a056 359
360TObjArray* 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