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