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