]>
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 | // // | |
22 | // Author: Stefan Rossegger // | |
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> | |
35 | ||
36 | #include <TMath.h> | |
37 | #include <TNtuple.h> | |
38 | #include <TEntryList.h> | |
39 | ||
40 | #include "AliRawReaderRoot.h" | |
41 | #include "AliTPCRawStream.h" | |
42 | #include "AliTPCROC.h" | |
43 | ||
44 | #include "AliTPCAltroEmulator.h" | |
45 | ||
46 | ClassImp(AliTPCCalibTCF) | |
47 | ||
48 | AliTPCCalibTCF::AliTPCCalibTCF() : | |
49 | TNamed(), | |
50 | fGateWidth(100), | |
51 | fSample(900), | |
52 | fPulseLength(500), | |
53 | fLowPulseLim(30), | |
54 | fUpPulseLim(1000), | |
55 | fRMSLim(4.) | |
56 | { | |
57 | // | |
58 | // AliTPCCalibTCF standard constructor | |
59 | // | |
60 | } | |
61 | ||
62 | //_____________________________________________________________________________ | |
63 | AliTPCCalibTCF::AliTPCCalibTCF(Int_t gateWidth, Int_t sample, Int_t pulseLength, Int_t lowPulseLim, Int_t upPulseLim, Double_t rmsLim) : | |
64 | TNamed(), | |
65 | fGateWidth(gateWidth), | |
66 | fSample(sample), | |
67 | fPulseLength(pulseLength), | |
68 | fLowPulseLim(lowPulseLim), | |
69 | fUpPulseLim(upPulseLim), | |
70 | fRMSLim(rmsLim) | |
71 | { | |
72 | // | |
73 | // AliTPCCalibTCF constructor with specific (non-standard) thresholds | |
74 | // | |
75 | } | |
76 | ||
77 | //_____________________________________________________________________________ | |
78 | AliTPCCalibTCF::AliTPCCalibTCF(const AliTPCCalibTCF &tcf) : | |
79 | TNamed(tcf), | |
80 | fGateWidth(tcf.fGateWidth), | |
81 | fSample(tcf.fSample), | |
82 | fPulseLength(tcf.fPulseLength), | |
83 | fLowPulseLim(tcf.fLowPulseLim), | |
84 | fUpPulseLim(tcf.fUpPulseLim), | |
85 | fRMSLim(tcf.fRMSLim) | |
86 | { | |
87 | // | |
88 | // AliTPCCalibTCF copy constructor | |
89 | // | |
90 | } | |
91 | ||
92 | ||
93 | //_____________________________________________________________________________ | |
94 | AliTPCCalibTCF& AliTPCCalibTCF::operator = (const AliTPCCalibTCF &source) | |
95 | { | |
96 | // | |
97 | // AliTPCCalibTCF assignment operator | |
98 | // | |
99 | ||
100 | if (&source == this) return *this; | |
101 | new (this) AliTPCCalibTCF(source); | |
102 | ||
103 | return *this; | |
104 | ||
105 | } | |
106 | ||
107 | //_____________________________________________________________________________ | |
108 | AliTPCCalibTCF::~AliTPCCalibTCF() | |
109 | { | |
110 | // | |
111 | // AliTPCCalibTCF destructor | |
112 | // | |
113 | } | |
114 | ||
115 | //_____________________________________________________________________________ | |
116 | void AliTPCCalibTCF::ProcessRawFile(const char *nameRawFile, const char *nameFileOut) { | |
117 | // | |
118 | // Loops over all events within one RawData file and collects proper pulses | |
119 | // (according to given tresholds) per pad | |
120 | // Histograms per pad are stored in 'nameFileOut' | |
121 | // | |
122 | ||
123 | AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile); | |
124 | rawReader->Reset(); | |
125 | ||
126 | while ( rawReader->NextEvent() ){ // loop | |
127 | printf("Reading next event ..."); | |
128 | AliTPCRawStream rawStream(rawReader); | |
129 | rawReader->Select("TPC"); | |
130 | ProcessRawEvent(&rawStream, nameFileOut); | |
131 | } | |
132 | ||
133 | rawReader->~AliRawReader(); | |
134 | ||
135 | } | |
136 | ||
137 | ||
138 | //_____________________________________________________________________________ | |
139 | void AliTPCCalibTCF::ProcessRawEvent(AliTPCRawStream *rawStream, const char *nameFileOut) { | |
140 | // | |
141 | // Extracts proper pulses (according the given tresholds) within one event | |
142 | // and accumulates them into one histogram per pad. All histograms are | |
143 | // saved in the file 'nameFileOut'. | |
144 | // The first bins of the histograms contain the following information: | |
145 | // bin 1: Number of accumulated pulses | |
146 | // bin 2;3;4: Sector; Row; Pad; | |
147 | // | |
148 | ||
149 | Int_t sector = rawStream->GetSector(); | |
150 | Int_t row = rawStream->GetRow(); | |
151 | Int_t prevTime = 999999; | |
152 | Int_t prevPad = 999999; | |
153 | ||
154 | TFile fileOut(nameFileOut,"UPDATE"); | |
155 | fileOut.cd(); | |
156 | ||
157 | TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth); | |
158 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000); | |
159 | ||
160 | rawStream->SetOldRCUFormat(1); | |
161 | ||
162 | while (rawStream->Next()) { | |
163 | ||
164 | // in case of a new row, get sector and row number | |
165 | if (rawStream->IsNewRow()){ | |
166 | sector = rawStream->GetSector(); | |
167 | row = rawStream->GetRow(); | |
168 | } | |
169 | ||
170 | Int_t pad = rawStream->GetPad(); | |
171 | Int_t time = rawStream->GetTime(); | |
172 | Int_t signal = rawStream->GetSignal(); | |
173 | ||
174 | if (!rawStream->IsNewPad()) { // Reading signal from one Pad | |
175 | if (time>prevTime) { | |
176 | printf("Wrong time: %d %d\n",rawStream->GetTime(),prevTime); | |
177 | rawStream->Dump(); | |
178 | } else { | |
179 | // still the same pad, save signal to temporary histogram | |
180 | if (time<=fSample+fGateWidth && time>fGateWidth) { | |
181 | tempHis->SetBinContent(time,signal); | |
182 | } | |
183 | } | |
184 | } else { | |
185 | // complete pulse found and stored into tempHis, now calculation | |
186 | // of it's properties and comparison to given thresholds | |
187 | ||
188 | Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX); | |
189 | Int_t maxpos = tempHis->GetMaximumBin(); | |
190 | ||
191 | Int_t first = (Int_t)TMath::Max(maxpos-10, 0); | |
192 | Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample); | |
193 | ||
194 | // simple baseline substraction ? better one needed ? (pedestalsubstr.?) | |
195 | // and RMS calculation with timebins before the pulse and at the end of | |
196 | // the signal | |
197 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
198 | // before the pulse | |
199 | tempRMSHis->Fill(tempHis->GetBinContent(first+ipos)); | |
200 | // at the end to get rid of pulses with serious baseline fluctuations | |
201 | tempRMSHis->Fill(tempHis->GetBinContent(last-ipos)); | |
202 | } | |
203 | Double_t baseline = tempRMSHis->GetMean(); | |
204 | Double_t rms = tempRMSHis->GetRMS(); | |
205 | tempRMSHis->Reset(); | |
206 | ||
207 | Double_t lowLim = fLowPulseLim+baseline; | |
208 | Double_t upLim = fUpPulseLim+baseline; | |
209 | ||
210 | // Decision if found pulse is a proper one according to given tresholds | |
211 | if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim){ | |
212 | char hname[100]; | |
213 | sprintf(hname,"sec%drow%dpad%d",sector,row,prevPad); | |
214 | ||
215 | TH1F *his = (TH1F*)fileOut.Get(hname); | |
216 | ||
217 | if (!his ) { // new entry (pulse in new pad found) | |
218 | ||
219 | his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4); | |
220 | his->SetBinContent(1,1); // pulse counter (1st pulse) | |
221 | his->SetBinContent(2,sector); // sector | |
222 | his->SetBinContent(3,row); // row | |
223 | his->SetBinContent(4,prevPad); // pad | |
224 | ||
225 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
226 | Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); | |
227 | his->SetBinContent(ipos+5,signal); | |
228 | } | |
229 | his->Write(hname); | |
230 | printf("new %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
231 | ||
232 | } else { // adding pulse to existing histogram (pad already found) | |
233 | ||
234 | his->AddBinContent(1,1); // pulse counter for each pad | |
235 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
236 | Int_t signal= (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); | |
237 | his->AddBinContent(ipos+5,signal); | |
238 | } | |
239 | printf("adding ... %s: Signal %d at bin %d \n", hname, max-(Int_t)baseline, maxpos+fGateWidth); | |
240 | his->Write(hname,kOverwrite); | |
241 | } | |
242 | } | |
243 | tempHis->Reset(); | |
244 | } | |
245 | prevTime = time; | |
246 | prevPad = pad; | |
247 | } | |
248 | ||
249 | tempHis->~TH1I(); | |
250 | tempRMSHis->~TH1I(); | |
251 | printf("Finished to read event ... \n"); | |
252 | fileOut.Close(); | |
253 | } | |
254 | ||
255 | //____________________________________________________________________________ | |
256 | void AliTPCCalibTCF::MergeHistoPerSector(const char *nameFileIn) { | |
257 | // | |
258 | // Merges all histograms within one sector, calculates the TCF parameters | |
259 | // of the 'histogram-per-sector' and stores (histo and parameters) into | |
260 | // seperated files ... | |
261 | // | |
262 | // note: first 4 timebins of a histogram hold specific informations | |
263 | // about number of collected pulses, sector, row and pad | |
264 | // | |
265 | // 'nameFileIn': root file produced with Process function which holds | |
266 | // one histogram per pad (sum of signals of proper pulses) | |
267 | // 'Sec+nameFileIn': root file with one histogram per sector | |
268 | // (information of row and pad are set to -1) | |
269 | // | |
270 | ||
271 | TFile fileIn(nameFileIn,"READ"); | |
272 | TH1F *hisPad = 0; | |
273 | TKey *key = 0; | |
274 | TIter next( fileIn.GetListOfKeys() ); | |
275 | ||
276 | char nameFileOut[100]; | |
277 | sprintf(nameFileOut,"Sec-%s",nameFileIn); | |
278 | ||
279 | TFile fileOut(nameFileOut,"RECREATE"); | |
280 | fileOut.cd(); | |
281 | ||
282 | Int_t nHist = fileIn.GetNkeys(); | |
283 | Int_t iHist = 0; // histogram counter for merge-status print | |
284 | ||
285 | while ( (key=(TKey*)next()) ) { | |
286 | ||
287 | iHist++; | |
288 | ||
289 | hisPad = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory | |
290 | Int_t pulseLength = hisPad->GetNbinsX() -4; | |
291 | // -4 because first four timebins contain pad specific informations | |
292 | Int_t npulse = (Int_t)hisPad->GetBinContent(1); | |
293 | Int_t sector = (Int_t)hisPad->GetBinContent(2); | |
294 | ||
295 | char hname[100]; | |
296 | sprintf(hname,"sector%d",sector); | |
297 | TH1F *his = (TH1F*)fileOut.Get(hname); | |
298 | ||
299 | if (!his ) { // new histogram (new sector) | |
300 | his = new TH1F(hname,hname, pulseLength+4, 0, pulseLength+4); | |
301 | his->SetBinContent(1,npulse); // pulse counter | |
302 | his->SetBinContent(2,sector); // set sector info | |
303 | his->SetBinContent(3,-1); // set to dummy value | |
304 | his->SetBinContent(4,-1); // set to dummy value | |
305 | for (Int_t ipos=0; ipos<pulseLength; ipos++){ | |
306 | his->SetBinContent(ipos+5,hisPad->GetBinContent(ipos+5)); | |
307 | } | |
308 | his->Write(hname); | |
309 | printf("found %s ...\n", hname); | |
310 | } else { // add to existing histogram for sector | |
311 | his->AddBinContent(1,npulse); // pulse counter | |
312 | for (Int_t ipos=0; ipos<pulseLength; ipos++){ | |
313 | his->AddBinContent(ipos+5,hisPad->GetBinContent(ipos+5)); | |
314 | } | |
315 | his->Write(hname,kOverwrite); | |
316 | } | |
317 | ||
318 | if (iHist%500==0) { | |
319 | printf("merging status: \t %d pads out of %d \n",iHist, nHist); | |
320 | } | |
321 | } | |
322 | printf("merging done ...\n"); | |
323 | fileIn.Close(); | |
324 | fileOut.Close(); | |
325 | ||
326 | // calculate TCF parameters on averaged pulse per Sector | |
327 | AnalyzeRootFile(nameFileOut); | |
328 | ||
329 | ||
330 | } | |
331 | ||
332 | ||
333 | //____________________________________________________________________________ | |
334 | void AliTPCCalibTCF::AnalyzeRootFile(const char *nameFileIn, Int_t minNumPulse) { | |
335 | // | |
336 | // This function takes a prepeared root file (accumulated histograms: output | |
337 | // of process function) and performs an analysis (fit and equalization) in | |
338 | // order to get the TCF parameters. These are stored in an TNtuple along with | |
339 | // the pad and creation infos. The tuple is written to the output file | |
340 | // "TCFparam+nameFileIn" | |
341 | // To reduce the analysis time, the minimum number of accumulated pulses within | |
342 | // one histogram 'minNumPulse' (to perform the analysis on) can be set | |
343 | // | |
344 | ||
345 | TFile fileIn(nameFileIn,"READ"); | |
346 | TH1F *hisIn; | |
347 | TKey *key; | |
348 | TIter next( fileIn.GetListOfKeys() ); | |
349 | ||
350 | char nameFileOut[100]; | |
351 | sprintf(nameFileOut,"TCFparam-%s",nameFileIn); | |
352 | ||
353 | TFile fileOut(nameFileOut,"RECREATE"); | |
354 | fileOut.cd(); | |
355 | ||
356 | TNtuple *paramTuple = new TNtuple("TCFparam","TCFparameter","sec:row:pad:npulse:Z0:Z1:Z2:P0:P1:P2"); | |
357 | ||
358 | Int_t nHist = fileIn.GetNkeys(); | |
359 | Int_t iHist = 0; // counter for print of analysis-status | |
360 | ||
9389f9a4 | 361 | while ((key = (TKey *) next())) { // loop over histograms |
eb7e0771 | 362 | |
363 | printf("Analyze histogramm %d out of %d\n",++iHist,nHist); | |
364 | hisIn = (TH1F*)fileIn.Get(key->GetName()); // copy object to memory | |
365 | ||
366 | Int_t numPulse = (Int_t)hisIn->GetBinContent(1); | |
367 | if ( numPulse >= minNumPulse ) { | |
368 | ||
369 | Double_t* coefP = new Double_t[3]; | |
370 | Double_t* coefZ = new Double_t[3]; | |
371 | for(Int_t i = 0; i < 3; i++){ | |
372 | coefP[i] = 0; | |
373 | coefZ[i] = 0; | |
374 | } | |
375 | // perform the analysis on the given histogram | |
376 | Int_t fitOk = AnalyzePulse(hisIn, coefZ, coefP); | |
377 | if (fitOk) { // Add found parameters to file | |
378 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
379 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
380 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
381 | paramTuple->Fill(sector,row,pad,numPulse,coefZ[0],coefZ[1],coefZ[2],coefP[0],coefP[1],coefP[2]); | |
382 | } | |
383 | coefP->~Double_t(); | |
384 | coefZ->~Double_t(); | |
385 | } | |
386 | ||
387 | } | |
388 | ||
389 | fileIn.Close(); | |
390 | paramTuple->Write(); | |
391 | fileOut.Close(); | |
392 | ||
393 | } | |
394 | ||
395 | ||
396 | //____________________________________________________________________________ | |
397 | Int_t AliTPCCalibTCF::AnalyzePulse(TH1F *hisIn, Double_t *coefZ, Double_t *coefP) { | |
398 | // | |
399 | // Performs the analysis on one specific pulse (histogram) by means of fitting | |
400 | // the pulse and equalization of the pulseheight. The found TCF parameters | |
401 | // are stored in the arrays coefZ and coefP | |
402 | // | |
403 | ||
404 | Int_t pulseLength = hisIn->GetNbinsX() -4; | |
405 | // -1 because the first four timebins usually contain pad specific informations | |
406 | Int_t npulse = (Int_t)hisIn->GetBinContent(1); | |
407 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
408 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
409 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
410 | ||
411 | // write pulseinformation to TNtuple and normalize to 100 ADC (because of | |
412 | // given upper and lower fit parameter limits) in order to pass the pulse | |
413 | // to TMinuit | |
414 | ||
415 | TNtuple *dataTuple = new TNtuple("ntupleFit","Pulse","timebin:sigNorm:error"); | |
416 | Double_t error = 0.05; | |
417 | Double_t max = hisIn->GetMaximum(FLT_MAX); | |
418 | for (Int_t ipos=0; ipos<pulseLength; ipos++) { | |
419 | Double_t errorz=error; | |
420 | if (ipos>100) { errorz = error*100; } // very simple weight: FIXME in case | |
421 | Double_t signal = hisIn->GetBinContent(ipos+5); | |
422 | Double_t signalNorm = signal/max*100; //pulseheight normaliz. to 100ADC | |
423 | dataTuple->Fill(ipos, signalNorm, errorz); | |
424 | } | |
425 | ||
426 | // Call fit function (TMinuit) to get the first 2 PZ Values for the | |
427 | // Tail Cancelation Filter | |
428 | Int_t fitOk = FitPulse(dataTuple, coefZ, coefP); | |
429 | ||
430 | if (fitOk) { | |
431 | // calculates the 3rd set (remaining 2 PZ values) in order to restore the | |
432 | // original height of the pulse | |
433 | Equalization(dataTuple, coefZ, coefP); | |
434 | ||
435 | printf("Calculated TCF parameters for: \n"); | |
436 | printf("Sector %d | Row %d | Pad %d |", sector, row, pad); | |
437 | printf(" Npulses: %d \n", npulse); | |
438 | for(Int_t i = 0; i < 3; i++){ | |
439 | printf("P[%d] = %f Z[%d] = %f \n",i,coefP[i],i,coefZ[i]); | |
440 | if (i==2) { printf("\n"); } | |
441 | } | |
442 | dataTuple->~TNtuple(); | |
443 | return 1; | |
444 | } else { // fit did not converge | |
445 | Error("FindFit", "TCF fit not converged - pulse abandoned "); | |
446 | printf("in Sector %d | Row %d | Pad %d |", sector, row, pad); | |
447 | printf(" Npulses: %d \n\n", npulse); | |
448 | coefP[2] = 0; coefZ[2] = 0; | |
449 | dataTuple->~TNtuple(); | |
450 | return 0; | |
451 | } | |
452 | ||
453 | } | |
454 | ||
455 | ||
456 | ||
457 | //____________________________________________________________________________ | |
458 | void AliTPCCalibTCF::TestTCFonRootFile(const char *nameFileIn, const char *nameFileTCF, Int_t plotFlag, Int_t lowKey, Int_t upKey) | |
459 | { | |
460 | // | |
461 | // Performs quality parameters evaluation of the calculated TCF parameters in | |
462 | // the file 'nameFileTCF' for every (accumulated) histogram within the | |
463 | // prepeared root file 'nameFileIn'. | |
464 | // The found quality parameters are stored in an TNtuple which will be saved | |
465 | // in a Root file 'Quality-*'. | |
466 | // If the parameter for the given pulse (given pad) was not found, the pulse | |
467 | // is rejected. | |
468 | // | |
469 | ||
470 | TFile fileIn(nameFileIn,"READ"); | |
471 | ||
472 | Double_t* coefP = new Double_t[3]; | |
473 | Double_t* coefZ = new Double_t[3]; | |
474 | for(Int_t i = 0; i < 3; i++){ | |
475 | coefP[i] = 0; | |
476 | coefZ[i] = 0; | |
477 | } | |
478 | ||
479 | char nameFileOut[100]; | |
480 | sprintf(nameFileOut,"Quality_%s_AT_%s",nameFileTCF, nameFileIn); | |
481 | TFile fileOut(nameFileOut,"RECREATE"); | |
482 | ||
483 | TNtuple *qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot"); | |
484 | ||
485 | TH1F *hisIn; | |
486 | TKey *key; | |
487 | TIter next( fileIn.GetListOfKeys() ); | |
488 | ||
489 | Int_t nHist = fileIn.GetNkeys(); | |
490 | Int_t iHist = 0; | |
491 | ||
492 | for(Int_t i=0;i<lowKey-1;i++){++iHist; key = (TKey *) next();} | |
2c632057 | 493 | while ((key = (TKey *) next())) { // loop over saved histograms |
eb7e0771 | 494 | |
495 | // loading pulse to memory; | |
496 | printf("validating pulse %d out of %d\n",++iHist,nHist); | |
497 | hisIn = (TH1F*)fileIn.Get(key->GetName()); | |
498 | ||
499 | // find the correct TCF parameter according to the his infos (first 4 bins) | |
500 | Int_t nPulse = FindCorTCFparam(hisIn, nameFileTCF, coefZ, coefP); | |
501 | if (nPulse) { // doing the TCF quality analysis | |
502 | Double_t *quVal = GetQualityOfTCF(hisIn,coefZ,coefP, plotFlag); | |
503 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
504 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
505 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
506 | qualityTuple->Fill(sector,row,pad,nPulse,quVal[0],quVal[1],quVal[2],quVal[3],quVal[4],quVal[5]); | |
507 | quVal->~Double_t(); | |
508 | } | |
509 | ||
510 | if (iHist>=upKey) {break;} | |
511 | ||
512 | } | |
513 | ||
514 | fileOut.cd(); | |
515 | qualityTuple->Write(); | |
516 | ||
517 | coefP->~Double_t(); | |
518 | coefZ->~Double_t(); | |
519 | ||
520 | fileOut.Close(); | |
521 | fileIn.Close(); | |
522 | ||
523 | } | |
524 | ||
525 | ||
526 | ||
527 | //_____________________________________________________________________________ | |
528 | void AliTPCCalibTCF::TestTCFonRawFile(const char *nameRawFile, const char *nameFileOut, const char *nameFileTCF, Int_t plotFlag) { | |
529 | // | |
530 | // Performs quality parameters evaluation of the calculated TCF parameters in | |
531 | // the file 'nameFileTCF' for every proper pulse (according to given thresholds) | |
532 | // within the RAW file 'nameRawFile'. | |
533 | // The found quality parameters are stored in a TNtuple which will be saved | |
534 | // in the Root file 'nameFileOut'. If the parameter for the given pulse | |
535 | // (given pad) was not found, the pulse is rejected. | |
536 | // | |
537 | ||
538 | // | |
539 | // Reads a RAW data file, extracts Pulses (according the given tresholds) | |
540 | // and test the found TCF parameters on them ... | |
541 | // | |
542 | ||
543 | AliRawReader *rawReader = new AliRawReaderRoot(nameRawFile); | |
544 | rawReader->Reset(); | |
545 | ||
546 | Double_t* coefP = new Double_t[3]; | |
547 | Double_t* coefZ = new Double_t[3]; | |
548 | for(Int_t i = 0; i < 3; i++){ | |
549 | coefP[i] = 0; | |
550 | coefZ[i] = 0; | |
551 | } | |
552 | ||
553 | while ( rawReader->NextEvent() ){ | |
554 | ||
555 | printf("Reading next event..."); | |
556 | AliTPCRawStream rawStream(rawReader); | |
557 | rawReader->Select("TPC"); | |
558 | ||
559 | Int_t sector = rawStream.GetSector(); | |
560 | Int_t row = rawStream.GetRow(); | |
561 | Int_t prevTime = 999999; | |
562 | Int_t prevPad = 999999; | |
563 | ||
564 | TH1I *tempHis = new TH1I("tempHis","tempHis",fSample+fGateWidth,fGateWidth,fSample+fGateWidth); | |
565 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",2000,0,2000); | |
566 | ||
567 | rawStream.SetOldRCUFormat(1); | |
568 | ||
569 | TFile fileOut(nameFileOut,"UPDATE"); // Quality Parameters storage | |
570 | TNtuple *qualityTuple = (TNtuple*)fileOut.Get("TCFquality"); | |
571 | if (!qualityTuple) { // no entry in file | |
572 | qualityTuple = new TNtuple("TCFquality","TCF quality Values","sec:row:pad:npulse:heightDev:areaRed:widthRed:undershot:maxUndershot:pulseRMS"); | |
573 | } | |
574 | ||
575 | while (rawStream.Next()) { | |
576 | ||
577 | if (rawStream.IsNewRow()){ | |
578 | sector = rawStream.GetSector(); | |
579 | row = rawStream.GetRow(); | |
580 | } | |
581 | ||
582 | Int_t pad = rawStream.GetPad(); | |
583 | Int_t time = rawStream.GetTime(); | |
584 | Int_t signal = rawStream.GetSignal(); | |
585 | ||
586 | if (!rawStream.IsNewPad()) { // Reading signal from one Pad | |
587 | if (time>prevTime) { | |
588 | printf("Wrong time: %d %d\n",rawStream.GetTime(),prevTime); | |
589 | rawStream.Dump(); | |
590 | } else { | |
591 | if (time<=fSample+fGateWidth && time>fGateWidth) { | |
592 | tempHis->SetBinContent(time,signal); | |
593 | } | |
594 | } | |
595 | } else { // Decision for saving pulse according to treshold settings | |
596 | ||
597 | Int_t max = (Int_t)tempHis->GetMaximum(FLT_MAX); | |
598 | Int_t maxpos = tempHis->GetMaximumBin(); | |
599 | ||
600 | Int_t first = (Int_t)TMath::Max(maxpos-10, 0); | |
601 | Int_t last = TMath::Min((Int_t)maxpos+fPulseLength-10, fSample); | |
602 | ||
603 | ||
604 | // simple baseline substraction ? better one needed ? (pedestalsubstr.?) | |
605 | // and RMS calculation with timebins before the pulse and at the end of | |
606 | // the signal | |
607 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
608 | // before the pulse | |
609 | tempRMSHis->Fill(tempHis->GetBinContent(first+ipos)); | |
610 | // at the end to get rid of pulses with serious baseline fluctuations | |
611 | tempRMSHis->Fill(tempHis->GetBinContent(last-ipos)); | |
612 | } | |
613 | Double_t baseline = tempRMSHis->GetMean(); | |
614 | Double_t rms = tempRMSHis->GetRMS(); | |
615 | tempRMSHis->Reset(); | |
616 | ||
617 | Double_t lowLim = fLowPulseLim+baseline; | |
618 | Double_t upLim = fUpPulseLim+baseline; | |
619 | ||
620 | // Decision if found pulse is a proper one according to given tresholds | |
621 | if (max>lowLim && max<upLim && !((last-first)<fPulseLength) && rms<fRMSLim){ | |
622 | // note: | |
623 | // assuming that lowLim is higher than the pedestal value! | |
624 | char hname[100]; | |
625 | sprintf(hname,"sec%drow%dpad%d",sector,row,prevPad); | |
626 | TH1F *his = new TH1F(hname,hname, fPulseLength+4, 0, fPulseLength+4); | |
627 | his->SetBinContent(1,1); // pulse counter (1st pulse) | |
628 | his->SetBinContent(2,sector); // sector | |
629 | his->SetBinContent(3,row); // row | |
630 | his->SetBinContent(4,prevPad); // pad | |
631 | for (Int_t ipos=0; ipos<last-first; ipos++){ | |
632 | Int_t signal = (Int_t)(tempHis->GetBinContent(ipos+first)-baseline); | |
633 | his->SetBinContent(ipos+5,signal); | |
634 | } | |
635 | ||
636 | printf("Pulse found in %s: ADC %d at bin %d \n", hname, max, maxpos+fGateWidth); | |
637 | ||
638 | // find the correct TCF parameter according to the his infos | |
639 | // (first 4 bins) | |
640 | Int_t nPulse = FindCorTCFparam(his, nameFileTCF, coefZ, coefP); | |
641 | ||
642 | if (nPulse) { // Parameters found - doing the TCF quality analysis | |
643 | Double_t *quVal = GetQualityOfTCF(his,coefZ,coefP, plotFlag); | |
644 | qualityTuple->Fill(sector,row,pad,nPulse,quVal[0],quVal[1],quVal[2],quVal[3],quVal[4],quVal[5]); | |
645 | quVal->~Double_t(); | |
646 | } | |
647 | his->~TH1F(); | |
648 | } | |
649 | tempHis->Reset(); | |
650 | } | |
651 | prevTime = time; | |
652 | prevPad = pad; | |
653 | } | |
654 | ||
655 | tempHis->~TH1I(); | |
656 | tempRMSHis->~TH1I(); | |
657 | ||
658 | printf("Finished to read event - close output file ... \n"); | |
659 | ||
660 | fileOut.cd(); | |
661 | qualityTuple->Write("TCFquality",kOverwrite); | |
662 | fileOut.Close(); | |
663 | ||
664 | ||
665 | ||
666 | } // event loop | |
667 | ||
668 | ||
669 | coefP->~Double_t(); | |
670 | coefZ->~Double_t(); | |
671 | ||
672 | rawReader->~AliRawReader(); | |
673 | ||
674 | } | |
675 | ||
676 | ||
677 | //____________________________________________________________________________ | |
678 | TNtuple *AliTPCCalibTCF::PlotOccupSummary(const char *nameFile, Int_t nPulseMin) { | |
679 | // | |
680 | // Plots the number of summed pulses per pad above a given minimum at the | |
681 | // pad position | |
682 | // 'nameFile': root-file created with the Process function | |
683 | // | |
684 | ||
685 | TFile *file = new TFile(nameFile,"READ"); | |
686 | ||
687 | TH1F *his; | |
688 | TKey *key; | |
689 | TIter next( file->GetListOfKeys() ); | |
690 | ||
691 | TNtuple *ntuple = new TNtuple("ntuple","ntuple","x:y:z:npulse"); | |
692 | ||
693 | Int_t nPads = 0; | |
694 | while ((key = (TKey *) next())) { // loop over histograms within the file | |
695 | ||
696 | his = (TH1F*)file->Get(key->GetName()); // copy object to memory | |
697 | ||
698 | Int_t npulse = (Int_t)his->GetBinContent(1); | |
699 | Int_t sec = (Int_t)his->GetBinContent(2); | |
700 | Int_t row = (Int_t)his->GetBinContent(3); | |
701 | Int_t pad = (Int_t)his->GetBinContent(4); | |
702 | ||
703 | if (row==-1 & pad==-1) { // summed pulses per sector | |
704 | row = 40; pad = 40; // set to approx middle row for better plot | |
705 | } | |
706 | ||
707 | Float_t *pos = new Float_t[3]; | |
708 | // find x,y,z position of the pad | |
709 | AliTPCROC::Instance()->GetPositionGlobal(sec,row,pad,pos); | |
710 | if (npulse>=nPulseMin) { | |
711 | ntuple->Fill(pos[0],pos[1],pos[2],npulse); | |
712 | printf("%d collected pulses in sector %d row %d pad %d\n",npulse,sec,row,pad); | |
713 | } | |
714 | pos->~Float_t(); | |
715 | nPads++; | |
716 | } | |
717 | ||
718 | TCanvas *c1 = new TCanvas("TCanvas","Number of pulses found",1000,500); | |
719 | c1->Divide(2,1); | |
720 | char cSel[100]; | |
721 | gStyle->SetPalette(1); | |
722 | gStyle->SetLabelOffset(-0.03,"Z"); | |
723 | ||
724 | if (nPads<72) { // pulse per pad | |
725 | ntuple->SetMarkerStyle(8); | |
726 | ntuple->SetMarkerSize(4); | |
727 | } else { // pulse per sector | |
728 | ntuple->SetMarkerStyle(7); | |
729 | } | |
730 | ||
731 | c1->cd(1); | |
732 | sprintf(cSel,"z>0&&npulse>=%d",nPulseMin); | |
733 | ntuple->Draw("y:x:npulse",cSel,"colz"); | |
734 | gPad->SetTitle("A side"); | |
735 | ||
736 | c1->cd(2); | |
737 | sprintf(cSel,"z<0&&npulse>%d",nPulseMin); | |
738 | ntuple->Draw("y:x:npulse",cSel,"colz"); | |
739 | gPad->SetTitle("C side"); | |
740 | ||
741 | file->Close(); | |
742 | return ntuple; | |
743 | ||
744 | } | |
745 | ||
746 | //____________________________________________________________________________ | |
747 | void AliTPCCalibTCF::PlotQualitySummary(const char *nameFileQuality, const char *plotSpec, const char *cut, const char *pOpt) | |
748 | { | |
749 | // | |
750 | // This function is an easy interface to load the QualityTuple (produced with | |
751 | // the function 'TestOn%File' and plots them according to the plot specifications | |
752 | // 'plotSpec' e.g. "widthRed:maxUndershot" | |
753 | // One may also set cut and plot options ("cut","pOpt") | |
754 | // | |
755 | // The stored quality parameters are ... | |
756 | // sec:row:pad:npulse: ... usual pad info | |
757 | // heightDev ... height deviation in percent | |
758 | // areaRed ... area reduction in percent | |
759 | // widthRed ... width reduction in percent | |
760 | // undershot ... mean undershot after the pulse in ADC | |
761 | // maxUndershot ... maximum of the undershot after the pulse in ADC | |
762 | // pulseRMS ... RMS of the pulse used to calculate the Quality parameters in ADC | |
763 | // | |
764 | ||
765 | TFile file(nameFileQuality,"READ"); | |
766 | TNtuple *qualityTuple = (TNtuple*)file.Get("TCFquality"); | |
767 | gStyle->SetPalette(1); | |
768 | qualityTuple->Draw(plotSpec,cut,pOpt); | |
769 | ||
770 | } | |
771 | ||
772 | //____________________________________________________________________________ | |
773 | void AliTPCCalibTCF::DumpTCFparamToFile(const char *nameFileTCF,const char *nameFileOut) | |
774 | { | |
775 | // | |
776 | // Writes the TCF parameters from file 'nameFileTCF' to a output file | |
777 | // | |
778 | ||
779 | // Note: currently just TCF parameters per Sector or TCF parameters for pad | |
780 | // which were analyzed. There is no method included so far to export | |
781 | // parameters for not analyzed pad, which means there are eventually | |
782 | // missing TCF parameters | |
783 | // TODO: carefull! Fill up missing pads with averaged (sector) values? | |
784 | ||
785 | ||
786 | // open file with TCF parameters | |
787 | TFile fileTCF(nameFileTCF,"READ"); | |
788 | TNtuple *paramTuple = (TNtuple*)fileTCF.Get("TCFparam"); | |
789 | ||
790 | // open output txt file ... | |
791 | FILE *output; | |
792 | output=fopen(nameFileOut,"w"); // open outfile. | |
793 | ||
794 | // Header line | |
795 | Int_t sectorWise = paramTuple->GetEntries("row==-1&&pad==-1"); | |
796 | if (sectorWise) { | |
797 | fprintf(output,"sector \t Z0 \t\t Z1 \t\t Z2 \t\t P0 \t\t P1 \t\t P2\n"); | |
798 | } else { | |
799 | fprintf(output,"sector \t row \t pad \t Z0 \t\t Z1 \t\t Z2 \t\t P0 \t\t P1 \t\t P2\n"); | |
800 | } | |
801 | ||
802 | for (Int_t i=0; i<paramTuple->GetEntries(); i++) { | |
803 | paramTuple->GetEntry(i); | |
804 | Float_t *p = paramTuple->GetArgs(); | |
805 | ||
806 | // _______________________________________________________________ | |
807 | // to Tuple to txt file - unsorted printout | |
808 | ||
809 | for (Int_t i=0; i<10; i++){ | |
810 | if (sectorWise) { | |
811 | if (i<1) fprintf(output,"%3.0f \t ",p[i]); // sector info | |
812 | if (i>3) fprintf(output,"%1.4f \t ",p[i]); // TCF param | |
813 | } else { | |
814 | if (i<3) fprintf(output,"%3.0f \t ",p[i]); // pad info | |
815 | if (i>3) fprintf(output,"%1.4f \t ",p[i]); // TCF param | |
816 | } | |
817 | } | |
818 | fprintf(output,"\n"); | |
819 | } | |
820 | ||
821 | // close output txt file | |
822 | fprintf(output,"\n"); | |
823 | fclose(output); | |
824 | ||
825 | fileTCF.Close(); | |
826 | ||
827 | ||
828 | } | |
829 | ||
830 | ||
831 | ||
832 | //_____________________________________________________________________________ | |
833 | Int_t AliTPCCalibTCF::FitPulse(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) { | |
834 | // | |
835 | // function to fit one pulse and to calculate the according pole-zero parameters | |
836 | // | |
837 | ||
838 | // initialize TMinuit with a maximum of 8 params | |
839 | TMinuit *gMinuit = new TMinuit(8); | |
840 | gMinuit->mncler(); // Reset Minuit's list of paramters | |
841 | gMinuit->SetPrintLevel(-1); // No Printout | |
842 | gMinuit->SetFCN(AliTPCCalibTCF::FitFcn); // To set the address of the | |
843 | // minimization function | |
844 | gMinuit->SetObjectFit(dataTuple); | |
845 | ||
846 | Double_t arglist[10]; | |
847 | Int_t ierflg = 0; | |
848 | ||
849 | arglist[0] = 1; | |
850 | gMinuit->mnexcm("SET ERR", arglist ,1,ierflg); | |
851 | ||
852 | // Set standard starting values and step sizes for each parameter | |
853 | // upper and lower limit (in a reasonable range) are set to improve | |
854 | // the stability of TMinuit | |
855 | static Double_t vstart[8] = {125, 4.0, 0.3, 0.5, 5.5, 100, 1, 2.24}; | |
856 | static Double_t step[8] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1}; | |
857 | static Double_t min[8] = {100, 3., 0.1, 0.2, 3., 60., 0., 2.0}; | |
858 | static Double_t max[8] = {200, 20., 5., 3., 30., 300., 20., 2.5}; | |
859 | ||
860 | gMinuit->mnparm(0, "A1", vstart[0], step[0], min[0], max[0], ierflg); | |
861 | gMinuit->mnparm(1, "A2", vstart[1], step[1], min[1], max[1], ierflg); | |
862 | gMinuit->mnparm(2, "A3", vstart[2], step[2], min[2], max[2], ierflg); | |
863 | gMinuit->mnparm(3, "T1", vstart[3], step[3], min[3], max[3], ierflg); | |
864 | gMinuit->mnparm(4, "T2", vstart[4], step[4], min[4], max[4], ierflg); | |
865 | gMinuit->mnparm(5, "T3", vstart[5], step[5], min[5], max[5], ierflg); | |
866 | gMinuit->mnparm(6, "T0", vstart[6], step[6], min[6], max[6], ierflg); | |
867 | gMinuit->mnparm(7, "TTP", vstart[7], step[7], min[7], max[7],ierflg); | |
868 | gMinuit->FixParameter(7); // 2.24 ... out of pulserRun Fit (->IRF) | |
869 | ||
870 | // Now ready for minimization step | |
871 | arglist[0] = 2000; // max num of iterations | |
872 | arglist[1] = 0.1; // tolerance | |
873 | ||
874 | gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg); | |
875 | ||
876 | Double_t p1 = 0.0 ; | |
877 | gMinuit->mnexcm("SET NOW", &p1 , 0, ierflg) ; // No Warnings | |
878 | ||
879 | if (ierflg == 4) { // Fit failed | |
880 | for (Int_t i=0;i<3;i++) { | |
881 | coefP[i] = 0; | |
882 | coefZ[i] = 0; | |
883 | } | |
884 | gMinuit->~TMinuit(); | |
885 | return 0; | |
886 | } else { // Fit successfull | |
887 | ||
888 | // Extract parameters from TMinuit | |
889 | Double_t *fitParam = new Double_t[6]; | |
890 | for (Int_t i=0;i<6;i++) { | |
891 | Double_t err = 0; | |
892 | Double_t val = 0; | |
893 | gMinuit->GetParameter(i,val,err); | |
894 | fitParam[i] = val; | |
895 | } | |
896 | ||
897 | // calculates the first 2 sets (4 PZ values) out of the fitted parameters | |
898 | Double_t *valuePZ = ExtractPZValues(fitParam); | |
899 | ||
900 | // TCF coefficients which are used for the equalisation step (stage) | |
901 | // ZERO/POLE Filter | |
9c2921ef | 902 | coefZ[0] = TMath::Exp(-1/valuePZ[2]); |
903 | coefZ[1] = TMath::Exp(-1/valuePZ[3]); | |
904 | coefP[0] = TMath::Exp(-1/valuePZ[0]); | |
905 | coefP[1] = TMath::Exp(-1/valuePZ[1]); | |
eb7e0771 | 906 | |
907 | fitParam->~Double_t(); | |
908 | valuePZ->~Double_t(); | |
909 | gMinuit->~TMinuit(); | |
910 | ||
911 | return 1; | |
912 | ||
913 | } | |
914 | ||
915 | } | |
916 | ||
917 | ||
918 | //____________________________________________________________________________ | |
919 | void AliTPCCalibTCF::FitFcn(Int_t &/*nPar*/, Double_t */*grad*/, Double_t &f, Double_t *par, Int_t /*iflag*/) | |
920 | { | |
921 | // | |
922 | // Minimization function needed for TMinuit with FitFunction included | |
923 | // Fit function: Sum of three convolution terms (IRF conv. with Exp.) | |
924 | // | |
925 | ||
926 | // Get Data ... | |
927 | TNtuple *dataTuple = (TNtuple *) gMinuit->GetObjectFit(); | |
928 | ||
929 | //calculate chisquare | |
930 | Double_t chisq = 0; | |
931 | Double_t delta = 0; | |
932 | for (Int_t i=0; i<dataTuple->GetEntries(); i++) { // loop over data points | |
933 | dataTuple->GetEntry(i); | |
934 | Float_t *p = dataTuple->GetArgs(); | |
935 | Double_t t = p[0]; | |
936 | Double_t signal = p[1]; // Normalized signal | |
937 | Double_t error = p[2]; | |
938 | ||
939 | // definition and evaluation if the IonTail specific fit function | |
940 | Double_t sigFit = 0; | |
941 | ||
942 | Double_t ttp = par[7]; // signal shaper raising time | |
943 | t=t-par[6]; // time adjustment | |
944 | ||
945 | if (t<0) { | |
946 | sigFit = 0; | |
947 | } else { | |
9c2921ef | 948 | 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 | 949 | |
9c2921ef | 950 | 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 | 951 | |
9c2921ef | 952 | 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 | 953 | |
954 | sigFit = par[0]*f1 + par[1]*f2 +par[2]*f3; | |
955 | } | |
956 | ||
957 | // chisqu calculation | |
958 | delta = (signal-sigFit)/error; | |
959 | chisq += delta*delta; | |
960 | } | |
961 | ||
962 | f = chisq; | |
963 | ||
964 | } | |
965 | ||
966 | ||
967 | ||
968 | //____________________________________________________________________________ | |
969 | Double_t* AliTPCCalibTCF::ExtractPZValues(Double_t *param) { | |
970 | // | |
971 | // Calculation of Pole and Zero values out of fit parameters | |
972 | // | |
973 | ||
974 | Double_t vA1, vA2, vA3, vTT1, vTT2, vTT3, vTa, vTb; | |
975 | vA1 = 0; vA2 = 0; vA3 = 0; | |
976 | vTT1 = 0; vTT2 = 0; vTT3 = 0; | |
977 | vTa = 0; vTb = 0; | |
978 | ||
979 | // nasty method of sorting the fit parameters to avoid wrong mapping | |
980 | // to the different stages of the TCF filter | |
981 | // (e.g. first 2 fit parameters represent the electron signal itself!) | |
982 | ||
983 | if (param[3]==param[4]) {param[3]=param[3]+0.0001;} | |
984 | if (param[5]==param[4]) {param[5]=param[5]+0.0001;} | |
985 | ||
986 | if ((param[5]>param[4])&(param[5]>param[3])) { | |
987 | if (param[4]>=param[3]) { | |
988 | vA1 = param[0]; vA2 = param[1]; vA3 = param[2]; | |
989 | vTT1 = param[3]; vTT2 = param[4]; vTT3 = param[5]; | |
990 | } else { | |
991 | vA1 = param[1]; vA2 = param[0]; vA3 = param[2]; | |
992 | vTT1 = param[4]; vTT2 = param[3]; vTT3 = param[5]; | |
993 | } | |
994 | } else if ((param[4]>param[5])&(param[4]>param[3])) { | |
995 | if (param[5]>=param[3]) { | |
996 | vA1 = param[0]; vA2 = param[2]; vA3 = param[1]; | |
997 | vTT1 = param[3]; vTT2 = param[5]; vTT3 = param[4]; | |
998 | } else { | |
999 | vA1 = param[2]; vA2 = param[0]; vA3 = param[1]; | |
1000 | vTT1 = param[5]; vTT2 = param[3]; vTT3 = param[4]; | |
1001 | } | |
1002 | } else if ((param[3]>param[4])&(param[3]>param[5])) { | |
1003 | if (param[5]>=param[4]) { | |
1004 | vA1 = param[1]; vA2 = param[2]; vA3 = param[0]; | |
1005 | vTT1 = param[4]; vTT2 = param[5]; vTT3 = param[3]; | |
1006 | } else { | |
1007 | vA1 = param[2]; vA2 = param[1]; vA3 = param[0]; | |
1008 | vTT1 = param[5]; vTT2 = param[4]; vTT3 = param[3]; | |
1009 | } | |
1010 | } | |
1011 | ||
1012 | ||
1013 | // Transformation of fit parameters into PZ values (needed by TCF) | |
1014 | Double_t beq = (vA1/vTT2+vA1/vTT3+vA2/vTT1+vA2/vTT3+vA3/vTT1+vA3/vTT2)/(vA1+vA2+vA3); | |
1015 | Double_t ceq = (vA1/(vTT2*vTT3)+vA2/(vTT1*vTT3)+vA3/(vTT1*vTT2))/(vA1+vA2+vA3); | |
1016 | ||
1017 | Double_t s1 = -beq/2-sqrt((beq*beq-4*ceq)/4); | |
1018 | Double_t s2 = -beq/2+sqrt((beq*beq-4*ceq)/4); | |
1019 | ||
1020 | if (vTT2<vTT3) {// not necessary but avoids significant undershots in first PZ | |
1021 | vTa = -1/s1; | |
1022 | vTb = -1/s2; | |
1023 | }else{ | |
1024 | vTa = -1/s2; | |
1025 | vTb = -1/s1; | |
1026 | } | |
1027 | ||
1028 | Double_t *valuePZ = new Double_t[4]; | |
1029 | valuePZ[0]=vTa; | |
1030 | valuePZ[1]=vTb; | |
1031 | valuePZ[2]=vTT2; | |
1032 | valuePZ[3]=vTT3; | |
1033 | ||
1034 | return valuePZ; | |
1035 | ||
1036 | } | |
1037 | ||
1038 | ||
1039 | //____________________________________________________________________________ | |
1040 | void AliTPCCalibTCF::Equalization(TNtuple *dataTuple, Double_t *coefZ, Double_t *coefP) { | |
1041 | // | |
1042 | // calculates the 3rd set of TCF parameters (remaining 2 PZ values) in | |
1043 | // order to restore the original pulse height and adds them to the passed arrays | |
1044 | // | |
1045 | ||
1046 | Double_t *s0 = new Double_t[1000]; // original pulse | |
1047 | Double_t *s1 = new Double_t[1000]; // pulse after 1st PZ filter | |
1048 | Double_t *s2 = new Double_t[1000]; // pulse after 2nd PZ filter | |
1049 | ||
1050 | const Int_t kPulseLength = dataTuple->GetEntries(); | |
1051 | ||
1052 | for (Int_t ipos=0; ipos<kPulseLength; ipos++) { | |
1053 | dataTuple->GetEntry(ipos); | |
1054 | Float_t *p = dataTuple->GetArgs(); | |
1055 | s0[ipos] = p[1]; | |
1056 | } | |
1057 | ||
1058 | // non-discret implementation of the first two TCF stages (recursive formula) | |
1059 | // discrete Altro emulator is not used because of accuracy! | |
1060 | s1[0] = s0[0]; // 1st PZ filter | |
1061 | for(Int_t ipos = 1; ipos < kPulseLength ; ipos++){ | |
1062 | s1[ipos] = s0[ipos] + coefP[0]*s1[ipos-1] - coefZ[0]*s0[ipos-1]; | |
1063 | } | |
1064 | s2[0] = s1[0]; // 2nd PZ filter | |
1065 | for(Int_t ipos = 1; ipos < kPulseLength ; ipos++){ | |
1066 | s2[ipos] = s1[ipos] + coefP[1]*s2[ipos-1] - coefZ[1]*s1[ipos-1]; | |
1067 | } | |
1068 | ||
1069 | // find maximum amplitude and position of original pulse and pulse after | |
1070 | // the first two stages of the TCF | |
1071 | Int_t s0pos = 0, s2pos = 0; | |
1072 | Double_t s0ampl = s0[0], s2ampl = s2[0]; // start values | |
1073 | for(Int_t ipos = 1; ipos < kPulseLength; ipos++){ | |
1074 | if (s0[ipos] > s0ampl){ | |
1075 | s0ampl = s0[ipos]; | |
1076 | s0pos = ipos; // should be pos 11 ... check? | |
1077 | } | |
1078 | if (s2[ipos] > s2ampl){ | |
1079 | s2ampl = s2[ipos]; | |
1080 | s2pos = ipos; | |
1081 | } | |
1082 | } | |
1083 | // calculation of 3rd set ... | |
1084 | if(s0ampl > s2ampl){ | |
1085 | coefZ[2] = 0; | |
1086 | coefP[2] = (s0ampl - s2ampl)/s0[s0pos-1]; | |
1087 | } else if (s0ampl < s2ampl) { | |
1088 | coefP[2] = 0; | |
1089 | coefZ[2] = (s2ampl - s0ampl)/s0[s0pos-1]; | |
1090 | } else { // same height ? will most likely not happen ? | |
1091 | coefP[2] = 0; | |
1092 | coefZ[2] = 0; | |
1093 | } | |
1094 | ||
1095 | s0->~Double_t(); | |
1096 | s1->~Double_t(); | |
1097 | s2->~Double_t(); | |
1098 | ||
1099 | } | |
1100 | ||
1101 | ||
1102 | ||
1103 | //____________________________________________________________________________ | |
1104 | Int_t AliTPCCalibTCF::FindCorTCFparam(TH1F *hisIn, const char *nameFileTCF, Double_t *coefZ, Double_t *coefP) { | |
1105 | // | |
1106 | // This function searches for the correct TCF parameters to the given | |
1107 | // histogram 'hisIn' within the file 'nameFileTCF' | |
1108 | // If no parameters for this pad (padinfo within the histogram!) where found | |
1109 | // the function returns 0 | |
1110 | ||
1111 | // Int_t numPulse = (Int_t)hisIn->GetBinContent(1); // number of pulses | |
1112 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
1113 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
1114 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
1115 | Int_t nPulse = 0; | |
1116 | ||
1117 | //-- searching for calculated TCF parameters for this pad/sector | |
1118 | TFile fileTCF(nameFileTCF,"READ"); | |
1119 | TNtuple *paramTuple = (TNtuple*)fileTCF.Get("TCFparam"); | |
1120 | ||
1121 | // create selection criteria to find the correct TCF params | |
1122 | char sel[100]; | |
1123 | if ( paramTuple->GetEntries("row==-1&&pad==-1") ) { | |
1124 | // parameters per SECTOR | |
1125 | sprintf(sel,"sec==%d&&row==-1&&pad==-1",sector); | |
1126 | } else { | |
1127 | // parameters per PAD | |
1128 | sprintf(sel,"sec==%d&&row==%d&&pad==%d",sector,row,pad); | |
1129 | } | |
1130 | ||
1131 | // list should contain just ONE entry! ... otherwise there is a mistake! | |
1132 | Long64_t entry = paramTuple->Draw(">>list",sel,"entrylist"); | |
1133 | TEntryList *list = (TEntryList*)gDirectory->Get("list"); | |
1134 | ||
1135 | if (entry) { // TCF set was found for this pad | |
1136 | Long64_t pos = list->GetEntry(0); | |
1137 | paramTuple->GetEntry(pos); // get specific TCF parameters | |
1138 | Float_t *p = paramTuple->GetArgs(); | |
1139 | // check ... | |
1140 | if(sector==p[0]) {printf("sector ok ... "); } | |
1141 | if(row==p[1]) {printf("row ok ... "); } | |
1142 | if(pad==p[2]) {printf("pad ok ... \n"); } | |
1143 | ||
1144 | // number of averaged pulses used to produce TCF params | |
1145 | nPulse = (Int_t)p[3]; | |
1146 | // TCF parameters | |
1147 | coefZ[0] = p[4]; coefP[0] = p[7]; | |
1148 | coefZ[1] = p[5]; coefP[1] = p[8]; | |
1149 | coefZ[2] = p[6]; coefP[2] = p[9]; | |
1150 | ||
1151 | } else { // no specific TCF parameters found for this pad | |
1152 | ||
1153 | printf("no specific TCF paramaters found for pad in ...\n"); | |
1154 | printf("in Sector %d | Row %d | Pad %d |\n", sector, row, pad); | |
1155 | nPulse = 0; | |
1156 | coefZ[0] = 0; coefP[0] = 0; | |
1157 | coefZ[1] = 0; coefP[1] = 0; | |
1158 | coefZ[2] = 0; coefP[2] = 0; | |
1159 | ||
1160 | } | |
1161 | ||
1162 | fileTCF.Close(); | |
1163 | ||
1164 | return nPulse; // number of averaged pulses for producing the TCF params | |
1165 | ||
1166 | } | |
1167 | ||
1168 | ||
1169 | //____________________________________________________________________________ | |
1170 | Double_t *AliTPCCalibTCF::GetQualityOfTCF(TH1F *hisIn, Double_t *coefZ, Double_t *coefP, Int_t plotFlag) { | |
1171 | // | |
1172 | // This function evaluates the quality parameters of the given TCF parameters | |
1173 | // tested on the passed pulse (hisIn) | |
1174 | // The quality parameters are stored in an array. They are ... | |
1175 | // height deviation [ADC] | |
1176 | // area reduction [percent] | |
1177 | // width reduction [percent] | |
1178 | // mean undershot [ADC] | |
1179 | // maximum of undershot after pulse [ADC] | |
1180 | // Pulse RMS [ADC] | |
1181 | ||
1182 | // perform ALTRO emulator | |
1183 | TNtuple *pulseTuple = ApplyTCFilter(hisIn, coefZ, coefP, plotFlag); | |
1184 | ||
1185 | printf("calculate quality val. for pulse in ... "); | |
1186 | printf(" Sector %d | Row %d | Pad %d |\n", (Int_t)hisIn->GetBinContent(2), (Int_t)hisIn->GetBinContent(3), (Int_t)hisIn->GetBinContent(4)); | |
1187 | ||
1188 | // Reasonable limit for the calculation of the quality values | |
1189 | Int_t binLimit = 80; | |
1190 | ||
1191 | // ============== Variable preparation | |
1192 | ||
1193 | // -- height difference in percent of orginal pulse | |
1194 | Double_t maxSig = pulseTuple->GetMaximum("sig"); | |
1195 | Double_t maxSigTCF = pulseTuple->GetMaximum("sigAfterTCF"); | |
1196 | // -- area reduction (above zero!) | |
1197 | Double_t area = 0; | |
1198 | Double_t areaTCF = 0; | |
1199 | // -- width reduction at certain ADC treshold | |
1200 | // TODO: set treshold at ZS treshold? (3 sigmas of noise?) | |
1201 | Int_t threshold = 3; // treshold in percent | |
1202 | Int_t threshADC = (Int_t)(maxSig/100*threshold); | |
1203 | Int_t startOfPulse = 0; Int_t startOfPulseTCF = 0; | |
1204 | Int_t posOfStart = 0; Int_t posOfStartTCF = 0; | |
1205 | Int_t widthFound = 0; Int_t widthFoundTCF = 0; | |
1206 | Int_t width = 0; Int_t widthTCF = 0; | |
1207 | // -- Calcluation of Undershot (mean of negavive signal after the first | |
1208 | // undershot) | |
1209 | Double_t undershotTCF = 0; | |
1210 | Double_t undershotStart = 0; | |
1211 | // -- Calcluation of Undershot (Sum of negative signal after the pulse) | |
1212 | Double_t maxUndershot = 0; | |
1213 | ||
1214 | ||
1215 | // === loop over timebins to calculate quality parameters | |
1216 | for (Int_t i=0; i<binLimit; i++) { | |
1217 | ||
1218 | // Read signal values | |
1219 | pulseTuple->GetEntry(i); | |
1220 | Float_t *p = pulseTuple->GetArgs(); | |
1221 | Double_t sig = p[1]; | |
1222 | Double_t sigTCF = p[2]; | |
1223 | ||
1224 | // calculation of area (above zero) | |
1225 | if (sig>0) {area += sig; } | |
1226 | if (sigTCF>0) {areaTCF += sigTCF; } | |
1227 | ||
1228 | ||
1229 | // Search for width at certain ADC treshold | |
1230 | // -- original signal | |
1231 | if (widthFound == 0) { | |
1232 | if( (sig > threshADC) && (startOfPulse == 0) ){ | |
1233 | startOfPulse = 1; | |
1234 | posOfStart = i; | |
1235 | } | |
1236 | if( (sig < threshADC) && (startOfPulse == 1) ){ | |
1237 | widthFound = 1; | |
1238 | width = i - posOfStart + 1; | |
1239 | } | |
1240 | } | |
1241 | // -- signal after TCF | |
1242 | if (widthFoundTCF == 0) { | |
1243 | if( (sigTCF > threshADC) && (startOfPulseTCF == 0) ){ | |
1244 | startOfPulseTCF = 1; | |
1245 | posOfStartTCF = i; | |
1246 | } | |
1247 | if( (sigTCF < threshADC) && (startOfPulseTCF == 1) ){ | |
1248 | widthFoundTCF = 1; | |
1249 | widthTCF = i -posOfStartTCF + 1; | |
1250 | } | |
1251 | ||
1252 | } | |
1253 | ||
1254 | // finds undershot start | |
1255 | if ( (widthFoundTCF==1) && (sigTCF<0) ) { | |
1256 | undershotStart = 1; | |
1257 | } | |
1258 | ||
1259 | // Calculation of undershot sum (after pulse) | |
1260 | if ( widthFoundTCF==1 ) { | |
1261 | undershotTCF += sigTCF; | |
1262 | } | |
1263 | ||
1264 | // Search for maximal undershot (is equal to minimum after the pulse) | |
1265 | if ( (undershotStart==1)&&(i<(posOfStartTCF+widthTCF+20)) ) { | |
1266 | if (maxUndershot>sigTCF) { maxUndershot = sigTCF; } | |
1267 | } | |
1268 | ||
1269 | } | |
1270 | ||
1271 | // == Calculation of Quality parameters | |
1272 | ||
1273 | // -- height difference in ADC | |
1274 | Double_t heightDev = maxSigTCF-maxSig; | |
1275 | ||
1276 | // Area reduction of the pulse in percent | |
1277 | Double_t areaReduct = 100-areaTCF/area*100; | |
1278 | ||
1279 | // Width reduction in percent | |
1280 | Double_t widthReduct = 0; | |
1281 | if ((widthFound==1)&&(widthFoundTCF==1)) { // in case of not too big IonTail | |
1282 | widthReduct = 100-(Double_t)widthTCF/(Double_t)width*100; | |
1283 | if (widthReduct<0) { widthReduct = 0;} | |
1284 | } | |
1285 | ||
1286 | // Undershot - mean of neg.signals after pulse | |
1287 | Double_t length = 1; | |
1288 | if (binLimit-widthTCF-posOfStartTCF) { length = (binLimit-widthTCF-posOfStartTCF);} | |
1289 | Double_t undershot = undershotTCF/length; | |
1290 | ||
1291 | ||
1292 | // calculation of pulse RMS with timebins before and at the end of the pulse | |
1293 | TH1I *tempRMSHis = new TH1I("tempRMSHis","tempRMSHis",100,-50,50); | |
1294 | for (Int_t ipos = 0; ipos<6; ipos++) { | |
1295 | // before the pulse | |
1296 | tempRMSHis->Fill(hisIn->GetBinContent(ipos+5)); | |
1297 | // at the end | |
1298 | tempRMSHis->Fill(hisIn->GetBinContent(hisIn->GetNbinsX()-ipos)); | |
1299 | } | |
1300 | Double_t pulseRMS = tempRMSHis->GetRMS(); | |
1301 | tempRMSHis->~TH1I(); | |
1302 | ||
1303 | if (plotFlag) { | |
1304 | // == Output | |
1305 | printf("height deviation [ADC]:\t\t\t %3.1f\n", heightDev); | |
1306 | printf("area reduction [percent]:\t\t %3.1f\n", areaReduct); | |
1307 | printf("width reduction [percent]:\t\t %3.1f\n", widthReduct); | |
1308 | printf("mean undershot [ADC]:\t\t\t %3.1f\n", undershot); | |
1309 | printf("maximum of undershot after pulse [ADC]: %3.1f\n", maxUndershot); | |
1310 | printf("RMS of the original pulse [ADC]: \t %3.2f\n\n", pulseRMS); | |
1311 | ||
1312 | } | |
1313 | ||
1314 | Double_t *qualityParam = new Double_t[6]; | |
1315 | qualityParam[0] = heightDev; | |
1316 | qualityParam[1] = areaReduct; | |
1317 | qualityParam[2] = widthReduct; | |
1318 | qualityParam[3] = undershot; | |
1319 | qualityParam[4] = maxUndershot; | |
1320 | qualityParam[5] = pulseRMS; | |
1321 | ||
1322 | pulseTuple->~TNtuple(); | |
1323 | ||
1324 | return qualityParam; | |
1325 | } | |
1326 | ||
1327 | ||
1328 | //____________________________________________________________________________ | |
1329 | TNtuple *AliTPCCalibTCF::ApplyTCFilter(TH1F *hisIn, Double_t *coefZ, Double_t *coefP, Int_t plotFlag) { | |
1330 | // | |
1331 | // Applies the given TCF parameters on the given pulse via the ALTRO emulator | |
1332 | // class (discret values) and stores both pulses into a returned TNtuple | |
1333 | // | |
1334 | ||
1335 | Int_t nbins = hisIn->GetNbinsX() -4; | |
1336 | // -1 because the first four timebins usually contain pad specific informations | |
1337 | Int_t nPulse = (Int_t)hisIn->GetBinContent(1); // Number of summed pulses | |
1338 | Int_t sector = (Int_t)hisIn->GetBinContent(2); | |
1339 | Int_t row = (Int_t)hisIn->GetBinContent(3); | |
1340 | Int_t pad = (Int_t)hisIn->GetBinContent(4); | |
1341 | ||
1342 | // redirect histogram values to arrays (discrete for altro emulator) | |
1343 | Double_t *signalIn = new Double_t[nbins]; | |
1344 | Double_t *signalOut = new Double_t[nbins]; | |
1345 | short *signalInD = new short[nbins]; | |
1346 | short *signalOutD = new short[nbins]; | |
1347 | for (Int_t ipos=0;ipos<nbins;ipos++) { | |
1348 | Double_t signal = hisIn->GetBinContent(ipos+5); // summed signal | |
1349 | signalIn[ipos]=signal/nPulse; // mean signal | |
1350 | signalInD[ipos]=(short)(TMath::Nint(signalIn[ipos])); //discrete mean signal | |
1351 | signalOutD[ipos]=signalInD[ipos]; // will be overwritten by AltroEmulator | |
1352 | } | |
1353 | ||
1354 | // transform TCF parameters into ALTRO readable format (Integer) | |
1355 | Int_t* valP = new Int_t[3]; | |
1356 | Int_t* valZ = new Int_t[3]; | |
1357 | for (Int_t i=0; i<3; i++) { | |
9c2921ef | 1358 | valP[i] = (Int_t)(coefP[i]*(TMath::Power(2,16)-1)); |
1359 | valZ[i] = (Int_t)(coefZ[i]*(TMath::Power(2,16)-1)); | |
eb7e0771 | 1360 | } |
1361 | ||
1362 | // discret ALTRO EMULATOR ____________________________ | |
1363 | AliTPCAltroEmulator *altro = new AliTPCAltroEmulator(nbins, signalOutD); | |
1364 | altro->ConfigAltro(0,1,0,0,0,0); // perform just the TailCancelation | |
1365 | altro->ConfigTailCancellationFilter(valP[0],valP[1],valP[2],valZ[0],valZ[1],valZ[2]); | |
1366 | altro->RunEmulation(); | |
1367 | delete altro; | |
1368 | ||
1369 | // non-discret implementation of the (recursive formula) | |
1370 | // discrete Altro emulator is not used because of accuracy! | |
1371 | Double_t *s1 = new Double_t[1000]; // pulse after 1st PZ filter | |
1372 | Double_t *s2 = new Double_t[1000]; // pulse after 2nd PZ filter | |
1373 | s1[0] = signalIn[0]; // 1st PZ filter | |
1374 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1375 | s1[ipos] = signalIn[ipos] + coefP[0]*s1[ipos-1] - coefZ[0]*signalIn[ipos-1]; | |
1376 | } | |
1377 | s2[0] = s1[0]; // 2nd PZ filter | |
1378 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1379 | s2[ipos] = s1[ipos] + coefP[1]*s2[ipos-1] - coefZ[1]*s1[ipos-1]; | |
1380 | } | |
1381 | signalOut[0] = s2[0]; // 3rd PZ filter | |
1382 | for(Int_t ipos = 1; ipos<nbins; ipos++){ | |
1383 | signalOut[ipos] = s2[ipos] + coefP[2]*signalOut[ipos-1] - coefZ[2]*s2[ipos-1]; | |
1384 | } | |
1385 | s1->~Double_t(); | |
1386 | s2->~Double_t(); | |
1387 | ||
1388 | // writing pulses to tuple | |
1389 | TNtuple *pulseTuple = new TNtuple("ntupleTCF","PulseTCF","timebin:sig:sigAfterTCF:sigND:sigNDAfterTCF"); | |
1390 | for (Int_t ipos=0;ipos<nbins;ipos++) { | |
1391 | pulseTuple->Fill(ipos,signalInD[ipos],signalOutD[ipos],signalIn[ipos],signalOut[ipos]); | |
1392 | } | |
1393 | ||
1394 | if (plotFlag) { | |
1395 | char hname[100]; | |
1396 | sprintf(hname,"sec%drow%dpad%d",sector,row,pad); | |
1397 | new TCanvas(hname,hname,600,400); | |
1398 | //just plotting non-discret pulses | they look pretties in case of mean sig ;-) | |
1399 | pulseTuple->Draw("sigND:timebin","","L"); | |
1400 | // pulseTuple->Draw("sig:timebin","","Lsame"); | |
1401 | pulseTuple->SetLineColor(3); | |
1402 | pulseTuple->Draw("sigNDAfterTCF:timebin","","Lsame"); | |
1403 | // pulseTuple->Draw("sigAfterTCF:timebin","","Lsame"); | |
1404 | } | |
1405 | ||
1406 | valP->~Int_t(); | |
1407 | valZ->~Int_t(); | |
1408 | ||
1409 | signalIn->~Double_t(); | |
1410 | signalOut->~Double_t(); | |
1411 | delete signalIn; | |
1412 | delete signalOut; | |
1413 | ||
1414 | return pulseTuple; | |
1415 | ||
1416 | } | |
1417 | ||
1418 | ||
1419 | ||
1420 | ||
1421 | //____________________________________________________________________________ | |
1422 | void AliTPCCalibTCF::PrintPulseThresholds() { | |
1423 | // | |
1424 | // Prints the pulse threshold settings | |
1425 | // | |
1426 | ||
1427 | printf(" %4.0d [ADC] ... expected Gate fluctuation length \n", fGateWidth); | |
1428 | printf(" %4.0d [ADC] ... expected usefull signal length \n", fSample); | |
1429 | printf(" %4.0d [ADC] ... needed pulselength for TC characterisation \n", fPulseLength); | |
1430 | printf(" %4.0d [ADC] ... lower pulse height limit \n", fLowPulseLim); | |
1431 | printf(" %4.0d [ADC] ... upper pulse height limit \n", fUpPulseLim); | |
1432 | printf(" %4.1f [ADC] ... maximal pulse RMS \n", fRMSLim); | |
1433 | ||
1434 | } |