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