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