TPC/Attic/AliTPCCalibTCF.cxx

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