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