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