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