Line data Source code
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 0 : if(!fgObjArray){
75 0 : fgObjArray = new TObjArray(8);
76 0 : }
77 :
78 0 : return fgObjArray;
79 0 : }
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 104 : if(!fgObjArray){
87 0 : printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
88 : }
89 :
90 104 : return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
91 : }
92 :
93 : void AliTRDdEdxCalibUtils::DeleteObjArray()
94 : {
95 : //
96 : //delete calib obj
97 : //
98 0 : if(fgObjArray){
99 0 : fgObjArray->SetOwner();
100 0 : delete fgObjArray;
101 0 : fgObjArray = 0x0;
102 0 : }
103 0 : }
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 0 : if(fgObjArray){
116 : arr8 = fgObjArray;
117 0 : }
118 : else{
119 0 : arr8 = new TObjArray(8);
120 0 : arr8->SetOwner();
121 :
122 0 : for(Int_t ii=0; ii<8; ii++){
123 0 : TObjArray * arr1 = new TObjArray(1);
124 0 : arr1->SetOwner();
125 0 : TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
126 0 : objn.ReplaceAll("Hist","Obj");
127 0 : arr1->SetName(objn);
128 :
129 0 : const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
130 0 : TVectorD * vec = new TVectorD(nbins);
131 0 : for(Int_t jj=0; jj<nbins; jj++){
132 0 : (*vec)[jj] = 1;
133 : }
134 0 : arr1->Add(vec);
135 0 : arr8->Add(arr1);
136 0 : }
137 : }
138 :
139 0 : AliCDBMetaData *metaData= new AliCDBMetaData();
140 0 : metaData->SetObjectClassName("TObjArray");
141 0 : metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
142 :
143 0 : AliCDBId id1("TRD/Calib/PHQ", run<0? 0 : run , run<0 ? 999999999 : run);
144 0 : AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
145 0 : gStorage->Put(arr8, id1, metaData);
146 :
147 0 : delete metaData;
148 0 : delete arr8;
149 :
150 : return kTRUE;
151 0 : }
152 :
153 : //============================================================
154 : // CalibHist
155 : //============================================================
156 : void AliTRDdEdxCalibUtils::DeleteHistArray()
157 : {
158 : //
159 : //delete calib hist
160 : //
161 0 : delete fgHistArray;
162 0 : fgHistArray = 0x0;
163 0 : }
164 :
165 : THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
166 : {
167 : //
168 : //
169 : //
170 0 : if(iter<fgHistArray->GetSize())
171 0 : return (THnBase*)fgHistArray->At(iter);
172 : else
173 0 : return 0x0;
174 0 : }
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 0 : delete fgHistArray;
183 0 : fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
184 0 : list->Add(fgHistArray);
185 0 : }
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 0 : TFile fcalib(filename);
196 0 : TObjArray * array = (TObjArray*)fcalib.Get(listname);
197 0 : const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
198 :
199 : AliTRDdEdxCalibHistArray * tmph=0x0;
200 0 : if(array){
201 0 : tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
202 0 : }
203 : else{
204 0 : tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
205 : }
206 0 : if(!tmph){
207 0 : printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
208 0 : fcalib.ls();
209 0 : if(array){
210 0 : array->ls();
211 0 : delete array;
212 : }
213 0 : return kFALSE;
214 : }
215 :
216 0 : delete fgHistArray;
217 0 : fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
218 :
219 0 : return kTRUE;
220 0 : }
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 0 : if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
228 :
229 0 : for(Int_t ii=0; ii<ncls; ii++){
230 0 : Double_t dq = (*arrayQ)[ii]/scale;
231 0 : const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
232 0 : const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
233 :
234 0 : if(dq<qmin)
235 0 : dq = qmin;
236 0 : if(dq>qmax)
237 0 : dq = qmax;
238 :
239 0 : const Double_t xx = (*arrayX)[ii];
240 0 : const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
241 0 : const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
242 :
243 0 : if(xx>=xmax || xx<xmin){
244 0 : printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(), xx, xmin, xmax); exit(1);
245 : }
246 :
247 0 : const Double_t var[]={xx, dq};
248 0 : hcalib->Fill(var);
249 0 : }
250 0 : }
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 0 : if(!fgHistArray){
258 0 : printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
259 : }
260 :
261 0 : THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
262 :
263 0 : TVectorD arrayQ(200), arrayX(200);
264 0 : const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
265 0 : FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
266 :
267 : static Int_t kprint = 100;
268 0 : if(kprint<0){
269 0 : printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
270 0 : printf("\nkinvq= %d;\n", kinvq);
271 0 : for(Int_t iq=0; iq<ncls; iq++){
272 0 : printf("arrayX[%3d] = %15.0f; arrayQ[%3d] = %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
273 : }
274 0 : printf("\n");
275 : }
276 0 : kprint++;
277 0 : }
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 0 : if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
290 0 : return -999;
291 :
292 0 : if(tpcsig<EPSILON)
293 0 : return -999;
294 :
295 0 : return tpcsig/45;
296 :
297 0 : }
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 0 : TObjArray *obj=new TObjArray(ndet);
306 0 : obj->SetOwner();
307 0 : for(Int_t ii=0; ii<ndet; ii++){
308 0 : obj->Add(new TVectorD(AliTRDseedV1::kNtb));
309 : }
310 :
311 : //ibin = binlowedge of bin(ibin+1) = the number fills this bin
312 0 : for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
313 0 : const Double_t stat = hnor->GetBinContent(ibin+1);
314 0 : if(stat<EPSILON){
315 0 : continue;
316 : }
317 :
318 0 : const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
319 0 : const Int_t itb = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
320 0 : TVectorD *vec=(TVectorD *)obj->At(idet);
321 0 : (*vec)[itb] = stat;
322 0 : }
323 :
324 0 : hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
325 0 : for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
326 0 : const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
327 0 : const TVectorD *vec=(TVectorD *)obj->At(idet);
328 :
329 : Int_t nonzero = 0;
330 0 : for(Int_t ii=0; ii<vec->GetNrows(); ii++){
331 0 : if((*vec)[ii]>EPSILON){
332 0 : nonzero++;
333 0 : }
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 0 : 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 0 : mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
342 0 : }
343 :
344 0 : hmean->SetBinContent(ibin+1, mean);
345 : }
346 :
347 0 : delete obj;
348 0 : }
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 0 : const TString hname = hh->GetName();
367 0 : const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
368 :
369 : //----------------------------------------
370 : // Define nbin, tag, cobj0
371 : //----------------------------------------
372 0 : const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
373 :
374 :
375 0 : TString tag(hname);
376 0 : tag.ReplaceAll("Hist","Obj");
377 :
378 0 : TObjArray * cobj0 = new TObjArray(1);
379 0 : cobj0->SetOwner();
380 0 : cobj0->SetName(tag);
381 0 : cobj0->Add(new TVectorD(nbin));
382 :
383 : //----------------------------------------
384 : // Define lowFrac, highFrac
385 : //----------------------------------------
386 : Double_t lowFrac = -999, highFrac = -999;
387 0 : if(!kinvq) {
388 0 : lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
389 0 : }
390 : else{
391 0 : lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
392 : }
393 :
394 : //----------------------------------------
395 : // Get analysis result
396 : //----------------------------------------
397 0 : TH1::AddDirectory(kFALSE);//important!
398 0 : TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
399 0 : TH2D *hpj = hh->Projection(1,0);
400 0 : AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
401 0 : GetPHCountMeanRMS(hnor, htrdphmean);
402 0 : if(lout){
403 0 : lout->Add(htrdphmean);
404 : }
405 0 : delete hpj;
406 :
407 0 : if(lout){
408 0 : lout->Add(hnor);
409 0 : lout->Add(hmpv);
410 0 : lout->Add(hwid);
411 0 : lout->Add(hres);
412 : }
413 :
414 : //----------------------------------------
415 : // Define Counter
416 : //----------------------------------------
417 : TVectorD *countDet=0x0;
418 : TObjArray *countSSL=0x0;
419 :
420 0 : if(!kinvq){
421 0 : countDet=new TVectorD(540);
422 0 : countSSL=new TObjArray(90);//SectorStackLayer
423 0 : countSSL->SetOwner();
424 0 : for(Int_t ii=0; ii<90; ii++){
425 0 : countSSL->Add(new TVectorD(6));
426 : }
427 0 : }
428 :
429 : //----------------------------------------
430 : // Fill cobj0
431 : //----------------------------------------
432 :
433 : //ibin = binlowedge of bin(ibin+1) = the number fills this bin
434 0 : for(Int_t ibin=0; ibin<nbin; ibin++){
435 0 : Double_t gnor = hnor->GetBinContent(ibin+1);
436 0 : Double_t gmpv = hmpv->GetBinContent(ibin+1);
437 0 : Double_t gwid = hwid->GetBinContent(ibin+1);
438 0 : Double_t gres = hres->GetBinContent(ibin+1);
439 :
440 : //--- set additional cut by kpass
441 : Bool_t kpass = kTRUE;
442 0 : Double_t gtrdphmean = -999;
443 0 : if(htrdphmean){
444 0 : gtrdphmean = htrdphmean->GetBinContent(ibin+1);
445 : //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
446 0 : if(gtrdphmean<EPSILON){
447 : kpass = kFALSE;
448 0 : }
449 0 : if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
450 : kpass = kFALSE;
451 0 : }
452 : }
453 :
454 : //--- set calibration constant p0
455 0 : 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 0 : if(gmpv>EPSILON && kpass){
462 0 : if(tag.Contains("T0")){
463 0 : p0 = gmpv;
464 0 : }
465 : else{
466 0 : p0 = 1/gmpv;
467 : }
468 : //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
469 : }
470 :
471 0 : (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
472 :
473 : //--- save optional record
474 0 : if(p0>EPSILON && countDet && countSSL){
475 0 : const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
476 0 : (*countDet)[idet]=1;
477 :
478 0 : const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
479 0 : const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
480 0 : const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
481 0 : TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
482 0 : (*vecsectorstack)[ilayer]=1;
483 0 : }
484 :
485 0 : if(calibStream){
486 0 : (*calibStream)<<tag<<
487 0 : "run="<<run<<
488 0 : "p0="<<p0<<
489 :
490 0 : "nor="<<gnor<<
491 0 : "mpv="<<gmpv<<
492 0 : "wid="<<gwid<<
493 0 : "res="<<gres<<
494 0 : "gtrdphmean="<<gtrdphmean<<
495 :
496 0 : "ibin="<<ibin<<
497 : "\n";
498 : }
499 0 : }
500 :
501 : //----------------------------------------
502 : // Status Report
503 : //----------------------------------------
504 0 : if(countDet && countSSL){
505 0 : TVectorD count2Dstack(90);
506 0 : for(Int_t ii=0; ii<90; ii++){
507 0 : TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
508 0 : const Int_t nlayer = (Int_t)vecsectorstack->Sum();
509 0 : if(nlayer==6){
510 0 : count2Dstack[ii]=1;
511 0 : }
512 : }
513 :
514 0 : 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 0 : }
516 :
517 : //----------------------------------------
518 : // Clean Up
519 : //----------------------------------------
520 :
521 0 : TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
522 : const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
523 0 : for(Int_t ihh=0; ihh<nhh; ihh++){
524 0 : if(!lout){
525 0 : delete (*hhs[ihh]);
526 : }
527 : }
528 :
529 0 : delete countDet;
530 0 : delete countSSL;
531 :
532 : //----------------------------------------
533 :
534 : return cobj0;
535 0 : }
536 :
|
