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