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