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