TPC/Attic/AliTPCCalibTCF.cxx

   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 /// \class AliTPCCalibTCF
  18 ///
  19 /// Class for Evaluation and Validation of the ALTRO Tail Cancelation Filter
  20 /// (TCF) parameters out of TPC Raw data
  21 ///
  22 /// \author Stefan Rossegger, Simon Feigl
  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
  49 /// \cond CLASSIMP
  50 ClassImp(AliTPCCalibTCF)
  51 /// \endcond
  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 {
  80   ///  AliTPCCalibTCF constructor with specific (non-standard) thresholds
  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 {
  95   ///  AliTPCCalibTCF copy constructor
  96
  97 }
  98
  99
 100 //_____________________________________________________________________________
 101 AliTPCCalibTCF& AliTPCCalibTCF::operator = (const AliTPCCalibTCF &source)
 102 {
 103   /// AliTPCCalibTCF assignment operator
 104
 105   if (&source == this) return *this;
 106   new (this) AliTPCCalibTCF(source);
 107
 108   return *this;
 109
 110 }
 111
 112 //_____________________________________________________________________________
 113 AliTPCCalibTCF::~AliTPCCalibTCF()
 114 {
 115   /// AliTPCCalibTCF destructor
 116
 117 }
 118
 119
 120 //_____________________________________________________________________________
 121 void AliTPCCalibTCF::ProcessRawFileV3(const char *nameRawFile, const char *nameFileOut) {
 122   /// New RCU data format!: Standard middle of 2009
 123   ///
 124   /// Loops over all events within one RawData file and collects proper pulses
 125   /// (according to given tresholds) per pad
 126   /// Histograms per pad are stored in '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) {
 163   /// New RCU data format!: Standard middle of 2009
 164   ///
 165   /// Extracts proper pulses (according the given tresholds) within one event
 166   /// and accumulates them into one histogram per pad. All histograms are
 167   /// saved in the file 'nameFileOut'.
 168   /// The first bins of the histograms contain the following information:
 169   ///   bin 1: Number of accumulated pulses
 170   ///   bin 2;3;4: Sector; Row; Pad;
 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) {
 333   /// Merges all histograms within one sector, calculates the TCF parameters
 334   /// of the 'histogram-per-sector' and stores (histo and parameters) into
 335   /// seperated files ...
 336   ///
 337   /// note: first 4 timebins of a histogram hold specific informations
 338   ///       about number of collected pulses, sector, row and pad
 339   ///
 340   /// 'nameFileIn':  root file produced with Process function which holds
 341   ///                one histogram per pad (sum of signals of proper pulses)
 342   /// 'Sec+nameFileIn': root file with one histogram per sector
 343   ///                   (information of row and pad are set to -1)
 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) {
 410   /// This function takes a prepeared root file (accumulated histograms: output
 411   /// of process function) and performs an analysis (fit and equalization) in
 412   /// order to get the TCF parameters. These are stored in an TNtuple along with
 413   /// the pad and creation infos. The tuple is written to the output file
 414   /// "TCFparam+nameFileIn"
 415   /// To reduce the analysis time, the minimum number of accumulated pulses within
 416   /// one histogram 'minNumPulse' (to perform the analysis on) can be set
 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) {
 477   /// Performs the analysis on one specific pulse (histogram) by means of fitting
 478   /// the pulse and equalization of the pulseheight. The found TCF parameters
 479   /// are stored in the arrays coefZ and 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 {
 544   /// Performs quality parameters evaluation of the calculated TCF parameters in
 545   /// the file 'nameFileTCF' for every (accumulated) histogram within the
 546   /// prepeared root file 'nameFileIn'.
 547   /// The found quality parameters are stored in an TNtuple which will be saved
 548   /// in a Root file 'Quality-*'.
 549   /// If the parameter for the given pulse (given pad) was not found, the pulse
 550   /// is rejected.
 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) {
 615   /// Performs quality parameters evaluation of the calculated TCF parameters in
 616   /// the file 'nameFileTCF' for every proper pulse (according to given thresholds)
 617   /// within the RAW file 'nameRawFile'.
 618   /// The found quality parameters are stored in a TNtuple which will be saved
 619   /// in the Root file 'nameFileOut'. If the parameter for the given pulse
 620   /// (given pad) was not found, the pulse is rejected.
 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) {
 817   /// Plots the number of summed pulses per pad on a given TPC side
 818   /// 'nameFileIn': root-file created with the Process function
 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) {
 896   /// Plots the number of summed pulses per pad above a given minimum at the
 897   /// pad position at a given TPC side
 898   /// 'nameFile': root-file created with the Process function
 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 {
 978   /// This function is an easy interface to load the QualityTuple (produced with
 979   /// the function 'TestOn%File' and plots them according to the plot specifications
 980   /// 'plotSpec' e.g. "widthRed:maxUndershot"
 981   /// One may also set cut and plot options ("cut","pOpt")
 982   ///
 983   /// The stored quality parameters are ...
 984   ///   sec:row:pad:npulse: ... usual pad info
 985   ///   heightDev ... height deviation in percent
 986   ///   areaRed ... area reduction in percent
 987   ///   widthRed ... width reduction in percent
 988   ///   undershot ... mean undershot after the pulse in ADC
 989   ///   maxUndershot ... maximum of the undershot after the pulse in ADC
 990   ///   pulseRMS ... RMS of the pulse used to calculate the Quality parameters in ADC
 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) {
1021   /// function to fit one pulse and to calculate the according pole-zero parameters
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 {
1106   /// Minimization function needed for TMinuit with FitFunction included
1107   /// Fit function: Sum of three convolution terms (IRF conv. with Exp.)
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) {
1153   /// Calculation of Pole and Zero values out of fit parameters
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) {
1222   /// calculates the 3rd set of TCF parameters (remaining 2 PZ values) in
1223   /// order to restore the original pulse height and adds them to the passed arrays
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) {
1297   /// This function searches for the correct TCF parameters to the given
1298   /// histogram 'hisIn' within the file 'nameFileTCF'
1299   /// If no parameters for this pad (padinfo within the histogram!) where found
1300   /// the function returns 0
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) {
1362   /// This function evaluates the quality parameters of the given TCF parameters
1363   /// tested on the passed pulse (hisIn)
1364   /// The quality parameters are stored in an array. They are ...
1365   ///    height deviation [ADC]
1366   ///    area reduction [percent]
1367   ///    width reduction [percent]
1368   ///    mean undershot [ADC]
1369   ///    maximum of undershot after pulse [ADC]
1370   ///    Pulse RMS [ADC]
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) {
1520   /// Applies the given TCF parameters on the given pulse via the ALTRO emulator
1521   /// class (discret values) and stores both pulses into a returned TNtuple
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() {
1606   /// Prints the pulse threshold settings
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 {
1622   /// Gets histograms from fileNameIn and adds contents to fileSum
1623   ///
1624   /// If fileSum doesn't exist, fileSum is created
1625   ///   mode = 0, just ONE BIG FILE ('fileSum') will be used
1626   ///   mode = 1, one file per sector ('fileSum-Sec#.root') will be used
1627   /// mode=1 is much faster, but the additional function 'MergeToOneFile' has to be used in order to
1628   /// get one big and complete collection file again ...
1629   ///
1630   /// !Make sure not to add the same file more than once!
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
1706   /// Merges all Sec-files together ...
1707   /// this is an additional functionality for the function MergeHistsPerFile
1708   /// if for example mode=1
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) {
1762   /// Writes TCF parameters per PAD to .data file
1763   ///
1764   /// from now on: "roc" refers to the offline sector numbering
1765   ///              "sector" refers to the 18 sectors per side
1766   ///
1767   /// Gets TCF parameters of single pads from nameFileTCFPerPad and writes them to
1768   /// the file 'tpcTCFparamPAD.data'
1769   ///
1770   /// If there are parameters for a pad missing, then the parameters of the roc,
1771   /// in which the pad is located, are used as the pad parameters. The parameters for
1772   /// the roc are retreived from nameFileTCFPerSec. If there are parameters for
1773   /// a roc missing, then the parameters are set to -1.
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) {
1915   /// Writes TCF parameters per SECTOR (=ROC) to .data file
1916   ///
1917   /// from now on: "roc" refers to the offline sector numbering
1918   ///              "sector" refers to the 18 sectors per side
1919   ///
1920   /// Gets TCF parameters of a roc from nameFileTCFPerSec and writes them to
1921   /// the file 'tpcTCFparamSector.data'
1922   ///
1923   /// If there are parameters for a roc missing, then the parameters are set to -1
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 }