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