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