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 : //Root includes
18 : #include <TH1F.h>
19 : #include <TH1D.h>
20 : #include <TH2F.h>
21 : #include <TH3F.h>
22 : #include <TString.h>
23 : #include <TMath.h>
24 : #include <TF1.h>
25 : #include <TRandom.h>
26 : #include <TDirectory.h>
27 : #include <TFile.h>
28 : #include <TAxis.h>
29 : //AliRoot includes
30 : #include "AliRawReader.h"
31 : #include "AliRawReaderRoot.h"
32 : #include "AliRawReaderDate.h"
33 : #include "AliTPCCalROC.h"
34 : #include "AliTPCCalPad.h"
35 : #include "AliTPCROC.h"
36 : #include "AliMathBase.h"
37 : #include "TTreeStream.h"
38 :
39 : //date
40 : #include "event.h"
41 :
42 : //header file
43 : #include "AliTPCCalibKr.h"
44 :
45 : //----------------------------------------------------------------------------
46 : // The AliTPCCalibKr class description (TPC Kr calibration).
47 : //
48 : //
49 : // The AliTPCCalibKr keeps the array of TH3F histograms (TPC_max_padraw,TPC_max_pad,TPC_ADC_cluster),
50 : // its data memebers and is filled by AliTPCCalibKrTask under conditions (Accept()).
51 : //
52 : // The ouput TH3F histograms are later used to determine the calibration parameters of TPC chambers.
53 : // These calculations are done by using AliTPCCalibKr::Analyse() function. The ouput calibration
54 : // parameters (details in AliTPCCalibKr::Analyse()) are stored in the calibKr.root file for each TPC pad.
55 : // In addition the debugCalibKr.root file with debug information is created.
56 : //
57 :
58 : /*
59 :
60 : // Usage example:
61 : //
62 :
63 : // 1. Analyse output histograms:
64 : TFile f("outHistFile.root");
65 : AliTPCCalibKr *obj = (AliTPCCalibKr*) cOutput.FindObject("AliTPCCalibKr");
66 : obj->SetRadius(0,0);
67 : obj->SetStep(1,1);
68 : obj->Analyse();
69 :
70 : // 2. See calibration parameters e.g.:
71 : TFile f("calibKr.root");
72 : spectrMean->GetCalROC(70)->GetValue(40,40);
73 : fitMean->GetCalROC(70)->GetValue(40,40);
74 :
75 : // 3. See debug information e.g.:
76 : TFile f("debugCalibKr.root");
77 : .ls;
78 :
79 : // -- Print calibKr TTree content
80 : calibKr->Print();
81 :
82 : // -- Draw calibKr TTree variables
83 : calibKr.Draw("fitMean");
84 :
85 : */
86 :
87 :
88 : //
89 : // Author: Jacek Otwinowski (J.Otwinowski@gsi.de) and Stafan Geartner (S.Gaertner@gsi.de)
90 : //-----------------------------------------------------------------------------
91 :
92 6 : ClassImp(AliTPCCalibKr)
93 :
94 : AliTPCCalibKr::AliTPCCalibKr() :
95 0 : TObject(),
96 0 : fASide(kTRUE),
97 0 : fCSide(kTRUE),
98 0 : fHistoKrArray(72),
99 0 : fADCOverClustSizeMin(0.0),
100 0 : fADCOverClustSizeMax(1.0e9),
101 0 : fMaxADCOverClustADCMin(0.0),
102 0 : fMaxADCOverClustADCMax(1.0e9),
103 0 : fTimeMin(0.0),
104 0 : fTimeMax(1.0e9),
105 0 : fClustSizeMin(0.0),
106 0 : fClustSizeMax(1.0e9),
107 0 : fTimebinRmsIrocMin(0.0),
108 0 : fPadRmsIrocMin(0.0),
109 0 : fRowRmsIrocMin(0.0),
110 0 : fClusterPadSize1DIrocMax(200),
111 0 : fCurveCoefficientIroc(1.0e9),
112 0 : fTimebinRmsOrocMin(0.0),
113 0 : fPadRmsOrocMin(0.0),
114 0 : fRowRmsOrocMin(0.0),
115 0 : fClusterPadSize1DOrocMax(200),
116 0 : fCurveCoefficientOroc(1.0e9),
117 0 : fIrocHistogramMin(100.),
118 0 : fIrocHistogramMax(6000.),
119 0 : fIrocHistogramNbins(200),
120 0 : fOrocHistogramMin(100.),
121 0 : fOrocHistogramMax(5500.),
122 0 : fOrocHistogramNbins(200),
123 0 : fRowRadius(0),
124 0 : fPadRadius(0),
125 0 : fRowStep(1),
126 0 : fPadStep(1)
127 :
128 0 : {
129 : //
130 : // default constructor
131 : //
132 :
133 : // TObjArray with histograms
134 0 : fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
135 0 : fHistoKrArray.Clear();
136 :
137 : // init cuts (by Stefan)
138 : // SetADCOverClustSizeRange(7,200);
139 : // SetMaxADCOverClustADCRange(0.01,0.4);
140 : // SetTimeRange(200,600);
141 : // SetClustSizeRange(6,200);
142 :
143 : //init cuts (by Adam)
144 0 : SetTimebinRmsMin(1.6,0.8);
145 0 : SetPadRmsMin(0.825,0.55);
146 0 : SetRowRmsMin(0.1,0.1);
147 0 : SetClusterPadSize1DMax(15,11);
148 0 : SetCurveCoefficient(1.,2.);
149 :
150 : //set histograms settings
151 0 : SetIrocHistogram(200,100,6000);
152 0 : SetOrocHistogram(200,100,5500);
153 :
154 : // init histograms
155 : //Init();
156 0 : }
157 :
158 : //_____________________________________________________________________
159 : AliTPCCalibKr::AliTPCCalibKr(const AliTPCCalibKr& pad) :
160 0 : TObject(pad),
161 :
162 0 : fASide(pad.fASide),
163 0 : fCSide(pad.fCSide),
164 0 : fHistoKrArray(72),
165 0 : fADCOverClustSizeMin(pad.fADCOverClustSizeMin),
166 0 : fADCOverClustSizeMax(pad.fADCOverClustSizeMax),
167 0 : fMaxADCOverClustADCMin(pad.fMaxADCOverClustADCMin),
168 0 : fMaxADCOverClustADCMax(pad.fMaxADCOverClustADCMax),
169 0 : fTimeMin(pad.fTimeMin),
170 0 : fTimeMax(pad.fTimeMax),
171 0 : fClustSizeMin(pad.fClustSizeMin),
172 0 : fClustSizeMax(pad.fClustSizeMax),
173 0 : fTimebinRmsIrocMin(pad.fTimebinRmsIrocMin),
174 0 : fPadRmsIrocMin(pad.fPadRmsIrocMin),
175 0 : fRowRmsIrocMin(pad.fRowRmsIrocMin),
176 0 : fClusterPadSize1DIrocMax(pad.fClusterPadSize1DIrocMax),
177 0 : fCurveCoefficientIroc(pad.fCurveCoefficientIroc),
178 0 : fTimebinRmsOrocMin(pad.fTimebinRmsOrocMin),
179 0 : fPadRmsOrocMin(pad.fPadRmsOrocMin),
180 0 : fRowRmsOrocMin(pad.fRowRmsOrocMin),
181 0 : fClusterPadSize1DOrocMax(pad.fClusterPadSize1DOrocMax),
182 0 : fCurveCoefficientOroc(pad.fCurveCoefficientOroc),
183 0 : fIrocHistogramMin(pad.fIrocHistogramMin),
184 0 : fIrocHistogramMax(pad.fIrocHistogramMax),
185 0 : fIrocHistogramNbins(pad.fIrocHistogramNbins),
186 0 : fOrocHistogramMin(pad.fOrocHistogramMin),
187 0 : fOrocHistogramMax(pad.fOrocHistogramMax),
188 0 : fOrocHistogramNbins(pad.fOrocHistogramNbins),
189 0 : fRowRadius(pad.fRowRadius),
190 0 : fPadRadius(pad.fPadRadius),
191 0 : fRowStep(pad.fRowStep),
192 0 : fPadStep(pad.fPadStep)
193 :
194 0 : {
195 : // copy constructor
196 :
197 : // TObjArray with histograms
198 0 : fHistoKrArray.SetOwner(kTRUE); // is owner of histograms
199 0 : fHistoKrArray.Clear();
200 :
201 0 : for (Int_t iSec = 0; iSec < 72; ++iSec)
202 : {
203 0 : TH3F *hOld = pad.GetHistoKr(iSec);
204 0 : if(hOld) {
205 0 : TH3F *hNew = new TH3F( *pad.GetHistoKr(iSec) );
206 0 : fHistoKrArray.AddAt(hNew,iSec);
207 0 : }
208 : }
209 0 : }
210 :
211 : //_____________________________________________________________________
212 : AliTPCCalibKr::~AliTPCCalibKr()
213 0 : {
214 : //
215 : // destructor
216 : //
217 :
218 : // for (Int_t iSec = 0; iSec < 72; ++iSec) {
219 : // if (fHistoKrArray.At(iSec)) delete fHistoKrArray.RemoveAt(iSec);
220 : // }
221 0 : fHistoKrArray.Delete();
222 0 : }
223 :
224 : //_____________________________________________________________________
225 : AliTPCCalibKr& AliTPCCalibKr::operator = (const AliTPCCalibKr &source)
226 : {
227 : // assignment operator
228 :
229 0 : if (&source == this) return *this;
230 0 : new (this) AliTPCCalibKr(source);
231 :
232 0 : return *this;
233 0 : }
234 :
235 : //_____________________________________________________________________
236 : void AliTPCCalibKr::Init()
237 : {
238 : //
239 : // init output histograms
240 : //
241 :
242 : // add histograms to the TObjArray
243 0 : for(Int_t i=0; i<72; ++i) {
244 :
245 : // C - side
246 0 : if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
247 0 : TH3F *hist = CreateHisto(i);
248 0 : if(hist) fHistoKrArray.AddAt(hist,i);
249 0 : }
250 :
251 : // A - side
252 0 : if(IsCSide(i) == kFALSE && fASide == kTRUE) {
253 0 : TH3F *hist = CreateHisto(i);
254 0 : if(hist) fHistoKrArray.AddAt(hist,i);
255 0 : }
256 : }
257 0 : }
258 :
259 : //_____________________________________________________________________
260 : Bool_t AliTPCCalibKr::Process(AliTPCclusterKr *cluster)
261 : {
262 : //
263 : // process events
264 : // call event by event
265 : //
266 :
267 0 : if(cluster) Update(cluster);
268 0 : else return kFALSE;
269 :
270 0 : return kTRUE;
271 0 : }
272 :
273 : //_____________________________________________________________________
274 : TH3F* AliTPCCalibKr::CreateHisto(Int_t chamber)
275 : {
276 : //
277 : // create new histogram
278 : //
279 0 : char name[256];
280 : TH3F *h;
281 :
282 0 : snprintf(name,256,"ADCcluster_ch%d",chamber);
283 :
284 0 : if( IsIROC(chamber) == kTRUE )
285 : {
286 0 : h = new TH3F(name,name,63,0,63,108,0,108,fIrocHistogramNbins,fIrocHistogramMin,fIrocHistogramMax);
287 0 : } else {
288 0 : h = new TH3F(name,name,96,0,96,140,0,140,fOrocHistogramNbins,fOrocHistogramMin,fOrocHistogramMax);
289 : }
290 0 : h->SetXTitle("padrow");
291 0 : h->SetYTitle("pad");
292 0 : h->SetZTitle("fADC");
293 :
294 0 : return h;
295 0 : }
296 :
297 : //_____________________________________________________________________
298 : Bool_t AliTPCCalibKr::IsIROC(Int_t chamber)
299 : {
300 : // check if IROCs
301 : // returns kTRUE if IROCs and kFALSE if OROCs
302 :
303 0 : if(chamber>=0 && chamber<36) return kTRUE;
304 :
305 0 : return kFALSE;
306 0 : }
307 :
308 : //_____________________________________________________________________
309 : Bool_t AliTPCCalibKr::IsCSide(Int_t chamber)
310 : {
311 : // check if C side
312 : // returns kTRUE if C side and kFALSE if A side
313 :
314 0 : if((chamber>=18 && chamber<36) || (chamber>=54 && chamber<72)) return kTRUE;
315 :
316 0 : return kFALSE;
317 0 : }
318 :
319 : //_____________________________________________________________________
320 : Bool_t AliTPCCalibKr::Update(AliTPCclusterKr *cl)
321 : {
322 : //
323 : // fill existing histograms
324 : //
325 :
326 0 : if (!Accept(cl)) return kFALSE;
327 0 : TH3F *h = (TH3F*)fHistoKrArray.At(cl->GetSec());
328 0 : if(!h) return kFALSE;
329 :
330 0 : h->Fill(cl->GetMax().GetRow(),cl->GetMax().GetPad(),cl->GetADCcluster());
331 :
332 0 : return kTRUE;
333 0 : }
334 :
335 : //_____________________________________________________________________
336 : Bool_t AliTPCCalibKr::Accept(AliTPCclusterKr *cl){
337 : //
338 : // cuts
339 : //
340 : /*
341 : TCut cutR0("cutR0","fADCcluster/fSize<200"); // adjust it according v seetings -
342 : TCut cutR1("cutR1","fADCcluster/fSize>7"); // cosmic tracks and noise removal
343 : TCut cutR2("cutR2","fMax.fAdc/fADCcluster<0.4"); // digital noise removal
344 : TCut cutR3("cutR3","fMax.fAdc/fADCcluster>0.01"); // noise removal
345 : TCut cutR4("cutR4","fMax.fTime>200"); // noise removal
346 : TCut cutR5("cutR5","fMax.fTime<600"); // noise removal
347 : TCut cutS1("cutS1","fSize<200"); // adjust it according v seetings - cosmic tracks
348 : TCut cutAll = cutR0+cutR1+cutR2+cutR3+cutR4+cutR5+cutS1;
349 : */
350 : /*
351 : //R0
352 : if ((float)cl->GetADCcluster()/ cl->GetSize() >200) return kFALSE;
353 : // R1
354 : if ((float)cl->GetADCcluster()/ cl->GetSize() <7) return kFALSE;
355 : //R2
356 : if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >0.4) return kFALSE;
357 : //R3
358 : if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <0.01) return kFALSE;
359 : //R4
360 : if (cl->GetMax().GetTime() < 200) return kFALSE;
361 : //R5
362 : if (cl->GetMax().GetTime() > 600) return kFALSE;
363 : //S1
364 : if (cl->GetSize()>200) return kFALSE;
365 : if (cl->GetSize()<6) return kFALSE;
366 :
367 : SetADCOverClustSizeRange(7,200);
368 : SetMaxADCOverClustADCRange(0.01,0.4);
369 : SetTimeRange(200,600);
370 : SetClustSizeRange(6,200);
371 : */
372 :
373 : //R0
374 0 : if ((float)cl->GetADCcluster()/ cl->GetSize() >fADCOverClustSizeMax) return kFALSE;
375 : // R1
376 0 : if ((float)cl->GetADCcluster()/ cl->GetSize() <fADCOverClustSizeMin) return kFALSE;
377 : //R2
378 0 : if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() >fMaxADCOverClustADCMax) return kFALSE;
379 : //R3
380 0 : if ((float)cl->GetMax().GetAdc()/ cl->GetADCcluster() <fMaxADCOverClustADCMin) return kFALSE;
381 : //R4
382 0 : if (cl->GetMax().GetTime() > fTimeMax) return kFALSE;
383 : //R5
384 0 : if (cl->GetMax().GetTime() < fTimeMin) return kFALSE;
385 : //S1
386 0 : if (cl->GetSize()>fClustSizeMax) return kFALSE;
387 0 : if (cl->GetSize()<fClustSizeMin) return kFALSE;
388 :
389 : //
390 : // cuts by Adam
391 : //
392 : /*
393 : TCut cutAI0("cutAI0","fTimebinRMS>1.6");
394 : TCut cutAI1("cutAI1","fPadRMS>0.825");
395 : TCut cutAI2("cutAI2","fRowRMS>0.1");
396 : TCut cutAI3("cutAI3","fPads1D<15");
397 : TCut cutAI4("cutAI4","fTimebinRMS+11.9-2.15*TMath::Log(1.*fADCcluster)<0");
398 :
399 : TCut cutAIAll = cutAI0+cutAI1+cutAI2+cutAI3+cutAI4;
400 :
401 : TCut cutAO0("cutAO0","fTimebinRMS>0.8");
402 : TCut cutAO1("cutAO1","fPadRMS>0.55");
403 : TCut cutAO2("cutAO2","fRowRMS>0.1");
404 : TCut cutAO3("cutAO3","fPads1D<11");
405 : TCut cutAO4("cutAO4","fTimebinRMS+11.9-2.15*TMath::Log(2.*fADCcluster)<0");
406 :
407 : TCut cutAOAll = cutAO0+cutAO1+cutAO2+cutAO3+cutAO4;
408 : TCut cutAll("cutAll","(fSec<36&&cutAIAll)||(fSec>=36&&cutAOAll)");
409 :
410 : */
411 0 : if(cl->GetSec()<36){ //IROCs
412 : //AI0
413 0 : if((float)cl->GetTimebinRMS() <= fTimebinRmsIrocMin) return kFALSE;
414 : //AI1
415 0 : if((float)cl->GetPadRMS() <= fPadRmsIrocMin) return kFALSE;
416 : //AI2
417 0 : if((float)cl->GetRowRMS() <= fRowRmsIrocMin) return kFALSE;
418 : //AI3
419 0 : if(cl->GetPads1D() >= fClusterPadSize1DIrocMax) return kFALSE;
420 : //AI4
421 0 : if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientIroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
422 : }else{ //OROCs
423 : //AO0
424 0 : if((float)cl->GetTimebinRMS() <= fTimebinRmsOrocMin) return kFALSE;
425 : //AO1
426 0 : if((float)cl->GetPadRMS() <= fPadRmsOrocMin) return kFALSE;
427 : //AO2
428 0 : if((float)cl->GetRowRMS() <= fRowRmsOrocMin) return kFALSE;
429 : //AO3
430 0 : if(cl->GetPads1D() >= fClusterPadSize1DOrocMax) return kFALSE;
431 : //AO4
432 0 : if((float)cl->GetTimebinRMS()+11.9-2.15*TMath::Log(fCurveCoefficientOroc*(float)cl->GetADCcluster()) >= 0) return kFALSE;
433 : }
434 :
435 0 : return kTRUE;
436 :
437 0 : }
438 :
439 : //_____________________________________________________________________
440 : TH3F* AliTPCCalibKr::GetHistoKr(Int_t chamber) const
441 : {
442 : // get histograms from fHistoKrArray
443 0 : return (TH3F*) fHistoKrArray.At(chamber);
444 : }
445 :
446 : //_____________________________________________________________________
447 : void AliTPCCalibKr::Analyse()
448 : {
449 : //
450 : // analyse the histograms and extract krypton calibration parameters
451 : //
452 :
453 : // AliTPCCalPads that will contain the calibration parameters
454 0 : AliTPCCalPad* spectrMeanCalPad = new AliTPCCalPad("spectrMean", "spectrMean");
455 0 : AliTPCCalPad* spectrRMSCalPad = new AliTPCCalPad("spectrRMS", "spectrRMS");
456 0 : AliTPCCalPad* fitMeanCalPad = new AliTPCCalPad("fitMean", "fitMean");
457 0 : AliTPCCalPad* fitRMSCalPad = new AliTPCCalPad("fitRMS", "fitRMS");
458 0 : AliTPCCalPad* fitNormChi2CalPad = new AliTPCCalPad("fitNormChi2", "fitNormChi2");
459 0 : AliTPCCalPad* entriesCalPad = new AliTPCCalPad("entries", "entries");
460 :
461 : // file stream for debugging purposes
462 0 : TTreeSRedirector* debugStream = new TTreeSRedirector("debugCalibKr.root");
463 :
464 : // if entries in spectrum less than minEntries, then the fit won't be performed
465 : Int_t minEntries = 1; //300;
466 :
467 : Double_t windowFrac = 0.12;
468 : // the 3d histogram will be projected on the pads given by the following window size
469 : // set the numbers to 0 if you want to do a pad-by-pad calibration
470 0 : UInt_t rowRadius = fRowRadius;//4;
471 0 : UInt_t padRadius = fPadRadius;//4;
472 :
473 : // the step size by which pad and row are incremented is given by the following numbers
474 : // set them to 1 if you want the finest granularity
475 0 : UInt_t rowStep = fRowStep;//1;//2 // formerly: 2*rowRadius
476 0 : UInt_t padStep = fPadStep;//1;//4 // formerly: 2*padRadius
477 :
478 0 : for (Int_t chamber = 0; chamber < 72; chamber++) {
479 : //if (chamber != 71) continue;
480 0 : AliTPCCalROC roc(chamber); // I need this only for GetNrows() and GetNPads()
481 :
482 : // Usually I would traverse each pad, take the spectrum for its neighbourhood and
483 : // obtain the calibration parameters. This takes very long, so instead I assign the same
484 : // calibration values to the whole neighbourhood and then go on to the next neighbourhood.
485 0 : UInt_t nRows = roc.GetNrows();
486 0 : for (UInt_t iRow = 0; iRow < nRows; iRow += rowStep) {
487 0 : UInt_t nPads = roc.GetNPads(iRow);
488 : //if (iRow >= 10) break;
489 0 : for (UInt_t iPad = 0; iPad < nPads; iPad += padStep) {
490 : //if (iPad >= 20) break;
491 0 : TH3F* h = GetHistoKr(chamber);
492 0 : if (!h) continue;
493 :
494 : // the 3d histogram will be projected on the pads given by the following bounds
495 : // for rows and pads
496 0 : Int_t rowLow = iRow - rowRadius;
497 0 : UInt_t rowUp = iRow + rowRadius + rowStep-1;
498 0 : Int_t padLow = iPad - padRadius;
499 0 : UInt_t padUp = iPad + padRadius + padStep-1;
500 : // if window goes out of chamber
501 0 : if (rowLow < 0) rowLow = 0;
502 0 : if (rowUp >= nRows) rowUp = nRows - 1;
503 0 : if (padLow < 0) padLow = 0;
504 0 : if (padUp >= nPads) padUp = nPads - 1;
505 :
506 : // project the histogram
507 : //TH1D* projH = h->ProjectionZ("projH", rowLow+1, rowUp+1, padLow+1, padUp+1); // SLOW
508 0 : TH1D* projH = ProjectHisto(h, "projH", rowLow+1, rowUp+1, padLow+1, padUp+1);
509 :
510 : // get the number of entries in the spectrum
511 0 : Double_t entries = projH->GetEntries();
512 0 : if (entries < minEntries) { delete projH; continue; }
513 :
514 : // get the two calibration parameters mean of spectrum and RMS of spectrum
515 0 : Double_t histMean = projH->GetMean();
516 0 : Double_t histRMS = (histMean != 0) ? projH->GetRMS() / histMean : 0.;
517 :
518 : // find maximum in spectrum to define a range (given by windowFrac) for which a Gauss is fitted
519 0 : Double_t maxEntries = projH->GetBinCenter(projH->GetMaximumBin());
520 0 : Int_t minBin = projH->FindBin((1.-windowFrac) * maxEntries);
521 0 : Int_t maxBin = projH->FindBin((1.+windowFrac) * maxEntries);
522 0 : Double_t integCharge = projH->Integral(minBin,maxBin);
523 :
524 0 : Int_t fitResult = projH->Fit("gaus", "Q0", "", (1.-windowFrac) * maxEntries, (1.+windowFrac) * maxEntries);
525 :
526 0 : if (fitResult != 0) {
527 0 : Error("Analyse", "Error while fitting spectrum for chamber %i, rows %i - %i, pads %i - %i, integrated charge %f.", chamber, rowLow, rowUp, padLow, padUp, integCharge);
528 : //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, %f\n", chamber, iRow, iPad, entries, maxEntries);
529 :
530 0 : delete projH;
531 0 : continue;
532 : }
533 :
534 : // get the two calibration parameters mean of gauss fit and sigma of gauss fit
535 0 : TF1* gausFit = projH->GetFunction("gaus");
536 0 : Double_t fitMean = gausFit->GetParameter(1);
537 0 : Double_t fitRMS = gausFit->GetParameter(2);
538 0 : Int_t numberFitPoints = gausFit->GetNumberFitPoints();
539 0 : if (numberFitPoints == 0) continue;
540 0 : Double_t fitNormChi2 = gausFit->GetChisquare() / numberFitPoints;
541 0 : delete gausFit;
542 0 : delete projH;
543 0 : if (fitMean <= 0) continue;
544 : //printf("[ch%i r%i, p%i] entries = %f, maxEntries = %f, fitMean = %f, fitRMS = %f\n", chamber, iRow, iPad, entries, maxEntries, fitMean, fitRMS);
545 :
546 : // write the calibration parameters for each pad that the 3d histogram was projected onto
547 : // (with considering the step size) to the CalPads
548 : // rowStep (padStep) odd: round down s/2 and fill this # of rows (pads) in both directions
549 : // rowStep (padStep) even: fill s/2 rows (pads) in ascending direction, s/2-1 in descending direction
550 0 : for (Int_t r = iRow - (rowStep/2 - (rowStep+1)%2); r <= (Int_t)(iRow + rowStep/2); r++) {
551 0 : if (r < 0 || r >= (Int_t)nRows) continue;
552 0 : UInt_t nPadsR = roc.GetNPads(r);
553 0 : for (Int_t p = iPad - (padStep/2 - (padStep+1)%2); p <= (Int_t)(iPad + padStep/2); p++) {
554 0 : if (p < 0 || p >= (Int_t)nPadsR) continue;
555 0 : spectrMeanCalPad->GetCalROC(chamber)->SetValue(r, p, histMean);
556 0 : spectrRMSCalPad->GetCalROC(chamber)->SetValue(r, p, histRMS);
557 0 : fitMeanCalPad->GetCalROC(chamber)->SetValue(r, p, fitMean);
558 0 : fitRMSCalPad->GetCalROC(chamber)->SetValue(r, p, fitRMS);
559 0 : fitNormChi2CalPad->GetCalROC(chamber)->SetValue(r, p, fitNormChi2);
560 0 : entriesCalPad->GetCalROC(chamber)->SetValue(r, p, entries);
561 :
562 0 : (*debugStream) << "calibKr" <<
563 0 : "sector=" << chamber << // chamber number
564 0 : "row=" << r << // row number
565 0 : "pad=" << p << // pad number
566 0 : "histMean=" << histMean << // mean of the spectrum
567 0 : "histRMS=" << histRMS << // RMS of the spectrum divided by the mean
568 0 : "fitMean=" << fitMean << // Gauss fitted mean of the 41.6 keV Kr peak
569 0 : "fitRMS=" << fitRMS << // Gauss fitted sigma of the 41.6 keV Kr peak
570 0 : "fitNormChi2" << fitNormChi2 << // normalized chi square of the Gauss fit
571 0 : "entries=" << entries << // number of entries for the spectrum
572 : "\n";
573 : }
574 0 : }
575 0 : }
576 : }
577 0 : }
578 :
579 0 : TFile f("calibKr.root", "recreate");
580 0 : spectrMeanCalPad->Write();
581 0 : spectrRMSCalPad->Write();
582 0 : fitMeanCalPad->Write();
583 0 : fitRMSCalPad->Write();
584 0 : fitNormChi2CalPad->Write();
585 0 : entriesCalPad->Write();
586 0 : f.Close();
587 0 : delete spectrMeanCalPad;
588 0 : delete spectrRMSCalPad;
589 0 : delete fitMeanCalPad;
590 0 : delete fitRMSCalPad;
591 0 : delete fitNormChi2CalPad;
592 0 : delete entriesCalPad;
593 0 : delete debugStream;
594 0 : }
595 :
596 : //_____________________________________________________________________
597 : TH1D* AliTPCCalibKr::ProjectHisto(TH3F* histo3D, const char* name, Int_t xMin, Int_t xMax, Int_t yMin, Int_t yMax)
598 : {
599 : // project the z-axis of a 3d histo to a specific range of the x- and y-axes,
600 : // replaces TH3F::ProjectZ() to gain more speed
601 :
602 0 : TAxis* xAxis = histo3D->GetXaxis();
603 0 : TAxis* yAxis = histo3D->GetYaxis();
604 0 : TAxis* zAxis = histo3D->GetZaxis();
605 0 : Double_t zMinVal = zAxis->GetXmin();
606 0 : Double_t zMaxVal = zAxis->GetXmax();
607 :
608 0 : Int_t nBinsZ = zAxis->GetNbins();
609 0 : TH1D* projH = new TH1D(name, name, nBinsZ, zMinVal, zMaxVal);
610 :
611 0 : Int_t nx = xAxis->GetNbins()+2;
612 0 : Int_t ny = yAxis->GetNbins()+2;
613 : Int_t bin = 0;
614 : Double_t entries = 0.;
615 0 : for (Int_t x = xMin; x <= xMax; x++) {
616 0 : for (Int_t y = yMin; y <= yMax; y++) {
617 0 : for (Int_t z = 0; z <= nBinsZ+1; z++) {
618 0 : bin = x + nx * (y + ny * z);
619 0 : Double_t val = histo3D->GetBinContent(bin);
620 0 : projH->Fill(zAxis->GetBinCenter(z), val);
621 0 : entries += val;
622 : }
623 : }
624 : }
625 0 : projH->SetEntries((Long64_t)entries);
626 0 : return projH;
627 0 : }
628 :
629 : //_____________________________________________________________________
630 : Long64_t AliTPCCalibKr::Merge(TCollection* list) {
631 : // merge function
632 : //
633 :
634 0 : if (!list)
635 0 : return 0;
636 :
637 0 : if (list->IsEmpty())
638 0 : return 1;
639 :
640 0 : TIterator* iter = list->MakeIterator();
641 : TObject* obj = 0;
642 :
643 0 : iter->Reset();
644 : Int_t count=0;
645 0 : while((obj = iter->Next()) != 0)
646 : {
647 0 : AliTPCCalibKr* entry = dynamic_cast<AliTPCCalibKr*>(obj);
648 0 : if (entry == 0) continue;
649 :
650 0 : for(int i=0; i<72; ++i) {
651 0 : if(IsCSide(i) == kTRUE && fCSide == kTRUE) {
652 0 : ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
653 0 : }
654 :
655 0 : if(IsCSide(i) == kFALSE && fASide == kTRUE) {
656 0 : ((TH3F*)fHistoKrArray.At(i))->Add( ((TH3F*)entry->fHistoKrArray.At(i)) );
657 0 : }
658 : }
659 :
660 0 : count++;
661 0 : }
662 :
663 0 : return count;
664 0 : }
|
