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