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