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