]>
Commit | Line | Data |
---|---|---|
eb7e0771 | 1 | /************************************************************************** |
2 | * Copyright(c) 2007-08, 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 | /////////////////////////////////////////////////////////////////////////////// | |
18 | // // | |
19 | // Class for Evaluation and Validation of the ALTRO Tail Cancelation Filter // | |
20 | // (TCF) parameters out of TPC Raw data // | |
21 | // // | |
d0bd4fcc | 22 | // Author: Stefan Rossegger, Simon Feigl // |
eb7e0771 | 23 | // // |
24 | /////////////////////////////////////////////////////////////////////////////// | |
25 | ||
26 | #include "AliTPCCalibTCF.h" | |
27 | ||
28 | #include <TObject.h> | |
29 | #include <TCanvas.h> | |
30 | #include <TFile.h> | |
31 | #include <TKey.h> | |
32 | #include <TStyle.h> | |
33 | #include <TMinuit.h> | |
34 | #include <TH1F.h> | |
d0bd4fcc | 35 | #include <TH2F.h> |
7409c8db | 36 | #include <AliSysInfo.h> |
eb7e0771 | 37 | |
38 | #include <TMath.h> | |
39 | #include <TNtuple.h> | |
40 | #include <TEntryList.h> | |
eb7e0771 | 41 | #include "AliRawReaderRoot.h" |
55f06b51 | 42 | #include "AliRawHLTManager.h" |
eb7e0771 | 43 | #include "AliTPCRawStream.h" |
752b0cc7 | 44 | #include "AliTPCRawStreamV3.h" |
eb7e0771 | 45 | #include "AliTPCROC.h" |
46 | ||
47 | #include "AliTPCAltroEmulator.h" | |
48 | ||
df4cfd77 | 49 | #include "AliTPCmapper.h" |
50 | #include <fstream> | |
51 | ||
eb7e0771 | 52 | ClassImp(AliTPCCalibTCF) |
53 | ||
54 | AliTPCCalibTCF::AliTPCCalibTCF() : | |
55 | TNamed(), | |
752b0cc7 | 56 | fGateWidth(50), |
eb7e0771 | 57 | fSample(900), |
752b0cc7 | 58 | fPulseLength(400), |
eb7e0771 | 59 | fLowPulseLim(30), |
752b0cc7 | 60 | fUpPulseLim(900), |
61 | fRMSLim(1.0), | |
62 | fRatioIntLim(2) | |
d0bd4fcc | 63 | |
eb7e0771 | 64 | { |
65 | // | |
66 | // AliTPCCalibTCF standard constructor | |
67 | // | |
68 | } | |
69 | ||
70 | //_____________________________________________________________________________ | |
d0bd4fcc | 71 | AliTPCCalibTCF::AliTPCCalibTCF(Int_t gateWidth, Int_t sample, Int_t pulseLength, Int_t lowPulseLim, Int_t upPulseLim, Double_t rmsLim, Double_t ratioIntLim) : |
eb7e0771 | 72 | TNamed(), |
73 | fGateWidth(gateWidth), | |
74 | fSample(sample), | |
75 | fPulseLength(pulseLength), | |
76 | fLowPulseLim(lowPulseLim), | |
77 | fUpPulseLim(upPulseLim), | |
d0bd4fcc | 78 | fRMSLim(rmsLim), |
79 | fRatioIntLim(ratioIntLim) | |
eb7e0771 | 80 | { |
81 | // | |
82 | // AliTPCCalibTCF constructor with specific (non-standard) thresholds | |
83 | // | |
84 | } | |
85 | ||
86 | //_____________________________________________________________________________ | |
87 | AliTPCCalibTCF::AliTPCCalibTCF(const AliTPCCalibTCF &tcf) : | |
88 | TNamed(tcf), | |
89 | fGateWidth(tcf.fGateWidth), | |
90 | fSample(tcf.fSample), | |
91 | fPulseLength(tcf.fPulseLength), | |
92 | fLowPulseLim(tcf.fLowPulseLim), | |
93 | fUpPulseLim(tcf.fUpPulseLim), | |
d0bd4fcc | 94 | fRMSLim(tcf.fRMSLim), |
95 | fRatioIntLim(tcf.fRatioIntLim) | |
eb7e0771 | 96 | { |
97 | // | |
98 | // AliTPCCalibTCF copy constructor | |
99 | // | |
100 | } | |
101 | ||
102 | ||
103 | //_____________________________________________________________________________ | |
104 | AliTPCCalibTCF& AliTPCCalibTCF::operator = (const AliTPCCalibTCF &source) | |
105 | { | |
106 | // | |
107 | // AliTPCCalibTCF assignment operator | |
108 | // | |
109 | ||
110 | if (&source == this) return *this; | |
111 | new (this) AliTPCCalibTCF(source); | |
112 | ||
113 | return *this; | |
114 | ||
115 | } | |
116 | ||
117 | //_____________________________________________________________________________ | |
118 | AliTPCCalibTCF::~AliTPCCalibTCF() | |
119 | { | |
120 | // | |
121 | // AliTPCCalibTCF destructor | |
122 | // | |
123 | } | |
124 | ||
752b0cc7 | 125 | |
126 | //_____________________________________________________________________________ | |
127 | void AliTPCCalibTCF::ProcessRawFileV3(const char *nameRawFile, const char *nameFileOut) { | |
128 | // | |
129 | // New RCU data format!: Standard middle of 2009 | |
130 | // | |
131 | // Loops over all events within one RawData file and collects proper pulses | |
132 | // (according to given tresholds) per pad | |
133 | // Histograms per pad are stored in 'nameFileOut' | |
134 | // | |
135 | ||
136 | AliRawReader *rawReader = AliRawReader::Create(nameRawFile); | |
137 | if (!rawReader) { | |
138 | printf("Could not create a raw reader for %s\n",nameRawFile); | |
139 | return; | |
140 | } | |
141 | ||
142 | rawReader->RewindEvents(); // just to make sure | |
143 | ||
144 | rawReader->Select("TPC"); | |
145 | ||
146 | if (!rawReader->NextEvent()) { | |
147 | printf("no events found in %s\n",nameRawFile); | |
148 | return; | |
149 | } | |
150 | ||
151 | // TPC stream reader | |
152 | AliTPCRawStreamV3 rawStream(rawReader); | |
153 | ||
154 | Int_t ievent=0; | |
155 | do { | |
156 | AliSysInfo::AddStamp(Form("start_event_%d",ievent), ievent,-1,-1); | |
157 | printf("Reading next event ... Nr: %d\n",ievent); | |
158 | // Start the basic data extraction | |
159 | ProcessRawEventV3(rawReader, &rawStream, nameFileOut); | |
160 | AliSysInfo::AddStamp(Form("end_event_%d",ievent), ievent,-1,-1); | |
161 | ievent++; | |
162 | ||
163 | } while (rawReader->NextEvent()); | |
164 | ||
165 | rawReader->~AliRawReader(); | |
166 | ||
167 | } | |
168 | ||
169 | ||
eb7e0771 | 170 | //_____________________________________________________________________________ |
55f06b51 | 171 | void AliTPCCalibTCF::ProcessRawFile(const char *nameRawFile, const char *nameFileOut, bool bUseHLTOUT) { |
eb7e0771 | 172 | // |
173 | // Loops over all events within one RawData file and collects proper pulses | |
174 | // (according to given tresholds) per pad | |
175 | // Histograms per pad are stored in 'nameFileOut' | |
176 | // | |
177 | ||
55f06b51 | 178 | // create the data reader |
eb7e0771 | 179 | AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile); |
55f06b51 | 180 | if (!rawReader) { |
181 | return; | |
182 | } | |
183 | ||
184 | // create HLT reader for redirection of TPC data from HLTOUT to TPC reconstruction | |
185 | AliRawReader *hltReader=AliRawHLTManager::AliRawHLTManager::CreateRawReaderHLT(rawReader, "TPC"); | |
186 | ||
187 | // now choose the data source | |
188 | if (bUseHLTOUT) rawReader=hltReader; | |
189 | ||
190 | // rawReader->Reset(); | |
191 | rawReader->RewindEvents(); | |
192 | ||
193 | if (!rawReader->NextEvent()) { | |
194 | printf("no events found in %s\n",nameRawFile); | |
195 | return; | |
196 | } | |
eb7e0771 | 197 | |
d0bd4fcc | 198 | Int_t ievent=0; |
55f06b51 | 199 | do { |
7409c8db | 200 | AliSysInfo::AddStamp(Form("start_event_%d",ievent), ievent,-1,-1); |
d0bd4fcc | 201 | printf("Reading next event ... Nr: %d\n",ievent); |
7409c8db | 202 | AliTPCRawStream *rawStream = new AliTPCRawStream(rawReader); |
eb7e0771 | 203 | rawReader->Select("TPC"); |
7409c8db | 204 | ProcessRawEvent(rawStream, nameFileOut); |
205 | delete rawStream; | |
206 | AliSysInfo::AddStamp(Form("end_event_%d",ievent), ievent,-1,-1); | |
d0bd4fcc | 207 | ievent++; |
55f06b51 | 208 | } while (rawReader->NextEvent()); |
eb7e0771 | 209 | |
210 | rawReader->~AliRawReader(); | |
211 | ||
212 | } | |
213 | ||
214 | ||
752b0cc7 | 215 | //_____________________________________________________________________________ |
216 | void AliTPCCalibTCF::ProcessRawEventV3( AliRawReader *rawReader, AliTPCRawStreamV3 *rawStream, const char *nameFileOut) { | |
217 | // | |
218 | // New RCU data format!: Standard middle of 2009 | |
219 | // | |
220 | // Extracts proper pulses (according the given tresholds) within one event | |
221 | // and accumulates them into one histogram per pad. All histograms are | |
222 | // saved in the file 'nameFileOut'. | |
223 | // The first bins of the histograms contain the following information: | |
224 | // bin 1: Number of accumulated pulses | |
225 | // bin 2;3;4: Sector; Row; Pad; | |
226 | // | |
227 | ||
228 | TFile fileOut(nameFileOut,"UPDATE"); | |
229 | fileOut.cd(); | |
230 | ||
231 | TH1I *tempHis = new TH1I("tempHis","tempHis",fSample,fGateWidth,fSample+fGateWidth); | |
232 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000); | |
233 | ||
234 | // loop over the data in this event | |
235 | ||
236 | while (rawStream->NextDDL() ) { | |
237 | ||
238 | Int_t ddl = rawReader->GetDDLID(); | |
239 | ||
240 | while (rawStream->NextChannel() ) { | |
241 | ||
242 | while (rawStream->NextBunch() ) { | |
243 | ||
244 | Int_t t0 = rawStream->GetStartTimeBin(); | |
245 | Int_t bl = rawStream->GetBunchLength(); | |
246 | ||
247 | if (bl<fSample+fGateWidth) continue; | |
248 | ||
249 | Int_t sector = rawStream->GetSector(); | |
250 | Int_t row = rawStream->GetRow(); | |
251 | Int_t pad = rawStream->GetPad(); | |
252 | ||
253 | UShort_t *signals=(UShort_t*)rawStream->GetSignals(); | |
254 | if (!signals) continue; | |
255 | ||
256 | // Write to temporary histogramm | |
257 | for (Int_t i=0;i<bl;++i) { | |
258 | UShort_t time=t0-i; | |
259 | UShort_t signal=signals[i]; | |
260 | if ( (fGateWidth<time) && (time<=fSample+fGateWidth) ) { | |
261 | tempHis->SetBinContent(time-fGateWidth,signal); | |
262 | } | |
263 | } | |
264 | ||
265 | // calculation of the pulse properties and comparison to thresholds settings | |
266 | ||
267 | Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX); | |
268 | Int_t maxpos = tempHis->GetMaximumBin(); | |
269 | ||
270 | Int_t first = (Int_t)TMath::Max(maxpos-10, 0); | |
271 | Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth); | |
272 | ||
273 | // simple baseline substraction ? better one needed ? (pedestalsubstr.?) | |
274 | // and RMS calculation with timebins before the pulse and at the end of | |
275 | // the signal | |
276 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
277 | // before the pulse | |
278 | tempRMSHis->Fill(tempHis->GetBinContent(first+ipos)); | |
279 | } | |
280 | for (Int_t ipos = 0; ipos<20; ipos++) { | |
281 | // at the end to get rid of pulses with serious baseline fluctuations | |
282 | tempRMSHis->Fill(tempHis->GetBinContent(last-ipos)); | |
283 | } | |
284 | ||
285 | Double_t baseline = tempRMSHis->GetMean(); | |
286 | Double_t rms = tempRMSHis->GetRMS(); | |
287 | tempRMSHis->Reset(); | |
288 | ||
289 | Double_t lowLim = fLowPulseLim+baseline; | |
290 | Double_t upLim = fUpPulseLim+baseline; | |
291 | ||
292 | // get rid of pulses which contain gate signal and/or too much noise | |
293 | // with the help of ratio of integrals | |
294 | Double_t intHist = 0; | |
295 | Double_t intPulse = 0; | |
296 | Double_t binValue; | |
297 | for(Int_t ipos=first; ipos<=last; ipos++) { | |
298 | binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline); | |
299 | intHist += binValue; | |
300 | if(ipos>=first+5 && ipos<=first+15) {intPulse += binValue;} | |
301 | } | |
302 | ||
303 | // gets rid of high frequency noise: | |
304 | // calculating ratio (value one to the right of maximum)/(maximum) | |
305 | // has to be >= 0.1; if maximum==0 set ratio to 0.1 | |
306 | Double_t maxCorr = max - baseline; | |
307 | Double_t binRatio = 0.1; | |
308 | if(TMath::Abs(maxCorr)>1e-5) { | |
309 | binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr; | |
310 | } | |
311 | ||
312 | // Decision if found pulse is a proper one according to given tresholds | |
313 | if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && (intHist/intPulse)<fRatioIntLim &&intPulse>10&& (binRatio >= 0.1) ) { | |
314 | ||
315 | // 1D histogramm for mean pulse per pad | |
316 | char hname[100]; | |
317 | snprintf(hname,100,"sec%drow%dpad%d",sector,row,pad); | |
318 | ||
319 | TH1F *his = (TH1F*)fileOut.Get(hname); | |
320 | ||
321 | if (!his ) { // new entry (pulse in new pad found) | |
322 | ||
323 | his = new TH1F(hname,hname, fPulseLength+5, 0, fPulseLength+5); | |
324 | his->SetBinContent(1,1); // pulse counter (1st pulse) | |
325 | his->SetBinContent(2,sector); // sector | |
326 | his->SetBinContent(3,row); // row | |
327 | his->SetBinContent(4,pad); // pad | |
328 | ||
329 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
330 | Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); | |
331 | his->SetBinContent(ipos+5,signal); | |
332 | } | |
333 | his->Write(hname); | |
334 | printf("new %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
335 | ||
336 | } else { // adding pulse to existing histogram (pad already found) | |
337 | ||
338 | his->AddBinContent(1,1); // pulse counter for each pad | |
339 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
340 | Int_t signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); | |
341 | his->AddBinContent(ipos+5,signal); | |
342 | } | |
343 | printf("adding ... %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
344 | his->Write(hname,kOverwrite); | |
345 | } | |
346 | ||
347 | ||
348 | // 2D histogramm for pulse spread within a DDL (normalized to one) | |
349 | char hname2d[100]; | |
350 | snprintf(hname2d,100,"2Dhisto_ddl%d",ddl); | |
351 | TH2F *his2d = (TH2F*)fileOut.Get(hname2d); | |
352 | if (!his2d ) { // new entry (ddl was not seen before) | |
353 | ||
354 | his2d = new TH2F(hname2d,hname2d, fPulseLength, 0., (Double_t)fPulseLength, 50,-0.02,0.02); | |
355 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
356 | Double_t signal = tempHis->GetBinContent(ipos+first)-baseline; | |
338e0dd9 | 357 | if (TMath::Abs(signal/maxCorr)>1e-10) // zero bins are biased |
358 | his2d->Fill(ipos,signal/maxCorr); | |
752b0cc7 | 359 | } |
360 | his2d->Write(hname2d); | |
361 | printf("new %s: \n", hname2d); | |
752b0cc7 | 362 | } else { // adding pulse to existing histogram |
363 | ||
364 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
365 | Double_t signal= tempHis->GetBinContent(ipos+first)-baseline; | |
338e0dd9 | 366 | if (TMath::Abs(signal/maxCorr)>1e-10) // zero bins are biased |
367 | his2d->Fill(ipos,signal/maxCorr); | |
752b0cc7 | 368 | } |
369 | his2d->Write(hname2d,kOverwrite); | |
370 | } | |
371 | ||
372 | tempHis->Reset(); | |
373 | ||
374 | } // pulse stored | |
375 | ||
376 | } // bunch loop | |
377 | }// channel loop | |
378 | } // ddl loop | |
379 | ||
380 | tempHis->~TH1I(); | |
381 | tempRMSHis->~TH1I(); | |
382 | printf("Finished to read event ... \n"); | |
383 | fileOut.Close(); | |
384 | ||
385 | } | |
386 | ||
387 | ||
eb7e0771 | 388 | //_____________________________________________________________________________ |
389 | void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nameFileOut) { | |
390 | // | |
391 | // Extracts proper pulses (according the given tresholds) within one event | |
392 | // and accumulates them into one histogram per pad. All histograms are | |
393 | // saved in the file 'nameFileOut'. | |
394 | // The first bins of the histograms contain the following information: | |
395 | // bin 1: Number of accumulated pulses | |
396 | // bin 2;3;4: Sector; Row; Pad; | |
397 | // | |
398 | ||
55f06b51 | 399 | rawStream->Reset(); |
400 | ||
eb7e0771 | 401 | Int_t sector = rawStream->GetSector(); |
402 | Int_t row = rawStream->GetRow(); | |
d0bd4fcc | 403 | |
404 | Int_t prevSec = 999999; | |
405 | Int_t prevRow = 999999; | |
406 | Int_t prevPad = 999999; | |
eb7e0771 | 407 | Int_t prevTime = 999999; |
eb7e0771 | 408 | |
409 | TFile fileOut(nameFileOut,"UPDATE"); | |
410 | fileOut.cd(); | |
411 | ||
412 | TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth); | |
413 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000); | |
414 | ||
55f06b51 | 415 | // printf("raw next: %d\n",rawStream->Next()); |
416 | ||
eb7e0771 | 417 | while (rawStream->Next()) { |
418 | ||
419 | // in case of a new row, get sector and row number | |
420 | if (rawStream->IsNewRow()){ | |
421 | sector = rawStream->GetSector(); | |
422 | row = rawStream->GetRow(); | |
7409c8db | 423 | // if (sector!=prevSec) AliSysInfo::AddStamp(Form("sector_%d_row_%d",sector,row), -1,sector,row); |
eb7e0771 | 424 | } |
425 | ||
426 | Int_t pad = rawStream->GetPad(); | |
427 | Int_t time = rawStream->GetTime(); | |
428 | Int_t signal = rawStream->GetSignal(); | |
55f06b51 | 429 | |
430 | //printf("%d: t:%d sig:%d\n",pad,time,signal); | |
d0bd4fcc | 431 | |
432 | // Reading signal from one Pad | |
433 | if (!rawStream->IsNewPad()) { | |
434 | ||
435 | // this pad always gave a useless signal, probably induced by the supply | |
436 | // voltage of the gate signal (date:2008-Aug-07) | |
437 | if(sector==51 && row==95 && pad==0) { | |
d0bd4fcc | 438 | continue; |
439 | } | |
eb7e0771 | 440 | |
d0bd4fcc | 441 | // only process pulses of pads with correct address |
442 | if(sector<0 || sector+1 > Int_t(AliTPCROC::Instance()->GetNSector())) { | |
d0bd4fcc | 443 | continue; |
444 | } | |
445 | if(row<0 || row+1 > Int_t(AliTPCROC::Instance()->GetNRows(sector))) { | |
d0bd4fcc | 446 | continue; |
447 | } | |
448 | if(pad<0 || pad+1 > Int_t(AliTPCROC::Instance()->GetNPads(sector,row))) { | |
d0bd4fcc | 449 | continue; |
450 | } | |
451 | ||
eb7e0771 | 452 | if (time>prevTime) { |
55f06b51 | 453 | //printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime); |
d0bd4fcc | 454 | continue; |
eb7e0771 | 455 | } else { |
456 | // still the same pad, save signal to temporary histogram | |
8ad79793 | 457 | if ( (time<=fSample+fGateWidth) && (time>=fGateWidth)) { |
eb7e0771 | 458 | tempHis->SetBinContent(time,signal); |
459 | } | |
d0bd4fcc | 460 | } |
461 | ||
eb7e0771 | 462 | } else { |
d0bd4fcc | 463 | |
eb7e0771 | 464 | // complete pulse found and stored into tempHis, now calculation |
d0bd4fcc | 465 | // of the properties and comparison to given thresholds |
55f06b51 | 466 | |
eb7e0771 | 467 | Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX); |
468 | Int_t maxpos = tempHis->GetMaximumBin(); | |
469 | ||
470 | Int_t first = (Int_t)TMath::Max(maxpos-10, 0); | |
8ad79793 | 471 | Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth); |
eb7e0771 | 472 | |
473 | // simple baseline substraction ? better one needed ? (pedestalsubstr.?) | |
474 | // and RMS calculation with timebins before the pulse and at the end of | |
475 | // the signal | |
476 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
477 | // before the pulse | |
478 | tempRMSHis->Fill(tempHis->GetBinContent(first+ipos)); | |
df4cfd77 | 479 | } |
480 | for (Int_t ipos = 0; ipos<20; ipos++) { | |
eb7e0771 | 481 | // at the end to get rid of pulses with serious baseline fluctuations |
482 | tempRMSHis->Fill(tempHis->GetBinContent(last-ipos)); | |
483 | } | |
d0bd4fcc | 484 | |
eb7e0771 | 485 | Double_t baseline = tempRMSHis->GetMean(); |
486 | Double_t rms = tempRMSHis->GetRMS(); | |
487 | tempRMSHis->Reset(); | |
488 | ||
489 | Double_t lowLim = fLowPulseLim+baseline; | |
490 | Double_t upLim = fUpPulseLim+baseline; | |
491 | ||
d0bd4fcc | 492 | // get rid of pulses which contain gate signal and/or too much noise |
493 | // with the help of ratio of integrals | |
494 | Double_t intHist = 0; | |
495 | Double_t intPulse = 0; | |
496 | Double_t binValue; | |
497 | for(Int_t ipos=first; ipos<=last; ipos++) { | |
498 | binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline); | |
499 | intHist += binValue; | |
500 | if(ipos>=first+5 && ipos<=first+25) {intPulse += binValue;} | |
501 | } | |
502 | ||
503 | // gets rid of high frequency noise: | |
504 | // calculating ratio (value one to the right of maximum)/(maximum) | |
505 | // has to be >= 0.1; if maximum==0 set ratio to 0.1 | |
506 | Double_t maxCorr = max - baseline; | |
507 | Double_t binRatio = 0.1; | |
56c85970 | 508 | if(TMath::Abs(maxCorr)>1e-5) { |
d0bd4fcc | 509 | binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr; |
510 | } | |
511 | ||
eb7e0771 | 512 | // Decision if found pulse is a proper one according to given tresholds |
8ad79793 | 513 | if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && (intHist/intPulse)<fRatioIntLim && (binRatio >= 0.1) ) { |
eb7e0771 | 514 | char hname[100]; |
56c85970 | 515 | snprintf(hname,100,"sec%drow%dpad%d",prevSec,prevRow,prevPad); |
eb7e0771 | 516 | |
517 | TH1F *his = (TH1F*)fileOut.Get(hname); | |
518 | ||
519 | if (!his ) { // new entry (pulse in new pad found) | |
520 | ||
521 | his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4); | |
7409c8db | 522 | his->SetBinContent(1,1); // pulse counter (1st pulse) |
d0bd4fcc | 523 | his->SetBinContent(2,prevSec); // sector |
524 | his->SetBinContent(3,prevRow); // row | |
525 | his->SetBinContent(4,prevPad); // pad | |
eb7e0771 | 526 | |
527 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
7409c8db | 528 | signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); |
eb7e0771 | 529 | his->SetBinContent(ipos+5,signal); |
530 | } | |
531 | his->Write(hname); | |
532 | printf("new %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
533 | ||
534 | } else { // adding pulse to existing histogram (pad already found) | |
535 | ||
536 | his->AddBinContent(1,1); // pulse counter for each pad | |
537 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
7409c8db | 538 | signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); |
eb7e0771 | 539 | his->AddBinContent(ipos+5,signal); |
540 | } | |
541 | printf("adding ... %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
542 | his->Write(hname,kOverwrite); | |
543 | } | |
544 | } | |
545 | tempHis->Reset(); | |
546 | } | |
547 | prevTime = time; | |
d0bd4fcc | 548 | prevSec = sector; |
549 | prevRow = row; | |
eb7e0771 | 550 | prevPad = pad; |
551 | } | |
552 | ||
553 | tempHis->~TH1I(); | |
554 | tempRMSHis->~TH1I(); | |
555 | printf("Finished to read event ... \n"); | |
556 | fileOut.Close(); | |
557 | } | |
558 | ||
559 | //____________________________________________________________________________ | |
560 | void AliTPCCalibTCF::MergeHistoPerSector(const char *nameFileIn) { | |
561 | // | |
562 | // Merges all histograms within one sector, calculates the TCF parameters | |
563 | // of the 'histogram-per-sector' and stores (histo and parameters) into | |
564 | // seperated files ... | |
565 | // | |
566 | // note: first 4 timebins of a histogram hold specific informations | |
567 | // about number of collected pulses, sector, row and pad | |
568 | // | |
569 | // 'nameFileIn': root file produced with Process function which holds | |
570 | // one histogram per pad (sum of signals of proper pulses) | |
571 | // 'Sec+nameFileIn': root file with one histogram per sector | |
572 | // (information of row and pad are set to -1) | |
573 | // | |
574 | ||
575 | TFile fileIn(nameFileIn,"READ"); | |
576 | TH1F *hisPad = 0; | |
577 | TKey *key = 0; | |
578 | TIter next( fileIn.GetListOfKeys() ); | |
579 | ||
580 | char nameFileOut[100]; | |
56c85970 | 581 | snprintf(nameFileOut,100,"Sec-%s",nameFileIn); |
eb7e0771 | 582 | |
583 | TFile fileOut(nameFileOut,"RECREATE"); | |
584 | fileOut.cd(); | |
585 | ||
586 | Int_t nHist = fileIn.GetNkeys(); | |
587 | Int_t iHist = 0; // histogram counter for merge-status print | |
588 | ||
589 | while ( (key=(TKey*)next()) ) { | |
590 | ||
591 | iHist++; | |
752b0cc7 | 592 | TString name(key->GetName()); |
593 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
594 | ||
595 | hisPad = (TH1F*)fileIn.Get(name.Data()); // copy object to memory | |
eb7e0771 | 596 | |
eb7e0771 | 597 | Int_t pulseLength = hisPad->GetNbinsX() -4; |
598 | // -4 because first four timebins contain pad specific informations | |
599 | Int_t npulse = (Int_t)hisPad->GetBinContent(1); | |
600 | Int_t sector = (Int_t)hisPad->GetBinContent(2); | |
601 | ||
602 | char hname[100]; | |
56c85970 | 603 | snprintf(hname,100,"sector%d",sector); |
eb7e0771 | 604 | TH1F *his = (TH1F*)fileOut.Get(hname); |
605 | ||
606 | if (!his ) { // new histogram (new sector) | |
607 | his = new TH1F(hname,hname, pulseLength+4, 0, pulseLength+4); | |
608 | his->SetBinContent(1,npulse); // pulse counter | |
609 | his->SetBinContent(2,sector); // set sector info | |
610 | his->SetBinContent(3,-1); // set to dummy value | |
611 | his->SetBinContent(4,-1); // set to dummy value | |
612 | for (Int_t ipos=0; ipos<pulseLength; ipos++){ | |
613 | his->SetBinContent(ipos+5,hisPad->GetBinContent(ipos+5)); | |
614 | } | |
615 | his->Write(hname); | |
616 | printf("found %s ...\n", hname); | |
617 | } else { // add to existing histogram for sector | |
618 | his->AddBinContent(1,npulse); // pulse counter | |
619 | for (Int_t ipos=0; ipos<pulseLength; ipos++){ | |
620 | his->AddBinContent(ipos+5,hisPad->GetBinContent(ipos+5)); | |
621 | } | |
752b0cc7 | 622 | his->Write(hname,kOverwrite); |
eb7e0771 | 623 | } |
624 | ||
625 | if (iHist%500==0) { | |
626 | printf("merging status: \t %d pads out of %d \n",iHist, nHist); | |
627 | } | |
628 | } | |
df4cfd77 | 629 | |
eb7e0771 | 630 | printf("merging done ...\n"); |
631 | fileIn.Close(); | |
632 | fileOut.Close(); | |
633 | ||
eb7e0771 | 634 | |
635 | } | |
636 | ||
637 | ||
638 | //____________________________________________________________________________ | |
d0bd4fcc | 639 | void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse, Int_t histStart, Int_t histEnd) { |
eb7e0771 | 640 | // |
641 | // This function takes a prepeared root file (accumulated histograms: output | |
642 | // of process function) and performs an analysis (fit and equalization) in | |
643 | // order to get the TCF parameters. These are stored in an TNtuple along with | |
644 | // the pad and creation infos. The tuple is written to the output file | |
645 | // "TCFparam+nameFileIn" | |
646 | // To reduce the analysis time, the minimum number of accumulated pulses within | |
647 | // one histogram 'minNumPulse' (to perform the analysis on) can be set | |
648 | // | |
649 | ||
650 | TFile fileIn(nameFileIn,"READ"); | |
651 | TH1F *hisIn; | |
652 | TKey *key; | |
653 | TIter next( fileIn.GetListOfKeys() ); | |
654 | ||
655 | char nameFileOut[100]; | |
56c85970 | 656 | snprintf(nameFileOut,100,"TCF-%s",nameFileIn); |
eb7e0771 | 657 | |
658 | TFile fileOut(nameFileOut,"RECREATE"); | |
659 | fileOut.cd(); | |
660 | ||
661 | TNtuple *paramTuple = new TNtuple("TCFparam","TCFparameter","sec:row:pad:npulse:Z0:Z1:Z2:P0:P1:P2"); | |
662 | ||
663 | Int_t nHist = fileIn.GetNkeys(); | |
664 | Int_t iHist = 0; // counter for print of analysis-status | |
665 | ||
9389f9a4 | 666 | while ((key = (TKey *) next())) { // loop over histograms |
d0bd4fcc | 667 | ++iHist; |
668 | if(iHist < histStart || iHist > histEnd) {continue;} | |
752b0cc7 | 669 | |
670 | TString name(key->GetName()); | |
671 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
672 | ||
eb7e0771 | 673 | hisIn = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory |
338e0dd9 | 674 | |
eb7e0771 | 675 | Int_t numPulse = (Int_t)hisIn->GetBinContent(1); |
676 | if ( numPulse >= minNumPulse ) { | |
df4cfd77 | 677 | printf("Analyze histogram %d out of %d\n",iHist,nHist); |
eb7e0771 | 678 | Double_t* coefP = new Double_t[3]; |
679 | Double_t* coefZ = new Double_t[3]; | |
680 | for(Int_t i = 0; i < 3; i++){ | |
681 | coefP[i] = 0; | |
682 | coefZ[i] = 0; | |
683 | } | |
684 | // perform the analysis on the given histogram | |
685 | Int_t fitOk = AnalyzePulse(hisIn, coefZ, coefP); | |
686 | if (fitOk) { // Add found parameters to file | |
687 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
688 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
689 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
690 | paramTuple->Fill(sector,row,pad,numPulse,coefZ[0],coefZ[1],coefZ[2],coefP[0],coefP[1],coefP[2]); | |
691 | } | |
692 | coefP->~Double_t(); | |
693 | coefZ->~Double_t(); | |
df4cfd77 | 694 | } else { |
695 | printf("Skip histogram %d out of %d | not enough accumulated pulses\n",iHist,nHist); | |
eb7e0771 | 696 | } |
df4cfd77 | 697 | |
eb7e0771 | 698 | } |
699 | ||
700 | fileIn.Close(); | |
701 | paramTuple->Write(); | |
702 | fileOut.Close(); | |
703 | ||
704 | } | |
705 | ||
706 | ||
707 | //____________________________________________________________________________ | |
56c85970 | 708 | Int_t AliTPCCalibTCF::AnalyzePulse(TH1F * const hisIn, Double_t *coefZ, Double_t *coefP) { |
eb7e0771 | 709 | // |
710 | // Performs the analysis on one specific pulse (histogram) by means of fitting | |
711 | // the pulse and equalization of the pulseheight. The found TCF parameters | |
712 | // are stored in the arrays coefZ and coefP | |
713 | // | |
714 | ||
715 | Int_t pulseLength = hisIn->GetNbinsX() -4; | |
d0bd4fcc | 716 | // -4 because the first four timebins usually contain pad specific informations |
eb7e0771 | 717 | Int_t npulse = (Int_t)hisIn->GetBinContent(1); |
718 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
719 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
720 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
721 | ||
722 | // write pulseinformation to TNtuple and normalize to 100 ADC (because of | |
723 | // given upper and lower fit parameter limits) in order to pass the pulse | |
724 | // to TMinuit | |
725 | ||
726 | TNtuple *dataTuple = new TNtuple("ntupleFit","Pulse","timebin:sigNorm:error"); | |
727 | Double_t error = 0.05; | |
728 | Double_t max = hisIn->GetMaximum(FLT_MAX); | |
729 | for (Int_t ipos=0; ipos<pulseLength; ipos++) { | |
730 | Double_t errorz=error; | |
731 | if (ipos>100) { errorz = error*100; } // very simple weight: FIXME in case | |
732 | Double_t signal = hisIn->GetBinContent(ipos+5); | |
733 | Double_t signalNorm = signal/max*100; //pulseheight normaliz. to 100ADC | |
734 | dataTuple->Fill(ipos, signalNorm, errorz); | |
735 | } | |
736 | ||
737 | // Call fit function (TMinuit) to get the first 2 PZ Values for the | |
738 | // Tail Cancelation Filter | |
739 | Int_t fitOk = FitPulse(dataTuple, coefZ, coefP); | |
740 | ||
741 | if (fitOk) { | |
742 | // calculates the 3rd set (remaining 2 PZ values) in order to restore the | |
743 | // original height of the pulse | |
d0bd4fcc | 744 | Int_t equOk = Equalization(dataTuple, coefZ, coefP); |
745 | if (!equOk) { | |
746 | Error("FindFit", "Pulse equalisation procedure failed - pulse abandoned "); | |
747 | printf("in Sector %d | Row %d | Pad %d |", sector, row, pad); | |
748 | printf(" Npulses: %d \n\n", npulse); | |
749 | coefP[2] = 0; coefZ[2] = 0; | |
750 | dataTuple->~TNtuple(); | |
751 | return 0; | |
752 | } | |
eb7e0771 | 753 | printf("Calculated TCF parameters for: \n"); |
754 | printf("Sector %d | Row %d | Pad %d |", sector, row, pad); | |
755 | printf(" Npulses: %d \n", npulse); | |
756 | for(Int_t i = 0; i < 3; i++){ | |
757 | printf("P[%d] = %f Z[%d] = %f \n",i,coefP[i],i,coefZ[i]); | |
758 | if (i==2) { printf("\n"); } | |
759 | } | |
760 | dataTuple->~TNtuple(); | |
761 | return 1; | |
762 | } else { // fit did not converge | |
763 | Error("FindFit", "TCF fit not converged - pulse abandoned "); | |
764 | printf("in Sector %d | Row %d | Pad %d |", sector, row, pad); | |
765 | printf(" Npulses: %d \n\n", npulse); | |
766 | coefP[2] = 0; coefZ[2] = 0; | |
767 | dataTuple->~TNtuple(); | |
768 | return 0; | |
769 | } | |
770 | ||
771 | } | |
772 | ||
773 | ||
774 | ||
775 | //____________________________________________________________________________ | |
df4cfd77 | 776 | void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameFileTCF, Int_t minNumPulse, Int_t plotFlag, Int_t lowKey, Int_t upKey) |
eb7e0771 | 777 | { |
778 | // | |
779 | // Performs quality parameters evaluation of the calculated TCF parameters in | |
780 | // the file 'nameFileTCF' for every (accumulated) histogram within the | |
781 | // prepeared root file 'nameFileIn'. | |
782 | // The found quality parameters are stored in an TNtuple which will be saved | |
783 | // in a Root file 'Quality-*'. | |
784 | // If the parameter for the given pulse (given pad) was not found, the pulse | |
785 | // is rejected. | |
786 | // | |
787 | ||
788 | TFile fileIn(nameFileIn,"READ"); | |
789 | ||
790 | Double_t* coefP = new Double_t[3]; | |
791 | Double_t* coefZ = new Double_t[3]; | |
792 | for(Int_t i = 0; i < 3; i++){ | |
793 | coefP[i] = 0; | |
794 | coefZ[i] = 0; | |
795 | } | |
796 | ||
797 | char nameFileOut[100]; | |
56c85970 | 798 | snprintf(nameFileOut,100,"Quality_%s_AT_%s",nameFileTCF, nameFileIn); |
eb7e0771 | 799 | TFile fileOut(nameFileOut,"RECREATE"); |
800 | ||
801 | TNtuple *qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot"); | |
802 | ||
803 | TH1F *hisIn; | |
804 | TKey *key; | |
805 | TIter next( fileIn.GetListOfKeys() ); | |
806 | ||
807 | Int_t nHist = fileIn.GetNkeys(); | |
808 | Int_t iHist = 0; | |
809 | ||
810 | for(Int_t i=0;i<lowKey-1;i++){++iHist; key = (TKey *) next();} | |
2c632057 | 811 | while ((key = (TKey *) next())) { // loop over saved histograms |
eb7e0771 | 812 | |
813 | // loading pulse to memory; | |
752b0cc7 | 814 | TString name(key->GetName()); |
815 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
816 | ||
eb7e0771 | 817 | printf("validating pulse %d out of %d\n",++iHist,nHist); |
818 | hisIn = (TH1F*)fileIn.Get(key->GetName()); | |
338e0dd9 | 819 | |
eb7e0771 | 820 | |
821 | // find the correct TCF parameter according to the his infos (first 4 bins) | |
822 | Int_t nPulse = FindCorTCFparam(hisIn, nameFileTCF, coefZ, coefP); | |
df4cfd77 | 823 | if (nPulse>=minNumPulse) { // doing the TCF quality analysis |
eb7e0771 | 824 | Double_t *quVal = GetQualityOfTCF(hisIn,coefZ,coefP, plotFlag); |
825 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
826 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
827 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
828 | qualityTuple->Fill(sector,row,pad,nPulse,quVal[0],quVal[1],quVal[2],quVal[3],quVal[4],quVal[5]); | |
829 | quVal->~Double_t(); | |
830 | } | |
831 | ||
832 | if (iHist>=upKey) {break;} | |
833 | ||
834 | } | |
835 | ||
836 | fileOut.cd(); | |
837 | qualityTuple->Write(); | |
838 | ||
839 | coefP->~Double_t(); | |
840 | coefZ->~Double_t(); | |
841 | ||
842 | fileOut.Close(); | |
843 | fileIn.Close(); | |
844 | ||
845 | } | |
846 | ||
847 | ||
848 | ||
849 | //_____________________________________________________________________________ | |
55f06b51 | 850 | void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameFileOut, const char *nameFileTCF, Int_t minNumPulse, Int_t plotFlag, bool bUseHLTOUT) { |
eb7e0771 | 851 | // |
852 | // Performs quality parameters evaluation of the calculated TCF parameters in | |
853 | // the file 'nameFileTCF' for every proper pulse (according to given thresholds) | |
854 | // within the RAW file 'nameRawFile'. | |
855 | // The found quality parameters are stored in a TNtuple which will be saved | |
856 | // in the Root file 'nameFileOut'. If the parameter for the given pulse | |
857 | // (given pad) was not found, the pulse is rejected. | |
858 | // | |
859 | ||
860 | // | |
861 | // Reads a RAW data file, extracts Pulses (according the given tresholds) | |
862 | // and test the found TCF parameters on them ... | |
863 | // | |
864 | ||
55f06b51 | 865 | |
866 | // create the data reader | |
eb7e0771 | 867 | AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile); |
55f06b51 | 868 | if (!rawReader) { |
869 | return; | |
870 | } | |
871 | ||
872 | // create HLT reader for redirection of TPC data from HLTOUT to TPC reconstruction | |
873 | AliRawReader *hltReader=AliRawHLTManager::AliRawHLTManager::CreateRawReaderHLT(rawReader, "TPC"); | |
874 | ||
875 | // now choose the data source | |
876 | if (bUseHLTOUT) rawReader=hltReader; | |
877 | ||
878 | // rawReader->Reset(); | |
879 | rawReader->RewindEvents(); | |
880 | ||
881 | if (!rawReader->NextEvent()) { | |
882 | printf("no events found in %s\n",nameRawFile); | |
883 | return; | |
884 | } | |
885 | ||
eb7e0771 | 886 | Double_t* coefP = new Double_t[3]; |
887 | Double_t* coefZ = new Double_t[3]; | |
888 | for(Int_t i = 0; i < 3; i++){ | |
889 | coefP[i] = 0; | |
890 | coefZ[i] = 0; | |
891 | } | |
892 | ||
d0bd4fcc | 893 | Int_t ievent = 0; |
894 | ||
895 | TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth); | |
896 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000); | |
897 | ||
898 | TFile fileOut(nameFileOut,"UPDATE"); // Quality Parameters storage | |
899 | TNtuple *qualityTuple = (TNtuple*)fileOut.Get("TCFquality"); | |
900 | if (!qualityTuple) { // no entry in file | |
901 | qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot:pulseRMS"); | |
902 | } | |
903 | ||
55f06b51 | 904 | do { |
eb7e0771 | 905 | |
d0bd4fcc | 906 | printf("Reading next event ... Nr:%d\n",ievent); |
7409c8db | 907 | AliTPCRawStream *rawStream = new AliTPCRawStream(rawReader); |
eb7e0771 | 908 | rawReader->Select("TPC"); |
d0bd4fcc | 909 | ievent++; |
eb7e0771 | 910 | |
7409c8db | 911 | Int_t sector = rawStream->GetSector(); |
912 | Int_t row = rawStream->GetRow(); | |
d0bd4fcc | 913 | |
914 | Int_t prevSec = 999999; | |
915 | Int_t prevRow = 999999; | |
916 | Int_t prevPad = 999999; | |
eb7e0771 | 917 | Int_t prevTime = 999999; |
eb7e0771 | 918 | |
7409c8db | 919 | while (rawStream->Next()) { |
eb7e0771 | 920 | |
7409c8db | 921 | if (rawStream->IsNewRow()){ |
922 | sector = rawStream->GetSector(); | |
923 | row = rawStream->GetRow(); | |
eb7e0771 | 924 | } |
925 | ||
7409c8db | 926 | Int_t pad = rawStream->GetPad(); |
927 | Int_t time = rawStream->GetTime(); | |
928 | Int_t signal = rawStream->GetSignal(); | |
eb7e0771 | 929 | |
7409c8db | 930 | if (!rawStream->IsNewPad()) { // Reading signal from one Pad |
d0bd4fcc | 931 | |
932 | // this pad always gave a useless signal, probably induced by the supply | |
933 | // voltage of the gate signal (date:2008-Aug-07) | |
934 | if(sector==51 && row==95 && pad==0) { | |
d0bd4fcc | 935 | continue; |
936 | } | |
937 | ||
938 | // only process pulses of pads with correct address | |
939 | if(sector<0 || sector+1 > Int_t(AliTPCROC::Instance()->GetNSector())) { | |
d0bd4fcc | 940 | continue; |
941 | } | |
942 | if(row<0 || row+1 > Int_t(AliTPCROC::Instance()->GetNRows(sector))) { | |
d0bd4fcc | 943 | continue; |
944 | } | |
945 | if(pad<0 || pad+1 > Int_t(AliTPCROC::Instance()->GetNPads(sector,row))) { | |
d0bd4fcc | 946 | continue; |
947 | } | |
948 | ||
eb7e0771 | 949 | if (time>prevTime) { |
7409c8db | 950 | // printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime); |
d0bd4fcc | 951 | continue; |
eb7e0771 | 952 | } else { |
d0bd4fcc | 953 | // still the same pad, save signal to temporary histogram |
eb7e0771 | 954 | if (time<=fSample+fGateWidth && time>fGateWidth) { |
955 | tempHis->SetBinContent(time,signal); | |
956 | } | |
957 | } | |
958 | } else { // Decision for saving pulse according to treshold settings | |
959 | ||
960 | Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX); | |
961 | Int_t maxpos = tempHis->GetMaximumBin(); | |
962 | ||
963 | Int_t first = (Int_t)TMath::Max(maxpos-10, 0); | |
8ad79793 | 964 | Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample+fGateWidth); |
eb7e0771 | 965 | |
966 | ||
967 | // simple baseline substraction ? better one needed ? (pedestalsubstr.?) | |
968 | // and RMS calculation with timebins before the pulse and at the end of | |
969 | // the signal | |
970 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
971 | // before the pulse | |
972 | tempRMSHis->Fill(tempHis->GetBinContent(first+ipos)); | |
df4cfd77 | 973 | } |
974 | for (Int_t ipos = 0; ipos<20; ipos++) { | |
eb7e0771 | 975 | // at the end to get rid of pulses with serious baseline fluctuations |
df4cfd77 | 976 | tempRMSHis->Fill(tempHis->GetBinContent(last-ipos)); |
eb7e0771 | 977 | } |
978 | Double_t baseline = tempRMSHis->GetMean(); | |
979 | Double_t rms = tempRMSHis->GetRMS(); | |
980 | tempRMSHis->Reset(); | |
981 | ||
982 | Double_t lowLim = fLowPulseLim+baseline; | |
983 | Double_t upLim = fUpPulseLim+baseline; | |
984 | ||
d0bd4fcc | 985 | // get rid of pulses which contain gate signal and/or too much noise |
986 | // with the help of ratio of integrals | |
987 | Double_t intHist = 0; | |
988 | Double_t intPulse = 0; | |
989 | Double_t binValue; | |
990 | for(Int_t ipos=first; ipos<=last; ipos++) { | |
991 | binValue = TMath::Abs(tempHis->GetBinContent(ipos) - baseline); | |
992 | intHist += binValue; | |
752b0cc7 | 993 | if(ipos>=first+5 && ipos<=first+15) {intPulse += binValue;} |
d0bd4fcc | 994 | } |
995 | ||
996 | // gets rid of high frequency noise: | |
997 | // calculating ratio (value one to the right of maximum)/(maximum) | |
998 | // has to be >= 0.1; if maximum==0 set ratio to 0.1 | |
999 | Double_t maxCorr = max - baseline; | |
1000 | Double_t binRatio = 0.1; | |
56c85970 | 1001 | if(TMath::Abs(maxCorr) > 1e-5 ) { |
d0bd4fcc | 1002 | binRatio = (tempHis->GetBinContent(maxpos+1) - baseline) / maxCorr; |
1003 | } | |
1004 | ||
1005 | ||
eb7e0771 | 1006 | // Decision if found pulse is a proper one according to given tresholds |
8ad79793 | 1007 | if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim && intHist/intPulse<fRatioIntLim && (binRatio >= 0.1) ){ |
eb7e0771 | 1008 | // note: |
1009 | // assuming that lowLim is higher than the pedestal value! | |
1010 | char hname[100]; | |
56c85970 | 1011 | snprintf(hname,100,"sec%drow%dpad%d",prevSec,prevRow,prevPad); |
eb7e0771 | 1012 | TH1F *his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4); |
1013 | his->SetBinContent(1,1); // pulse counter (1st pulse) | |
d0bd4fcc | 1014 | his->SetBinContent(2,prevSec); // sector |
1015 | his->SetBinContent(3,prevRow); // row | |
1016 | his->SetBinContent(4,prevPad); // pad | |
1017 | ||
eb7e0771 | 1018 | for (Int_t ipos=0; ipos<last-first; ipos++){ |
7409c8db | 1019 | signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); |
eb7e0771 | 1020 | his->SetBinContent(ipos+5,signal); |
1021 | } | |
1022 | ||
1023 | printf("Pulse found in %s: ADC %d at bin %d \n", hname, max, maxpos+fGateWidth); | |
1024 | ||
1025 | // find the correct TCF parameter according to the his infos | |
1026 | // (first 4 bins) | |
1027 | Int_t nPulse = FindCorTCFparam(his, nameFileTCF, coefZ, coefP); | |
1028 | ||
df4cfd77 | 1029 | if (nPulse>=minNumPulse) { // Parameters found - doing the TCF quality analysis |
eb7e0771 | 1030 | Double_t *quVal = GetQualityOfTCF(his,coefZ,coefP, plotFlag); |
1031 | qualityTuple->Fill(sector,row,pad,nPulse,quVal[0],quVal[1],quVal[2],quVal[3],quVal[4],quVal[5]); | |
1032 | quVal->~Double_t(); | |
1033 | } | |
1034 | his->~TH1F(); | |
1035 | } | |
1036 | tempHis->Reset(); | |
1037 | } | |
1038 | prevTime = time; | |
d0bd4fcc | 1039 | prevSec = sector; |
1040 | prevRow = row; | |
eb7e0771 | 1041 | prevPad = pad; |
d0bd4fcc | 1042 | |
eb7e0771 | 1043 | } |
1044 | ||
d0bd4fcc | 1045 | printf("Finished to read event ... \n"); |
eb7e0771 | 1046 | |
7409c8db | 1047 | delete rawStream; |
1048 | ||
1049 | ||
1050 | } while (rawReader->NextEvent()); // event loop | |
eb7e0771 | 1051 | |
d0bd4fcc | 1052 | printf("Finished to read file - close output file ... \n"); |
1053 | ||
1054 | fileOut.cd(); | |
1055 | qualityTuple->Write("TCFquality",kOverwrite); | |
1056 | fileOut.Close(); | |
1057 | ||
1058 | tempHis->~TH1I(); | |
1059 | tempRMSHis->~TH1I(); | |
eb7e0771 | 1060 | |
1061 | coefP->~Double_t(); | |
1062 | coefZ->~Double_t(); | |
1063 | ||
1064 | rawReader->~AliRawReader(); | |
1065 | ||
1066 | } | |
1067 | ||
df4cfd77 | 1068 | //____________________________________________________________________________ |
1069 | TH2F *AliTPCCalibTCF::PlotOccupSummary2Dhist(const char *nameFileIn, Int_t side) { | |
1070 | // | |
1071 | // Plots the number of summed pulses per pad on a given TPC side | |
1072 | // 'nameFileIn': root-file created with the Process function | |
1073 | // | |
1074 | ||
1075 | TFile fileIn(nameFileIn,"READ"); | |
1076 | TH1F *his; | |
1077 | TKey *key; | |
1078 | TIter next(fileIn.GetListOfKeys()); | |
1079 | ||
1080 | TH2F * his2D = new TH2F("his2D","his2D", 250,-250,250,250,-250,250); | |
8ad79793 | 1081 | |
df4cfd77 | 1082 | AliTPCROC * roc = AliTPCROC::Instance(); |
1083 | ||
1084 | Int_t nHist=fileIn.GetNkeys(); | |
8ad79793 | 1085 | if (!nHist) { return 0; } |
1086 | ||
df4cfd77 | 1087 | Int_t iHist = 0; |
1088 | Float_t xyz[3]; | |
1089 | ||
1090 | Int_t binx = 0; | |
1091 | Int_t biny = 0; | |
1092 | ||
1093 | Int_t npulse = 0; | |
1094 | Int_t sec = 0; | |
1095 | Int_t row = 0; | |
1096 | Int_t pad = 0; | |
1097 | ||
1098 | while ((key = (TKey *) next())) { // loop over histograms within the file | |
8ad79793 | 1099 | iHist++; |
752b0cc7 | 1100 | |
1101 | TString name(key->GetName()); | |
1102 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
1103 | ||
df4cfd77 | 1104 | his = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory |
1105 | ||
1106 | npulse = (Int_t)his->GetBinContent(1); | |
1107 | sec = (Int_t)his->GetBinContent(2); | |
1108 | row = (Int_t)his->GetBinContent(3); | |
1109 | pad = (Int_t)his->GetBinContent(4); | |
1110 | ||
8ad79793 | 1111 | if ( (side==0) && (sec%36>=18) ) continue; |
1112 | if ( (side>0) && (sec%36<18) ) continue; | |
df4cfd77 | 1113 | |
56c85970 | 1114 | if ( (row<0) && (pad<0) ) { // row and pad are equal to -1, then -> summed pulses per sector |
df4cfd77 | 1115 | // fill all pad with this values |
7409c8db | 1116 | for (UInt_t rowi=0; rowi<roc->GetNRows(sec); rowi++) { |
8ad79793 | 1117 | for (UInt_t padi=0; padi<roc->GetNPads(sec,rowi); padi++) { |
7409c8db | 1118 | roc->GetPositionGlobal(sec,rowi,padi,xyz); |
df4cfd77 | 1119 | binx = 1+TMath::Nint((xyz[0]+250.)*0.5); |
1120 | biny = 1+TMath::Nint((xyz[1]+250.)*0.5); | |
1121 | his2D->SetBinContent(binx,biny,npulse); | |
1122 | } | |
1123 | } | |
1124 | } else { | |
1125 | roc->GetPositionGlobal(sec,row,pad,xyz); | |
1126 | binx = 1+TMath::Nint((xyz[0]+250.)*0.5); | |
1127 | biny = 1+TMath::Nint((xyz[1]+250.)*0.5); | |
1128 | ||
1129 | his2D->SetBinContent(binx,biny,npulse); | |
1130 | } | |
1131 | if (iHist%100==0){ printf("hist %d out of %d\n",iHist,nHist);} | |
1132 | } | |
1133 | his2D->SetXTitle("x (cm)"); | |
1134 | his2D->SetYTitle("y (cm)"); | |
8ad79793 | 1135 | his2D->SetStats(0); |
1136 | ||
1137 | his2D->DrawCopy("colz"); | |
df4cfd77 | 1138 | |
1139 | if (!side) { | |
1140 | gPad->SetTitle("A side"); | |
1141 | } else { | |
1142 | gPad->SetTitle("C side"); | |
1143 | } | |
1144 | ||
df4cfd77 | 1145 | return his2D; |
1146 | } | |
1147 | ||
eb7e0771 | 1148 | |
1149 | //____________________________________________________________________________ | |
df4cfd77 | 1150 | void AliTPCCalibTCF::PlotOccupSummary(const char *nameFile, Int_t side, Int_t nPulseMin) { |
eb7e0771 | 1151 | // |
1152 | // Plots the number of summed pulses per pad above a given minimum at the | |
df4cfd77 | 1153 | // pad position at a given TPC side |
eb7e0771 | 1154 | // 'nameFile': root-file created with the Process function |
1155 | // | |
1156 | ||
1157 | TFile *file = new TFile(nameFile,"READ"); | |
eb7e0771 | 1158 | TH1F *his; |
1159 | TKey *key; | |
1160 | TIter next( file->GetListOfKeys() ); | |
1161 | ||
df4cfd77 | 1162 | |
1163 | char nameFileOut[100]; | |
56c85970 | 1164 | snprintf(nameFileOut,100,"Occup-%s",nameFile); |
df4cfd77 | 1165 | TFile fileOut(nameFileOut,"RECREATE"); |
8ad79793 | 1166 | // fileOut.cd(); |
df4cfd77 | 1167 | |
eb7e0771 | 1168 | TNtuple *ntuple = new TNtuple("ntuple","ntuple","x:y:z:npulse"); |
8ad79793 | 1169 | // ntuple->SetDirectory(0); // force to be memory resistent |
eb7e0771 | 1170 | |
df4cfd77 | 1171 | Int_t nHist=file->GetNkeys(); |
8ad79793 | 1172 | if (!nHist) { return; } |
df4cfd77 | 1173 | Int_t iHist = 0; |
7409c8db | 1174 | |
1175 | Int_t secWise = 0; | |
1176 | ||
eb7e0771 | 1177 | while ((key = (TKey *) next())) { // loop over histograms within the file |
752b0cc7 | 1178 | |
1179 | TString name(key->GetName()); | |
1180 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
1181 | ||
eb7e0771 | 1182 | his = (TH1F*)file->Get(key->GetName()); // copy object to memory |
8ad79793 | 1183 | iHist++; |
eb7e0771 | 1184 | Int_t npulse = (Int_t)his->GetBinContent(1); |
1185 | Int_t sec = (Int_t)his->GetBinContent(2); | |
1186 | Int_t row = (Int_t)his->GetBinContent(3); | |
1187 | Int_t pad = (Int_t)his->GetBinContent(4); | |
1188 | ||
56c85970 | 1189 | if ( (row<0) && (pad<0) ) { // row and pad are equal to -1, then -> summed pulses per sector |
eb7e0771 | 1190 | row = 40; pad = 40; // set to approx middle row for better plot |
7409c8db | 1191 | secWise=1; |
eb7e0771 | 1192 | } |
1193 | ||
1194 | Float_t *pos = new Float_t[3]; | |
1195 | // find x,y,z position of the pad | |
1196 | AliTPCROC::Instance()->GetPositionGlobal(sec,row,pad,pos); | |
1197 | if (npulse>=nPulseMin) { | |
1198 | ntuple->Fill(pos[0],pos[1],pos[2],npulse); | |
df4cfd77 | 1199 | if (iHist%100==0){ printf("hist %d out of %d\n",iHist,nHist);} |
eb7e0771 | 1200 | } |
1201 | pos->~Float_t(); | |
eb7e0771 | 1202 | } |
eb7e0771 | 1203 | |
7409c8db | 1204 | if (secWise) { // pulse per sector |
eb7e0771 | 1205 | ntuple->SetMarkerStyle(8); |
1206 | ntuple->SetMarkerSize(4); | |
df4cfd77 | 1207 | } else { // pulse per Pad |
eb7e0771 | 1208 | ntuple->SetMarkerStyle(7); |
1209 | } | |
1210 | ||
df4cfd77 | 1211 | char cSel[100]; |
1212 | if (!side) { | |
56c85970 | 1213 | snprintf(cSel,100,"z>0&&npulse>=%d",nPulseMin); |
df4cfd77 | 1214 | ntuple->Draw("y:x:npulse",cSel,"colz"); |
df4cfd77 | 1215 | } else { |
56c85970 | 1216 | snprintf(cSel,100,"z<0&&npulse>=%d",nPulseMin); |
df4cfd77 | 1217 | ntuple->Draw("y:x:npulse",cSel,"colz"); |
8ad79793 | 1218 | } |
1219 | ||
1220 | if (!side) { | |
1221 | gPad->SetTitle("A side"); | |
1222 | } else { | |
df4cfd77 | 1223 | gPad->SetTitle("C side"); |
1224 | } | |
eb7e0771 | 1225 | |
8ad79793 | 1226 | |
df4cfd77 | 1227 | ntuple->Write(); |
1228 | fileOut.Close(); | |
eb7e0771 | 1229 | file->Close(); |
eb7e0771 | 1230 | } |
1231 | ||
1232 | //____________________________________________________________________________ | |
1233 | void AliTPCCalibTCF::PlotQualitySummary(const char *nameFileQuality, const char *plotSpec, const char *cut, const char *pOpt) | |
1234 | { | |
1235 | // | |
1236 | // This function is an easy interface to load the QualityTuple (produced with | |
1237 | // the function 'TestOn%File' and plots them according to the plot specifications | |
1238 | // 'plotSpec' e.g. "widthRed:maxUndershot" | |
1239 | // One may also set cut and plot options ("cut","pOpt") | |
1240 | // | |
1241 | // The stored quality parameters are ... | |
1242 | // sec:row:pad:npulse: ... usual pad info | |
1243 | // heightDev ... height deviation in percent | |
1244 | // areaRed ... area reduction in percent | |
1245 | // widthRed ... width reduction in percent | |
1246 | // undershot ... mean undershot after the pulse in ADC | |
1247 | // maxUndershot ... maximum of the undershot after the pulse in ADC | |
1248 | // pulseRMS ... RMS of the pulse used to calculate the Quality parameters in ADC | |
1249 | // | |
1250 | ||
1251 | TFile file(nameFileQuality,"READ"); | |
1252 | TNtuple *qualityTuple = (TNtuple*)file.Get("TCFquality"); | |
df4cfd77 | 1253 | //gStyle->SetPalette(1); |
d0bd4fcc | 1254 | |
1255 | TH2F *his2D = new TH2F(plotSpec,nameFileQuality,11,-10,1,25,1,100); | |
1256 | char plSpec[100]; | |
56c85970 | 1257 | snprintf(plSpec,100,"%s>>%s",plotSpec,plotSpec); |
d0bd4fcc | 1258 | qualityTuple->Draw(plSpec,cut,pOpt); |
1259 | ||
1260 | gStyle->SetLabelSize(0.03,"X"); | |
1261 | gStyle->SetLabelSize(0.03,"Y"); | |
1262 | gStyle->SetLabelSize(0.03,"Z"); | |
1263 | gStyle->SetLabelOffset(-0.02,"X"); | |
1264 | gStyle->SetLabelOffset(-0.01,"Y"); | |
1265 | gStyle->SetLabelOffset(-0.03,"Z"); | |
1266 | ||
d0bd4fcc | 1267 | his2D->GetXaxis()->SetTitle("max. undershot [ADC]"); |
1268 | his2D->GetYaxis()->SetTitle("width Reduction [%]"); | |
1269 | ||
1270 | his2D->DrawCopy(pOpt); | |
8ad79793 | 1271 | |
1272 | gPad->SetPhi(0.1);gPad->SetTheta(90); | |
d0bd4fcc | 1273 | |
1274 | his2D->~TH2F(); | |
eb7e0771 | 1275 | |
1276 | } | |
1277 | ||
eb7e0771 | 1278 | //_____________________________________________________________________________ |
1279 | Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) { | |
1280 | // | |
1281 | // function to fit one pulse and to calculate the according pole-zero parameters | |
1282 | // | |
1283 | ||
1284 | // initialize TMinuit with a maximum of 8 params | |
7409c8db | 1285 | TMinuit *minuitFit = new TMinuit(8); |
1286 | minuitFit->mncler(); // Reset Minuit's list of paramters | |
1287 | minuitFit->SetPrintLevel(-1); // No Printout | |
1288 | minuitFit->SetFCN(AliTPCCalibTCF::FitFcn); // To set the address of the | |
eb7e0771 | 1289 | // minimization function |
7409c8db | 1290 | minuitFit->SetObjectFit(dataTuple); |
eb7e0771 | 1291 | |
1292 | Double_t arglist[10]; | |
1293 | Int_t ierflg = 0; | |
1294 | ||
1295 | arglist[0] = 1; | |
7409c8db | 1296 | minuitFit->mnexcm("SET ERR", arglist ,1,ierflg); |
eb7e0771 | 1297 | |
1298 | // Set standard starting values and step sizes for each parameter | |
1299 | // upper and lower limit (in a reasonable range) are set to improve | |
1300 | // the stability of TMinuit | |
1301 | static Double_t vstart[8] = {125, 4.0, 0.3, 0.5, 5.5, 100, 1, 2.24}; | |
1302 | static Double_t step[8] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1}; | |
1303 | static Double_t min[8] = {100, 3., 0.1, 0.2, 3., 60., 0., 2.0}; | |
1304 | static Double_t max[8] = {200, 20., 5., 3., 30., 300., 20., 2.5}; | |
1305 | ||
7409c8db | 1306 | minuitFit->mnparm(0, "A1", vstart[0], step[0], min[0], max[0], ierflg); |
1307 | minuitFit->mnparm(1, "A2", vstart[1], step[1], min[1], max[1], ierflg); | |
1308 | minuitFit->mnparm(2, "A3", vstart[2], step[2], min[2], max[2], ierflg); | |
1309 | minuitFit->mnparm(3, "T1", vstart[3], step[3], min[3], max[3], ierflg); | |
1310 | minuitFit->mnparm(4, "T2", vstart[4], step[4], min[4], max[4], ierflg); | |
1311 | minuitFit->mnparm(5, "T3", vstart[5], step[5], min[5], max[5], ierflg); | |
1312 | minuitFit->mnparm(6, "T0", vstart[6], step[6], min[6], max[6], ierflg); | |
1313 | minuitFit->mnparm(7, "TTP", vstart[7], step[7], min[7], max[7],ierflg); | |
1314 | minuitFit->FixParameter(7); // 2.24 ... out of pulserRun Fit (->IRF) | |
eb7e0771 | 1315 | |
1316 | // Now ready for minimization step | |
1317 | arglist[0] = 2000; // max num of iterations | |
1318 | arglist[1] = 0.1; // tolerance | |
1319 | ||
7409c8db | 1320 | minuitFit->mnexcm("MIGRAD", arglist ,2,ierflg); |
eb7e0771 | 1321 | |
1322 | Double_t p1 = 0.0 ; | |
7409c8db | 1323 | minuitFit->mnexcm("SET NOW", &p1 , 0, ierflg) ; // No Warnings |
eb7e0771 | 1324 | |
1325 | if (ierflg == 4) { // Fit failed | |
1326 | for (Int_t i=0;i<3;i++) { | |
1327 | coefP[i] = 0; | |
1328 | coefZ[i] = 0; | |
1329 | } | |
7409c8db | 1330 | minuitFit->~TMinuit(); |
eb7e0771 | 1331 | return 0; |
1332 | } else { // Fit successfull | |
1333 | ||
1334 | // Extract parameters from TMinuit | |
1335 | Double_t *fitParam = new Double_t[6]; | |
1336 | for (Int_t i=0;i<6;i++) { | |
1337 | Double_t err = 0; | |
1338 | Double_t val = 0; | |
7409c8db | 1339 | minuitFit->GetParameter(i,val,err); |
eb7e0771 | 1340 | fitParam[i] = val; |
1341 | } | |
1342 | ||
1343 | // calculates the first 2 sets (4 PZ values) out of the fitted parameters | |
1344 | Double_t *valuePZ = ExtractPZValues(fitParam); | |
1345 | ||
1346 | // TCF coefficients which are used for the equalisation step (stage) | |
1347 | // ZERO/POLE Filter | |
9c2921ef | 1348 | coefZ[0] = TMath::Exp(-1/valuePZ[2]); |
1349 | coefZ[1] = TMath::Exp(-1/valuePZ[3]); | |
1350 | coefP[0] = TMath::Exp(-1/valuePZ[0]); | |
1351 | coefP[1] = TMath::Exp(-1/valuePZ[1]); | |
eb7e0771 | 1352 | |
1353 | fitParam->~Double_t(); | |
1354 | valuePZ->~Double_t(); | |
7409c8db | 1355 | minuitFit->~TMinuit(); |
eb7e0771 | 1356 | |
1357 | return 1; | |
1358 | ||
1359 | } | |
1360 | ||
1361 | } | |
1362 | ||
1363 | ||
1364 | //____________________________________________________________________________ | |
56c85970 | 1365 | void AliTPCCalibTCF::FitFcn(Int_t &/*nPar*/, Double_t */*grad*/, Double_t &f, Double_t * const par, Int_t /*iflag*/) |
eb7e0771 | 1366 | { |
1367 | // | |
1368 | // Minimization function needed for TMinuit with FitFunction included | |
1369 | // Fit function: Sum of three convolution terms (IRF conv. with Exp.) | |
1370 | // | |
1371 | ||
1372 | // Get Data ... | |
1373 | TNtuple *dataTuple = (TNtuple *) gMinuit->GetObjectFit(); | |
1374 | ||
1375 | //calculate chisquare | |
1376 | Double_t chisq = 0; | |
1377 | Double_t delta = 0; | |
1378 | for (Int_t i=0; i<dataTuple->GetEntries(); i++) { // loop over data points | |
1379 | dataTuple->GetEntry(i); | |
1380 | Float_t *p = dataTuple->GetArgs(); | |
1381 | Double_t t = p[0]; | |
1382 | Double_t signal = p[1]; // Normalized signal | |
1383 | Double_t error = p[2]; | |
1384 | ||
1385 | // definition and evaluation if the IonTail specific fit function | |
1386 | Double_t sigFit = 0; | |
1387 | ||
1388 | Double_t ttp = par[7]; // signal shaper raising time | |
1389 | t=t-par[6]; // time adjustment | |
1390 | ||
1391 | if (t<0) { | |
1392 | sigFit = 0; | |
1393 | } else { | |
9c2921ef | 1394 | Double_t f1 = 1/TMath::Power((4-ttp/par[3]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[3]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[3])/ttp+TMath::Power(t*(4-ttp/par[3])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[3])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[3])/ttp,4)/24))); |
eb7e0771 | 1395 | |
9c2921ef | 1396 | Double_t f2 = 1/TMath::Power((4-ttp/par[4]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[4]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[4])/ttp+TMath::Power(t*(4-ttp/par[4])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[4])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[4])/ttp,4)/24))); |
eb7e0771 | 1397 | |
9c2921ef | 1398 | Double_t f3 = 1/TMath::Power((4-ttp/par[5]),5)*(24*ttp*TMath::Exp(4)*(TMath::Exp(-t/par[5]) - TMath::Exp(-4*t/ttp) * ( 1+t*(4-ttp/par[5])/ttp+TMath::Power(t*(4-ttp/par[5])/ttp,2)/2 + TMath::Power(t*(4-ttp/par[5])/ttp,3)/6 + TMath::Power(t*(4-ttp/par[5])/ttp,4)/24))); |
eb7e0771 | 1399 | |
1400 | sigFit = par[0]*f1 + par[1]*f2 +par[2]*f3; | |
1401 | } | |
1402 | ||
1403 | // chisqu calculation | |
1404 | delta = (signal-sigFit)/error; | |
1405 | chisq += delta*delta; | |
1406 | } | |
1407 | ||
1408 | f = chisq; | |
1409 | ||
1410 | } | |
1411 | ||
1412 | ||
1413 | ||
1414 | //____________________________________________________________________________ | |
1415 | Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { | |
1416 | // | |
1417 | // Calculation of Pole and Zero values out of fit parameters | |
1418 | // | |
1419 | ||
1420 | Double_t vA1, vA2, vA3, vTT1, vTT2, vTT3, vTa, vTb; | |
1421 | vA1 = 0; vA2 = 0; vA3 = 0; | |
1422 | vTT1 = 0; vTT2 = 0; vTT3 = 0; | |
1423 | vTa = 0; vTb = 0; | |
1424 | ||
1425 | // nasty method of sorting the fit parameters to avoid wrong mapping | |
1426 | // to the different stages of the TCF filter | |
1427 | // (e.g. first 2 fit parameters represent the electron signal itself!) | |
1428 | ||
56c85970 | 1429 | if ((param[3]-param[4]) <1e-5 ) {param[3]=param[3]+0.0001;} // if equal |
1430 | if ((param[5]-param[4]) <1e-5 ) {param[5]=param[5]+0.0001;} // if equal | |
eb7e0771 | 1431 | |
8ad79793 | 1432 | if ((param[5]>param[4])&&(param[5]>param[3])) { |
eb7e0771 | 1433 | if (param[4]>=param[3]) { |
1434 | vA1 = param[0]; vA2 = param[1]; vA3 = param[2]; | |
1435 | vTT1 = param[3]; vTT2 = param[4]; vTT3 = param[5]; | |
1436 | } else { | |
1437 | vA1 = param[1]; vA2 = param[0]; vA3 = param[2]; | |
1438 | vTT1 = param[4]; vTT2 = param[3]; vTT3 = param[5]; | |
1439 | } | |
8ad79793 | 1440 | } else if ((param[4]>param[5])&&(param[4]>param[3])) { |
eb7e0771 | 1441 | if (param[5]>=param[3]) { |
1442 | vA1 = param[0]; vA2 = param[2]; vA3 = param[1]; | |
1443 | vTT1 = param[3]; vTT2 = param[5]; vTT3 = param[4]; | |
1444 | } else { | |
1445 | vA1 = param[2]; vA2 = param[0]; vA3 = param[1]; | |
1446 | vTT1 = param[5]; vTT2 = param[3]; vTT3 = param[4]; | |
1447 | } | |
8ad79793 | 1448 | } else if ((param[3]>param[4])&&(param[3]>param[5])) { |
eb7e0771 | 1449 | if (param[5]>=param[4]) { |
1450 | vA1 = param[1]; vA2 = param[2]; vA3 = param[0]; | |
1451 | vTT1 = param[4]; vTT2 = param[5]; vTT3 = param[3]; | |
1452 | } else { | |
1453 | vA1 = param[2]; vA2 = param[1]; vA3 = param[0]; | |
1454 | vTT1 = param[5]; vTT2 = param[4]; vTT3 = param[3]; | |
1455 | } | |
1456 | } | |
1457 | ||
1458 | ||
1459 | // Transformation of fit parameters into PZ values (needed by TCF) | |
1460 | Double_t beq = (vA1/vTT2+vA1/vTT3+vA2/vTT1+vA2/vTT3+vA3/vTT1+vA3/vTT2)/(vA1+vA2+vA3); | |
1461 | Double_t ceq = (vA1/(vTT2*vTT3)+vA2/(vTT1*vTT3)+vA3/(vTT1*vTT2))/(vA1+vA2+vA3); | |
1462 | ||
1463 | Double_t s1 = -beq/2-sqrt((beq*beq-4*ceq)/4); | |
1464 | Double_t s2 = -beq/2+sqrt((beq*beq-4*ceq)/4); | |
1465 | ||
1466 | if (vTT2<vTT3) {// not necessary but avoids significant undershots in first PZ | |
1467 | vTa = -1/s1; | |
1468 | vTb = -1/s2; | |
1469 | }else{ | |
1470 | vTa = -1/s2; | |
1471 | vTb = -1/s1; | |
1472 | } | |
1473 | ||
1474 | Double_t *valuePZ = new Double_t[4]; | |
1475 | valuePZ[0]=vTa; | |
1476 | valuePZ[1]=vTb; | |
1477 | valuePZ[2]=vTT2; | |
1478 | valuePZ[3]=vTT3; | |
1479 | ||
1480 | return valuePZ; | |
1481 | ||
1482 | } | |
1483 | ||
1484 | ||
1485 | //____________________________________________________________________________ | |
d0bd4fcc | 1486 | Int_t AliTPCCalibTCF::Equalization(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) { |
eb7e0771 | 1487 | // |
1488 | // calculates the 3rd set of TCF parameters (remaining 2 PZ values) in | |
1489 | // order to restore the original pulse height and adds them to the passed arrays | |
1490 | // | |
1491 | ||
eb7e0771 | 1492 | const Int_t kPulseLength = dataTuple->GetEntries(); |
56c85970 | 1493 | |
1494 | if (kPulseLength<2) { | |
1495 | // prinft("PulseLength does not make sense\n"); | |
1496 | return 0; | |
1497 | } | |
1498 | ||
1499 | Double_t *s0 = new Double_t[kPulseLength]; // original pulse | |
1500 | Double_t *s1 = new Double_t[kPulseLength]; // pulse after 1st PZ filter | |
1501 | Double_t *s2 = new Double_t[kPulseLength]; // pulse after 2nd PZ filter | |
1502 | ||
eb7e0771 | 1503 | for (Int_t ipos=0; ipos<kPulseLength; ipos++) { |
1504 | dataTuple->GetEntry(ipos); | |
1505 | Float_t *p = dataTuple->GetArgs(); | |
1506 | s0[ipos] = p[1]; | |
1507 | } | |
1508 | ||
1509 | // non-discret implementation of the first two TCF stages (recursive formula) | |
1510 | // discrete Altro emulator is not used because of accuracy! | |
1511 | s1[0] = s0[0]; // 1st PZ filter | |
1512 | for(Int_t ipos = 1; ipos < kPulseLength ; ipos++){ | |
1513 | s1[ipos] = s0[ipos] + coefP[0]*s1[ipos-1] - coefZ[0]*s0[ipos-1]; | |
1514 | } | |
1515 | s2[0] = s1[0]; // 2nd PZ filter | |
1516 | for(Int_t ipos = 1; ipos < kPulseLength ; ipos++){ | |
1517 | s2[ipos] = s1[ipos] + coefP[1]*s2[ipos-1] - coefZ[1]*s1[ipos-1]; | |
1518 | } | |
1519 | ||
1520 | // find maximum amplitude and position of original pulse and pulse after | |
1521 | // the first two stages of the TCF | |
1522 | Int_t s0pos = 0, s2pos = 0; | |
1523 | Double_t s0ampl = s0[0], s2ampl = s2[0]; // start values | |
1524 | for(Int_t ipos = 1; ipos < kPulseLength; ipos++){ | |
1525 | if (s0[ipos] > s0ampl){ | |
1526 | s0ampl = s0[ipos]; | |
1527 | s0pos = ipos; // should be pos 11 ... check? | |
1528 | } | |
1529 | if (s2[ipos] > s2ampl){ | |
1530 | s2ampl = s2[ipos]; | |
1531 | s2pos = ipos; | |
1532 | } | |
1533 | } | |
1534 | // calculation of 3rd set ... | |
1535 | if(s0ampl > s2ampl){ | |
1536 | coefZ[2] = 0; | |
1537 | coefP[2] = (s0ampl - s2ampl)/s0[s0pos-1]; | |
1538 | } else if (s0ampl < s2ampl) { | |
1539 | coefP[2] = 0; | |
1540 | coefZ[2] = (s2ampl - s0ampl)/s0[s0pos-1]; | |
1541 | } else { // same height ? will most likely not happen ? | |
d0bd4fcc | 1542 | printf("No equalization because of identical height\n"); |
eb7e0771 | 1543 | coefP[2] = 0; |
1544 | coefZ[2] = 0; | |
1545 | } | |
1546 | ||
56c85970 | 1547 | delete [] s0; |
1548 | delete [] s1; | |
1549 | delete [] s2; | |
df4cfd77 | 1550 | |
d0bd4fcc | 1551 | // if equalization out of range (<0 or >=1) it failed! |
df4cfd77 | 1552 | // if ratio of amplitudes of fittet to original pulse < 0.9 it failed! |
1553 | if (coefP[2]<0 || coefZ[2]<0 || coefP[2]>=1 || coefZ[2]>=1 || TMath::Abs(s2ampl / s0ampl)<0.9) { | |
d0bd4fcc | 1554 | return 0; |
1555 | } else { | |
1556 | return 1; | |
1557 | } | |
df4cfd77 | 1558 | |
eb7e0771 | 1559 | } |
1560 | ||
1561 | ||
1562 | ||
1563 | //____________________________________________________________________________ | |
56c85970 | 1564 | Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F * const hisIn, const char *nameFileTCF, Double_t *coefZ, Double_t *coefP) { |
eb7e0771 | 1565 | // |
1566 | // This function searches for the correct TCF parameters to the given | |
1567 | // histogram 'hisIn' within the file 'nameFileTCF' | |
1568 | // If no parameters for this pad (padinfo within the histogram!) where found | |
1569 | // the function returns 0 | |
1570 | ||
1571 | // Int_t numPulse = (Int_t)hisIn->GetBinContent(1); // number of pulses | |
1572 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
1573 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
1574 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
1575 | Int_t nPulse = 0; | |
1576 | ||
1577 | //-- searching for calculated TCF parameters for this pad/sector | |
1578 | TFile fileTCF(nameFileTCF,"READ"); | |
1579 | TNtuple *paramTuple = (TNtuple*)fileTCF.Get("TCFparam"); | |
1580 | ||
1581 | // create selection criteria to find the correct TCF params | |
1582 | char sel[100]; | |
1583 | if ( paramTuple->GetEntries("row==-1&&pad==-1") ) { | |
1584 | // parameters per SECTOR | |
56c85970 | 1585 | snprintf(sel,100,"sec==%d&&row==-1&&pad==-1",sector); |
eb7e0771 | 1586 | } else { |
1587 | // parameters per PAD | |
56c85970 | 1588 | snprintf(sel,100,"sec==%d&&row==%d&&pad==%d",sector,row,pad); |
eb7e0771 | 1589 | } |
1590 | ||
1591 | // list should contain just ONE entry! ... otherwise there is a mistake! | |
1592 | Long64_t entry = paramTuple->Draw(">>list",sel,"entrylist"); | |
1593 | TEntryList *list = (TEntryList*)gDirectory->Get("list"); | |
1594 | ||
1595 | if (entry) { // TCF set was found for this pad | |
1596 | Long64_t pos = list->GetEntry(0); | |
1597 | paramTuple->GetEntry(pos); // get specific TCF parameters | |
1598 | Float_t *p = paramTuple->GetArgs(); | |
1599 | // check ... | |
56c85970 | 1600 | if((sector-p[0])<1e-5) {printf("sector ok ... "); } |
1601 | if((row-p[1])<1e-5) {printf("row ok ... "); } | |
1602 | if((pad-p[2])<1e-5) {printf("pad ok ... \n"); } | |
eb7e0771 | 1603 | |
1604 | // number of averaged pulses used to produce TCF params | |
1605 | nPulse = (Int_t)p[3]; | |
1606 | // TCF parameters | |
1607 | coefZ[0] = p[4]; coefP[0] = p[7]; | |
1608 | coefZ[1] = p[5]; coefP[1] = p[8]; | |
1609 | coefZ[2] = p[6]; coefP[2] = p[9]; | |
1610 | ||
1611 | } else { // no specific TCF parameters found for this pad | |
1612 | ||
df4cfd77 | 1613 | printf(" no specific TCF paramaters found for pad in ...\n"); |
1614 | printf(" Sector %d | Row %d | Pad %d |\n", sector, row, pad); | |
eb7e0771 | 1615 | nPulse = 0; |
1616 | coefZ[0] = 0; coefP[0] = 0; | |
1617 | coefZ[1] = 0; coefP[1] = 0; | |
1618 | coefZ[2] = 0; coefP[2] = 0; | |
1619 | ||
1620 | } | |
1621 | ||
1622 | fileTCF.Close(); | |
1623 | ||
1624 | return nPulse; // number of averaged pulses for producing the TCF params | |
1625 | ||
1626 | } | |
1627 | ||
1628 | ||
1629 | //____________________________________________________________________________ | |
1630 | Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t *coefP, Int_t plotFlag) { | |
1631 | // | |
1632 | // This function evaluates the quality parameters of the given TCF parameters | |
1633 | // tested on the passed pulse (hisIn) | |
1634 | // The quality parameters are stored in an array. They are ... | |
1635 | // height deviation [ADC] | |
1636 | // area reduction [percent] | |
1637 | // width reduction [percent] | |
1638 | // mean undershot [ADC] | |
1639 | // maximum of undershot after pulse [ADC] | |
1640 | // Pulse RMS [ADC] | |
1641 | ||
1642 | // perform ALTRO emulator | |
1643 | TNtuple *pulseTuple = ApplyTCFilter(hisIn, coefZ, coefP, plotFlag); | |
1644 | ||
1645 | printf("calculate quality val. for pulse in ... "); | |
1646 | printf(" Sector %d | Row %d | Pad %d |\n", (Int_t)hisIn->GetBinContent(2), (Int_t)hisIn->GetBinContent(3), (Int_t)hisIn->GetBinContent(4)); | |
1647 | ||
1648 | // Reasonable limit for the calculation of the quality values | |
1649 | Int_t binLimit = 80; | |
1650 | ||
1651 | // ============== Variable preparation | |
1652 | ||
1653 | // -- height difference in percent of orginal pulse | |
1654 | Double_t maxSig = pulseTuple->GetMaximum("sig"); | |
1655 | Double_t maxSigTCF = pulseTuple->GetMaximum("sigAfterTCF"); | |
1656 | // -- area reduction (above zero!) | |
1657 | Double_t area = 0; | |
1658 | Double_t areaTCF = 0; | |
1659 | // -- width reduction at certain ADC treshold | |
1660 | // TODO: set treshold at ZS treshold? (3 sigmas of noise?) | |
1661 | Int_t threshold = 3; // treshold in percent | |
1662 | Int_t threshADC = (Int_t)(maxSig/100*threshold); | |
1663 | Int_t startOfPulse = 0; Int_t startOfPulseTCF = 0; | |
1664 | Int_t posOfStart = 0; Int_t posOfStartTCF = 0; | |
1665 | Int_t widthFound = 0; Int_t widthFoundTCF = 0; | |
1666 | Int_t width = 0; Int_t widthTCF = 0; | |
1667 | // -- Calcluation of Undershot (mean of negavive signal after the first | |
1668 | // undershot) | |
1669 | Double_t undershotTCF = 0; | |
1670 | Double_t undershotStart = 0; | |
1671 | // -- Calcluation of Undershot (Sum of negative signal after the pulse) | |
1672 | Double_t maxUndershot = 0; | |
1673 | ||
1674 | ||
1675 | // === loop over timebins to calculate quality parameters | |
1676 | for (Int_t i=0; i<binLimit; i++) { | |
1677 | ||
1678 | // Read signal values | |
1679 | pulseTuple->GetEntry(i); | |
1680 | Float_t *p = pulseTuple->GetArgs(); | |
1681 | Double_t sig = p[1]; | |
1682 | Double_t sigTCF = p[2]; | |
1683 | ||
1684 | // calculation of area (above zero) | |
1685 | if (sig>0) {area += sig; } | |
1686 | if (sigTCF>0) {areaTCF += sigTCF; } | |
1687 | ||
1688 | ||
1689 | // Search for width at certain ADC treshold | |
1690 | // -- original signal | |
1691 | if (widthFound == 0) { | |
1692 | if( (sig > threshADC) && (startOfPulse == 0) ){ | |
1693 | startOfPulse = 1; | |
1694 | posOfStart = i; | |
1695 | } | |
d0bd4fcc | 1696 | if( (sig <= threshADC) && (startOfPulse == 1) ){ |
eb7e0771 | 1697 | widthFound = 1; |
1698 | width = i - posOfStart + 1; | |
1699 | } | |
1700 | } | |
1701 | // -- signal after TCF | |
1702 | if (widthFoundTCF == 0) { | |
1703 | if( (sigTCF > threshADC) && (startOfPulseTCF == 0) ){ | |
1704 | startOfPulseTCF = 1; | |
1705 | posOfStartTCF = i; | |
1706 | } | |
d0bd4fcc | 1707 | if( (sigTCF <= threshADC) && (startOfPulseTCF == 1) ){ |
eb7e0771 | 1708 | widthFoundTCF = 1; |
1709 | widthTCF = i -posOfStartTCF + 1; | |
1710 | } | |
1711 | ||
1712 | } | |
1713 | ||
1714 | // finds undershot start | |
1715 | if ( (widthFoundTCF==1) && (sigTCF<0) ) { | |
1716 | undershotStart = 1; | |
1717 | } | |
1718 | ||
1719 | // Calculation of undershot sum (after pulse) | |
1720 | if ( widthFoundTCF==1 ) { | |
1721 | undershotTCF += sigTCF; | |
1722 | } | |
1723 | ||
1724 | // Search for maximal undershot (is equal to minimum after the pulse) | |
56c85970 | 1725 | if ( ((undershotStart-1)<1e-7)&&(i<(posOfStartTCF+widthTCF+20)) ) { |
eb7e0771 | 1726 | if (maxUndershot>sigTCF) { maxUndershot = sigTCF; } |
1727 | } | |
1728 | ||
1729 | } | |
1730 | ||
1731 | // == Calculation of Quality parameters | |
1732 | ||
1733 | // -- height difference in ADC | |
1734 | Double_t heightDev = maxSigTCF-maxSig; | |
1735 | ||
1736 | // Area reduction of the pulse in percent | |
1737 | Double_t areaReduct = 100-areaTCF/area*100; | |
1738 | ||
1739 | // Width reduction in percent | |
1740 | Double_t widthReduct = 0; | |
1741 | if ((widthFound==1)&&(widthFoundTCF==1)) { // in case of not too big IonTail | |
1742 | widthReduct = 100-(Double_t)widthTCF/(Double_t)width*100; | |
1743 | if (widthReduct<0) { widthReduct = 0;} | |
1744 | } | |
1745 | ||
1746 | // Undershot - mean of neg.signals after pulse | |
1747 | Double_t length = 1; | |
1748 | if (binLimit-widthTCF-posOfStartTCF) { length = (binLimit-widthTCF-posOfStartTCF);} | |
1749 | Double_t undershot = undershotTCF/length; | |
1750 | ||
1751 | ||
1752 | // calculation of pulse RMS with timebins before and at the end of the pulse | |
1753 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",100,-50,50); | |
1754 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
1755 | // before the pulse | |
1756 | tempRMSHis->Fill(hisIn->GetBinContent(ipos+5)); | |
55f06b51 | 1757 | // at the end |
eb7e0771 | 1758 | tempRMSHis->Fill(hisIn->GetBinContent(hisIn->GetNbinsX()-ipos)); |
1759 | } | |
1760 | Double_t pulseRMS = tempRMSHis->GetRMS(); | |
1761 | tempRMSHis->~TH1I(); | |
1762 | ||
1763 | if (plotFlag) { | |
1764 | // == Output | |
1765 | printf("height deviation [ADC]:\t\t\t %3.1f\n", heightDev); | |
1766 | printf("area reduction [percent]:\t\t %3.1f\n", areaReduct); | |
1767 | printf("width reduction [percent]:\t\t %3.1f\n", widthReduct); | |
1768 | printf("mean undershot [ADC]:\t\t\t %3.1f\n", undershot); | |
1769 | printf("maximum of undershot after pulse [ADC]: %3.1f\n", maxUndershot); | |
55f06b51 | 1770 | printf("RMS of the original (or summed) pulse [ADC]: \t %3.2f\n\n", pulseRMS); |
eb7e0771 | 1771 | |
1772 | } | |
1773 | ||
1774 | Double_t *qualityParam = new Double_t[6]; | |
1775 | qualityParam[0] = heightDev; | |
1776 | qualityParam[1] = areaReduct; | |
1777 | qualityParam[2] = widthReduct; | |
1778 | qualityParam[3] = undershot; | |
1779 | qualityParam[4] = maxUndershot; | |
1780 | qualityParam[5] = pulseRMS; | |
1781 | ||
1782 | pulseTuple->~TNtuple(); | |
1783 | ||
1784 | return qualityParam; | |
1785 | } | |
1786 | ||
1787 | ||
1788 | //____________________________________________________________________________ | |
56c85970 | 1789 | TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F * const hisIn, Double_t * const coefZ, Double_t * const coefP, Int_t plotFlag) { |
eb7e0771 | 1790 | // |
1791 | // Applies the given TCF parameters on the given pulse via the ALTRO emulator | |
1792 | // class (discret values) and stores both pulses into a returned TNtuple | |
1793 | // | |
1794 | ||
1795 | Int_t nbins = hisIn->GetNbinsX() -4; | |
1796 | // -1 because the first four timebins usually contain pad specific informations | |
1797 | Int_t nPulse = (Int_t)hisIn->GetBinContent(1); // Number of summed pulses | |
1798 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
1799 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
1800 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
1801 | ||
1802 | // redirect histogram values to arrays (discrete for altro emulator) | |
1803 | Double_t *signalIn = new Double_t[nbins]; | |
1804 | Double_t *signalOut = new Double_t[nbins]; | |
1805 | short *signalInD = new short[nbins]; | |
1806 | short *signalOutD = new short[nbins]; | |
1807 | for (Int_t ipos=0;ipos<nbins;ipos++) { | |
1808 | Double_t signal = hisIn->GetBinContent(ipos+5); // summed signal | |
1809 | signalIn[ipos]=signal/nPulse; // mean signal | |
1810 | signalInD[ipos]=(short)(TMath::Nint(signalIn[ipos])); //discrete mean signal | |
1811 | signalOutD[ipos]=signalInD[ipos]; // will be overwritten by AltroEmulator | |
1812 | } | |
1813 | ||
1814 | // transform TCF parameters into ALTRO readable format (Integer) | |
6fefa5ce | 1815 | Int_t valK[3]; |
1816 | Int_t valL[3]; | |
eb7e0771 | 1817 | for (Int_t i=0; i<3; i++) { |
df4cfd77 | 1818 | valK[i] = (Int_t)(coefP[i]*(TMath::Power(2,16)-1)); |
1819 | valL[i] = (Int_t)(coefZ[i]*(TMath::Power(2,16)-1)); | |
eb7e0771 | 1820 | } |
1821 | ||
1822 | // discret ALTRO EMULATOR ____________________________ | |
1823 | AliTPCAltroEmulator *altro = new AliTPCAltroEmulator(nbins, signalOutD); | |
1824 | altro->ConfigAltro(0,1,0,0,0,0); // perform just the TailCancelation | |
df4cfd77 | 1825 | altro->ConfigTailCancellationFilter(valK[0],valK[1],valK[2],valL[0],valL[1],valL[2]); |
eb7e0771 | 1826 | altro->RunEmulation(); |
1827 | delete altro; | |
1828 | ||
1829 | // non-discret implementation of the (recursive formula) | |
1830 | // discrete Altro emulator is not used because of accuracy! | |
1831 | Double_t *s1 = new Double_t[1000]; // pulse after 1st PZ filter | |
1832 | Double_t *s2 = new Double_t[1000]; // pulse after 2nd PZ filter | |
1833 | s1[0] = signalIn[0]; // 1st PZ filter | |
1834 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1835 | s1[ipos] = signalIn[ipos] + coefP[0]*s1[ipos-1] - coefZ[0]*signalIn[ipos-1]; | |
1836 | } | |
1837 | s2[0] = s1[0]; // 2nd PZ filter | |
1838 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1839 | s2[ipos] = s1[ipos] + coefP[1]*s2[ipos-1] - coefZ[1]*s1[ipos-1]; | |
1840 | } | |
1841 | signalOut[0] = s2[0]; // 3rd PZ filter | |
1842 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1843 | signalOut[ipos] = s2[ipos] + coefP[2]*signalOut[ipos-1] - coefZ[2]*s2[ipos-1]; | |
1844 | } | |
1845 | s1->~Double_t(); | |
1846 | s2->~Double_t(); | |
1847 | ||
1848 | // writing pulses to tuple | |
1849 | TNtuple *pulseTuple = new TNtuple("ntupleTCF","PulseTCF","timebin:sig:sigAfterTCF:sigND:sigNDAfterTCF"); | |
1850 | for (Int_t ipos=0;ipos<nbins;ipos++) { | |
1851 | pulseTuple->Fill(ipos,signalInD[ipos],signalOutD[ipos],signalIn[ipos],signalOut[ipos]); | |
1852 | } | |
1853 | ||
1854 | if (plotFlag) { | |
1855 | char hname[100]; | |
56c85970 | 1856 | snprintf(hname,100,"sec%drow%dpad%d",sector,row,pad); |
eb7e0771 | 1857 | new TCanvas(hname,hname,600,400); |
1858 | //just plotting non-discret pulses | they look pretties in case of mean sig ;-) | |
1859 | pulseTuple->Draw("sigND:timebin","","L"); | |
1860 | // pulseTuple->Draw("sig:timebin","","Lsame"); | |
1861 | pulseTuple->SetLineColor(3); | |
1862 | pulseTuple->Draw("sigNDAfterTCF:timebin","","Lsame"); | |
1863 | // pulseTuple->Draw("sigAfterTCF:timebin","","Lsame"); | |
1864 | } | |
1865 | ||
4ce766eb | 1866 | delete [] signalIn; |
1867 | delete [] signalOut; | |
56c85970 | 1868 | delete [] signalInD; |
1869 | delete [] signalOutD; | |
1870 | ||
eb7e0771 | 1871 | return pulseTuple; |
1872 | ||
1873 | } | |
1874 | ||
1875 | ||
eb7e0771 | 1876 | //____________________________________________________________________________ |
1877 | void AliTPCCalibTCF::PrintPulseThresholds() { | |
1878 | // | |
1879 | // Prints the pulse threshold settings | |
1880 | // | |
1881 | ||
1882 | printf(" %4.0d [ADC] ... expected Gate fluctuation length \n", fGateWidth); | |
1883 | printf(" %4.0d [ADC] ... expected usefull signal length \n", fSample); | |
1884 | printf(" %4.0d [ADC] ... needed pulselength for TC characterisation \n", fPulseLength); | |
1885 | printf(" %4.0d [ADC] ... lower pulse height limit \n", fLowPulseLim); | |
1886 | printf(" %4.0d [ADC] ... upper pulse height limit \n", fUpPulseLim); | |
1887 | printf(" %4.1f [ADC] ... maximal pulse RMS \n", fRMSLim); | |
df4cfd77 | 1888 | printf(" %4.1f [ADC] ... pulse/tail integral ratio \n", fRatioIntLim); |
eb7e0771 | 1889 | |
1890 | } | |
d0bd4fcc | 1891 | |
1892 | ||
1893 | //____________________________________________________________________________ | |
df4cfd77 | 1894 | void AliTPCCalibTCF::MergeHistoPerFile(const char *fileNameIn, const char *fileNameSum, Int_t mode) |
d0bd4fcc | 1895 | { |
df4cfd77 | 1896 | // Gets histograms from fileNameIn and adds contents to fileSum |
1897 | // | |
1898 | // If fileSum doesn't exist, fileSum is created | |
1899 | // mode = 0, just ONE BIG FILE ('fileSum') will be used | |
1900 | // mode = 1, one file per sector ('fileSum-Sec#.root') will be used | |
1901 | // mode=1 is much faster, but the additional function 'MergeToOneFile' has to be used in order to | |
1902 | // get one big and complete collection file again ... | |
d0bd4fcc | 1903 | // |
df4cfd77 | 1904 | // !Make sure not to add the same file more than once! |
d0bd4fcc | 1905 | |
1906 | TFile fileIn(fileNameIn,"READ"); | |
1907 | TH1F *hisIn; | |
1908 | TKey *key; | |
1909 | TIter next(fileIn.GetListOfKeys()); | |
56c85970 | 1910 | // opens a file, although, it might not be uses (see "mode") |
1911 | TFile *fileOut = new TFile(fileNameSum,"UPDATE"); | |
d0bd4fcc | 1912 | //fileOut.cd(); |
1913 | ||
1914 | Int_t nHist=fileIn.GetNkeys(); | |
1915 | Int_t iHist=0; | |
df4cfd77 | 1916 | |
1917 | Int_t secPrev = -1; | |
1918 | char fileNameSumSec[100]; | |
1919 | ||
56c85970 | 1920 | |
d0bd4fcc | 1921 | while((key=(TKey*)next())) { |
1922 | const char *hisName = key->GetName(); | |
1923 | ||
752b0cc7 | 1924 | TString name(key->GetName()); |
1925 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
1926 | ||
d0bd4fcc | 1927 | hisIn=(TH1F*)fileIn.Get(hisName); |
338e0dd9 | 1928 | |
d0bd4fcc | 1929 | Int_t numPulse=(Int_t)hisIn->GetBinContent(1); |
df4cfd77 | 1930 | Int_t sec=(Int_t)hisIn->GetBinContent(2); |
d0bd4fcc | 1931 | Int_t pulseLength= hisIn->GetNbinsX()-4; |
1932 | ||
df4cfd77 | 1933 | // in case of mode 1, store histos in files per sector |
1934 | if (sec!=secPrev && mode != 0) { | |
1935 | if (secPrev>0) { // closing old file | |
1936 | fileOut->Close(); | |
1937 | } | |
1938 | // opening new file | |
56c85970 | 1939 | snprintf(fileNameSumSec,100,"%s-Sec%d.root",fileNameSum,sec); |
df4cfd77 | 1940 | fileOut = new TFile(fileNameSumSec,"UPDATE"); |
1941 | secPrev = sec; | |
1942 | } | |
d0bd4fcc | 1943 | |
df4cfd77 | 1944 | // search for existing histogram |
1945 | TH1F *his=(TH1F*)fileOut->Get(hisName); | |
1946 | if (iHist%100==0) { | |
1947 | printf("Histogram %d / %d, %s, Action: ",iHist,nHist,hisName); | |
1948 | if (!his) { | |
1949 | printf("NEW\n"); | |
1950 | } else { | |
1951 | printf("ADD\n"); | |
1952 | } | |
1953 | } | |
1954 | iHist++; | |
1955 | ||
d0bd4fcc | 1956 | if (!his) { |
d0bd4fcc | 1957 | his=hisIn; |
1958 | his->Write(hisName); | |
1959 | } else { | |
d0bd4fcc | 1960 | his->AddBinContent(1,numPulse); |
1961 | for (Int_t ii=5; ii<pulseLength+5; ii++) { | |
1962 | his->AddBinContent(ii,hisIn->GetBinContent(ii)); | |
1963 | } | |
1964 | his->Write(hisName,TObject::kOverwrite); | |
1965 | } | |
1966 | } | |
df4cfd77 | 1967 | |
d0bd4fcc | 1968 | printf("closing files (may take a while)...\n"); |
df4cfd77 | 1969 | fileOut->Close(); |
1970 | ||
1971 | ||
d0bd4fcc | 1972 | fileIn.Close(); |
1973 | printf("...DONE\n\n"); | |
1974 | } | |
1975 | ||
df4cfd77 | 1976 | |
1977 | //____________________________________________________________________________ | |
1978 | void AliTPCCalibTCF::MergeToOneFile(const char *nameFileSum) { | |
1979 | ||
1980 | // Merges all Sec-files together ... | |
1981 | // this is an additional functionality for the function MergeHistsPerFile | |
1982 | // if for example mode=1 | |
1983 | ||
1984 | TH1F *hisIn; | |
1985 | TKey *key; | |
1986 | ||
1987 | // just delete the file entries ... | |
7409c8db | 1988 | TFile fileSumD(nameFileSum,"RECREATE"); |
1989 | fileSumD.Close(); | |
df4cfd77 | 1990 | |
1991 | char nameFileSumSec[100]; | |
1992 | ||
1993 | for (Int_t sec=0; sec<72; sec++) { // loop over all possible filenames | |
1994 | ||
56c85970 | 1995 | snprintf(nameFileSumSec,100,"%s-Sec%d.root",nameFileSum,sec); |
df4cfd77 | 1996 | TFile *fileSumSec = new TFile(nameFileSumSec,"READ"); |
1997 | ||
1998 | Int_t nHist=fileSumSec->GetNkeys(); | |
1999 | Int_t iHist=0; | |
2000 | ||
2001 | if (nHist) { // file found \ NKeys not empty | |
2002 | ||
2003 | TFile fileSum(nameFileSum,"UPDATE"); | |
2004 | fileSum.cd(); | |
2005 | ||
2006 | printf("Sector file %s found\n",nameFileSumSec); | |
2007 | TIter next(fileSumSec->GetListOfKeys()); | |
71a41e6c | 2008 | while( (key=(TKey*)next()) ) { |
df4cfd77 | 2009 | const char *hisName = key->GetName(); |
338e0dd9 | 2010 | TString name(hisName); |
2011 | if (name.Contains("ddl") ) continue; // ignore the 2d histogramms per ddl | |
df4cfd77 | 2012 | hisIn=(TH1F*)fileSumSec->Get(hisName); |
2013 | ||
338e0dd9 | 2014 | |
df4cfd77 | 2015 | if (iHist%100==0) { |
2016 | printf("found histogram %d / %d, %s\n",iHist,nHist,hisName); | |
2017 | } | |
2018 | iHist++; | |
2019 | ||
2020 | // TH1F *his = (TH1F*)hisIn->Clone(hisName); | |
2021 | hisIn->Write(hisName); | |
2022 | ||
2023 | } | |
2024 | printf("Saving histograms from sector %d (may take a while) ...",sec); | |
2025 | fileSum.Close(); | |
2026 | ||
2027 | } | |
2028 | fileSumSec->Close(); | |
2029 | } | |
2030 | printf("...DONE\n\n"); | |
2031 | } | |
2032 | ||
2033 | ||
2034 | //____________________________________________________________________________ | |
2035 | Int_t AliTPCCalibTCF::DumpTCFparamToFilePerPad(const char *nameFileTCFPerPad,const char *nameFileTCFPerSec, const char *nameMappingFile) { | |
2036 | // | |
2037 | // Writes TCF parameters per PAD to .data file | |
2038 | // | |
2039 | // from now on: "roc" refers to the offline sector numbering | |
2040 | // "sector" refers to the 18 sectors per side | |
2041 | // | |
2042 | // Gets TCF parameters of single pads from nameFileTCFPerPad and writes them to | |
2043 | // the file 'tpcTCFparamPAD.data' | |
2044 | // | |
2045 | // If there are parameters for a pad missing, then the parameters of the roc, | |
2046 | // in which the pad is located, are used as the pad parameters. The parameters for | |
2047 | // the roc are retreived from nameFileTCFPerSec. If there are parameters for | |
2048 | // a roc missing, then the parameters are set to -1. | |
2049 | ||
56c85970 | 2050 | Float_t k0 = -1, k1 = -1, k2 = -1, l0 = -1, l1 = -1, l2 = -1; |
df4cfd77 | 2051 | Int_t roc, row, pad, side, sector, rcu, hwAddr; |
2052 | Int_t entryNum = 0; | |
2053 | Int_t checksum = 0; | |
2054 | Int_t tpcPadNum = 557568; | |
2055 | Int_t validFlag = 1; // 1 if parameters for pad exist, 0 if they are only inherited from the roc | |
2056 | ||
2057 | Bool_t *entryID = new Bool_t[7200000]; // helping vector | |
2058 | for (Int_t ii = 0; ii<7200000; ii++) { | |
2059 | entryID[ii]=0; | |
2060 | } | |
2061 | ||
2062 | // get file/tuple with parameters per pad | |
2063 | TFile fileTCFparam(nameFileTCFPerPad); | |
2064 | TNtuple *paramTuple = (TNtuple*)fileTCFparam.Get("TCFparam"); | |
2065 | ||
2066 | // get mapping file | |
2067 | // usual location of mapping file: $ALICE_ROOT/TPC/Calib/tpcMapping.root | |
2068 | TFile *fileMapping = new TFile(nameMappingFile, "read"); | |
2069 | AliTPCmapper *mapping = (AliTPCmapper*) fileMapping->Get("tpcMapping"); | |
2070 | delete fileMapping; | |
2071 | ||
2072 | if (mapping == 0) { | |
2073 | printf("Failed to get mapping object from %s. ...\n", nameMappingFile); | |
2074 | return -1; | |
2075 | } else { | |
2076 | printf("Got mapping object from %s\n", nameMappingFile); | |
2077 | } | |
2078 | ||
2079 | // creating outputfile | |
2080 | ofstream fileOut; | |
2081 | char nameFileOut[255]; | |
56c85970 | 2082 | snprintf(nameFileOut,255,"tpcTCFparamPAD.data"); |
df4cfd77 | 2083 | fileOut.open(nameFileOut); |
2084 | // following not used: | |
2085 | // char headerLine[255]; | |
56c85970 | 2086 | // snprintf(headerLine,255,"15\tside\tsector\tRCU\tHWadr\tk0\tk1\tk2\tl0\tl1\tl2\tValidFlag"); |
df4cfd77 | 2087 | // fileOut << headerLine << std::endl; |
2088 | fileOut << "15" << std::endl; | |
2089 | ||
2090 | // loop over nameFileTCFPerPad, write parameters into outputfile | |
2091 | // NOTE: NO SPECIFIC ORDER !!! | |
2092 | printf("\nstart assigning parameters to pad...\n"); | |
2093 | for (Int_t iParam = 0; iParam < paramTuple->GetEntries(); iParam++) { | |
2094 | paramTuple->GetEntry(iParam); | |
2095 | Float_t *paramArgs = paramTuple->GetArgs(); | |
2096 | roc = Int_t(paramArgs[0]); | |
2097 | row = Int_t(paramArgs[1]); | |
2098 | pad = Int_t(paramArgs[2]); | |
2099 | side = Int_t(mapping->GetSideFromRoc(roc)); | |
2100 | sector = Int_t(mapping->GetSectorFromRoc(roc)); | |
2101 | rcu = Int_t(mapping->GetRcu(roc,row,pad)); | |
2102 | hwAddr = Int_t(mapping->GetHWAddress(roc,row,pad)); | |
56c85970 | 2103 | k0 = TMath::Nint(paramArgs[7] * (TMath::Power(2,16) - 1)); |
2104 | k1 = TMath::Nint(paramArgs[8] * (TMath::Power(2,16) - 1)); | |
2105 | k2 = TMath::Nint(paramArgs[9] * (TMath::Power(2,16) - 1)); | |
2106 | l0 = TMath::Nint(paramArgs[4] * (TMath::Power(2,16) - 1)); | |
2107 | l1 = TMath::Nint(paramArgs[5] * (TMath::Power(2,16) - 1)); | |
2108 | l2 = TMath::Nint(paramArgs[6] * (TMath::Power(2,16) - 1)); | |
df4cfd77 | 2109 | if (entryNum%10000==0) { |
2110 | printf("assigned pad %i / %i\n",entryNum,tpcPadNum); | |
2111 | } | |
2112 | ||
2113 | fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << hwAddr << "\t"; | |
56c85970 | 2114 | fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl; |
df4cfd77 | 2115 | entryID[roc*100000 + row*1000 + pad] = 1; |
2116 | } | |
2117 | ||
2118 | // Wrote all found TCF params per pad into data file | |
2119 | // NOW FILLING UP THE REST WITH THE PARAMETERS FROM THE ROC MEAN | |
2120 | ||
2121 | // get file/tuple with parameters per roc | |
2122 | TFile fileSecTCFparam(nameFileTCFPerSec); | |
2123 | TNtuple *paramTupleSec = (TNtuple*)fileSecTCFparam.Get("TCFparam"); | |
2124 | ||
2125 | // loop over all pads and get/write parameters for pads which don't have | |
2126 | // parameters assigned yet | |
2127 | validFlag = 0; | |
2128 | for (roc = 0; roc<72; roc++) { | |
2129 | side = Int_t(mapping->GetSideFromRoc(roc)); | |
2130 | sector = Int_t(mapping->GetSectorFromRoc(roc)); | |
2131 | for (Int_t iParamSec = 0; iParamSec < paramTupleSec->GetEntries(); iParamSec++) { | |
2132 | paramTupleSec->GetEntry(iParamSec); | |
2133 | Float_t *paramArgsSec = paramTupleSec->GetArgs(); | |
56c85970 | 2134 | if ((paramArgsSec[0]-roc)<1e-7) { // if roc is found |
2135 | k0 = TMath::Nint(paramArgsSec[7] * (TMath::Power(2,16) - 1)); | |
2136 | k1 = TMath::Nint(paramArgsSec[8] * (TMath::Power(2,16) - 1)); | |
2137 | k2 = TMath::Nint(paramArgsSec[9] * (TMath::Power(2,16) - 1)); | |
2138 | l0 = TMath::Nint(paramArgsSec[4] * (TMath::Power(2,16) - 1)); | |
2139 | l1 = TMath::Nint(paramArgsSec[5] * (TMath::Power(2,16) - 1)); | |
2140 | l2 = TMath::Nint(paramArgsSec[6] * (TMath::Power(2,16) - 1)); | |
df4cfd77 | 2141 | break; |
2142 | } else { | |
56c85970 | 2143 | k0 = k1 = k2 = l0 = l1 = l2 = -1; |
df4cfd77 | 2144 | } |
2145 | } | |
2146 | for (row = 0; row<mapping->GetNpadrows(roc); row++) { | |
2147 | for (pad = 0; pad<mapping->GetNpads(roc,row); pad++) { | |
2148 | if (entryID[roc*100000 + row*1000 + pad]==1) { | |
2149 | continue; | |
2150 | } | |
2151 | ||
2152 | entryID[roc*100000 + row*1000 + pad] = 1; | |
2153 | rcu = Int_t(mapping->GetRcu(roc,row,pad)); | |
2154 | hwAddr = Int_t(mapping->GetHWAddress(roc,row,pad)); | |
2155 | if (entryNum%10000==0) { | |
2156 | printf("assigned pad %i / %i\n",entryNum,tpcPadNum); | |
2157 | } | |
2158 | ||
2159 | fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << hwAddr << "\t"; | |
56c85970 | 2160 | fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl; |
df4cfd77 | 2161 | } |
2162 | } | |
2163 | } | |
2164 | ||
2165 | printf("assigned pad %i / %i\ndone assigning\n",entryNum,tpcPadNum); | |
2166 | ||
2167 | // check if correct amount of sets of parameters were written | |
2168 | for (Int_t ii = 0; ii<7200000; ii++) { | |
2169 | checksum += entryID[ii]; | |
2170 | } | |
2171 | if (checksum == tpcPadNum) { | |
2172 | printf("checksum ok, sets of parameters written = %i\n",checksum); | |
2173 | } else { | |
2174 | printf("\nCHECKSUM WRONG, sets of parameters written = %i, should be %i\n\n",checksum,tpcPadNum); | |
2175 | } | |
2176 | ||
2177 | // closing & destroying | |
2178 | fileOut.close(); | |
2179 | fileTCFparam.Close(); | |
2180 | fileSecTCFparam.Close(); | |
56c85970 | 2181 | delete [] entryID; |
df4cfd77 | 2182 | printf("output written to file: %s\n",nameFileOut); |
2183 | return 0; | |
2184 | } | |
2185 | ||
2186 | ||
2187 | ||
2188 | //____________________________________________________________________________ | |
2189 | Int_t AliTPCCalibTCF::DumpTCFparamToFilePerSector(const char *nameFileTCFPerSec, const char *nameMappingFile) { | |
2190 | // | |
2191 | // Writes TCF parameters per SECTOR (=ROC) to .data file | |
2192 | // | |
2193 | // from now on: "roc" refers to the offline sector numbering | |
2194 | // "sector" refers to the 18 sectors per side | |
2195 | // | |
2196 | // Gets TCF parameters of a roc from nameFileTCFPerSec and writes them to | |
2197 | // the file 'tpcTCFparamSector.data' | |
2198 | // | |
2199 | // If there are parameters for a roc missing, then the parameters are set to -1 | |
2200 | ||
56c85970 | 2201 | Float_t k0 = -1, k1 = -1, k2 = -1, l0 = -1, l1 = -1, l2 = -1; |
df4cfd77 | 2202 | Int_t entryNum = 0; |
2203 | Int_t validFlag = 0; // 1 if parameters for roc exist | |
2204 | ||
2205 | // get file/tuple with parameters per roc | |
2206 | TFile fileTCFparam(nameFileTCFPerSec); | |
2207 | TNtuple *paramTupleSec = (TNtuple*)fileTCFparam.Get("TCFparam"); | |
2208 | ||
2209 | ||
2210 | // get mapping file | |
2211 | // usual location of mapping file: $ALICE_ROOT/TPC/Calib/tpcMapping.root | |
2212 | TFile *fileMapping = new TFile(nameMappingFile, "read"); | |
2213 | AliTPCmapper *mapping = (AliTPCmapper*) fileMapping->Get("tpcMapping"); | |
2214 | delete fileMapping; | |
2215 | ||
2216 | if (mapping == 0) { | |
2217 | printf("Failed to get mapping object from %s. ...\n", nameMappingFile); | |
2218 | return -1; | |
2219 | } else { | |
2220 | printf("Got mapping object from %s\n", nameMappingFile); | |
2221 | } | |
2222 | ||
2223 | ||
2224 | // creating outputfile | |
2225 | ||
2226 | ofstream fileOut; | |
2227 | char nameFileOut[255]; | |
56c85970 | 2228 | snprintf(nameFileOut,255,"tpcTCFparamSector.data"); |
df4cfd77 | 2229 | fileOut.open(nameFileOut); |
2230 | // following not used: | |
2231 | // char headerLine[255]; | |
56c85970 | 2232 | // snprintf(headerLine,255,"16\tside\tsector\tRCU\tHWadr\tk0\tk1\tk2\tl0\tl1\tl2\tValidFlag"); |
df4cfd77 | 2233 | // fileOut << headerLine << std::endl; |
2234 | fileOut << "16" << std::endl; | |
2235 | ||
2236 | // loop over all rcu's in the TPC (6 per sector) | |
2237 | printf("\nstart assigning parameters to rcu's...\n"); | |
2238 | for (Int_t side = 0; side<2; side++) { | |
2239 | for (Int_t sector = 0; sector<18; sector++) { | |
2240 | for (Int_t rcu = 0; rcu<6; rcu++) { | |
2241 | ||
2242 | validFlag = 0; | |
2243 | Int_t roc = Int_t(mapping->GetRocFromPatch(side, sector, rcu)); | |
2244 | ||
2245 | // get parameters (through loop search) for rcu from corresponding roc | |
2246 | for (Int_t iParam = 0; iParam < paramTupleSec->GetEntries(); iParam++) { | |
2247 | paramTupleSec->GetEntry(iParam); | |
2248 | Float_t *paramArgs = paramTupleSec->GetArgs(); | |
56c85970 | 2249 | if ((paramArgs[0]-roc)<1e-7) { // if roc is found |
df4cfd77 | 2250 | validFlag = 1; |
56c85970 | 2251 | k0 = TMath::Nint(paramArgs[7] * (TMath::Power(2,16) - 1)); |
2252 | k1 = TMath::Nint(paramArgs[8] * (TMath::Power(2,16) - 1)); | |
2253 | k2 = TMath::Nint(paramArgs[9] * (TMath::Power(2,16) - 1)); | |
2254 | l0 = TMath::Nint(paramArgs[4] * (TMath::Power(2,16) - 1)); | |
2255 | l1 = TMath::Nint(paramArgs[5] * (TMath::Power(2,16) - 1)); | |
2256 | l2 = TMath::Nint(paramArgs[6] * (TMath::Power(2,16) - 1)); | |
df4cfd77 | 2257 | break; |
2258 | } | |
2259 | } | |
2260 | if (!validFlag) { // No TCF parameters found for this roc | |
56c85970 | 2261 | k0 = k1 = k2 = l0 = l1 = l2 = -1; |
df4cfd77 | 2262 | } |
2263 | ||
2264 | fileOut << entryNum++ << "\t" << side << "\t" << sector << "\t" << rcu << "\t" << -1 << "\t"; | |
56c85970 | 2265 | fileOut << k0 << "\t" << k1 << "\t" << k2 << "\t" << l0 << "\t" << l1 << "\t" << l2 << "\t" << validFlag << std::endl; |
df4cfd77 | 2266 | } |
2267 | } | |
2268 | } | |
2269 | ||
2270 | printf("done assigning\n"); | |
2271 | ||
2272 | // closing files | |
2273 | fileOut.close(); | |
2274 | fileTCFparam.Close(); | |
2275 | printf("output written to file: %s\n",nameFileOut); | |
2276 | return 0; | |
2277 | ||
2278 | } |