TPC/AliTPCCalibTCF.cxx

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