Classes moved to STEERBase.
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraFit.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////////////////////////
19 //                                                                             
20 // AliTRDCalibraFit                                                               
21 //                                                                             
22 // This class is for the TRD calibration of the relative gain factor, the drift velocity,
23 // the time 0 and the pad response function. It fits the histos.       
24 // The 2D histograms or vectors (first converted in 1D histos) will be fitted  
25 // if they have enough entries, otherwise the (default) value of the choosen database 
26 // will be put. For the relative gain calibration the resulted factors will be globally 
27 // normalized to the gain factors of the choosen database. It unables to precise     
28 // previous calibration procedure.
29 // The function SetDebug enables the user to see:                                     
30 // _fDebug = 0: nothing, only the values are written in the tree if wanted
31 // _fDebug = 1: a comparaison of the coefficients found and the default values 
32 //              in the choosen database.
33 //              fCoef , histogram of the coefs as function of the calibration group number
34 //              fDelta , histogram of the relative difference of the coef with the default
35 //                        value in the database as function of the calibration group number
36 //              fError , dirstribution of this relative difference
37 // _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
38 // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
39 //              pad row and col number
40 // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
41 //              also the comparaison histograms of the 1 for this detector
42 //
43 //                            
44 // Author:
45 //   R. Bailhache (R.Bailhache@gsi.de)
46 //                            
47 //////////////////////////////////////////////////////////////////////////////////////
48
49 #include <TTree.h>
50 #include <TLine.h>
51 #include <TH1I.h>
52 #include <TStyle.h>
53 #include <TProfile2D.h>
54 #include <TFile.h>
55 #include <TCanvas.h>
56 #include <TGraphErrors.h>
57 #include <TGraph.h>
58 #include <TObjArray.h>
59 #include <TH1.h>
60 #include <TH1F.h>
61 #include <TF1.h>
62 #include <TH2F.h>
63 #include <TAxis.h>
64 #include <TStopwatch.h>
65 #include <TMath.h>
66 #include <TLegend.h>
67 #include <TDirectory.h>
68 #include <TROOT.h>
69 #include <TTreeStream.h>
70 #include <TLinearFitter.h>
71 #include <TVectorD.h>
72 #include <TArrayF.h>
73
74 #include "AliLog.h"
75 #include "AliCDBManager.h"
76 #include "AliMathBase.h"
77
78 #include "AliTRDCalibraFit.h"
79 #include "AliTRDCalibraMode.h"
80 #include "AliTRDCalibraVector.h"
81 #include "AliTRDCalibraVdriftLinearFit.h"
82 #include "AliTRDcalibDB.h"
83 #include "AliTRDgeometry.h"
84 #include "AliTRDpadPlane.h"
85 #include "AliTRDgeometry.h"
86 #include "./Cal/AliTRDCalROC.h"
87 #include "./Cal/AliTRDCalPad.h"
88 #include "./Cal/AliTRDCalDet.h"
89
90
91 ClassImp(AliTRDCalibraFit)
92
93 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
94 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
95
96 //_____________singleton implementation_________________________________________________
97 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
98 {
99   //
100   // Singleton implementation
101   //
102
103   if (fgTerminated != kFALSE) {
104     return 0;
105   }
106
107   if (fgInstance == 0) {
108     fgInstance = new AliTRDCalibraFit();
109   }
110
111   return fgInstance;
112
113 }
114
115 //______________________________________________________________________________________
116 void AliTRDCalibraFit::Terminate()
117 {
118   //
119   // Singleton implementation
120   // Deletes the instance of this class
121   //
122
123   fgTerminated = kTRUE;
124
125   if (fgInstance != 0) {
126     delete fgInstance;
127     fgInstance = 0;
128   }
129
130 }
131
132 //______________________________________________________________________________________
133 AliTRDCalibraFit::AliTRDCalibraFit()
134   :TObject()
135   ,fGeo(0)
136   ,fNumberOfBinsExpected(0)
137   ,fMethod(0)
138   ,fBeginFitCharge(3.5)
139   ,fFitPHPeriode(1)
140   ,fTakeTheMaxPH(kFALSE)
141   ,fT0Shift(0.124797)
142   ,fRangeFitPRF(1.0)
143   ,fAccCDB(kFALSE)
144   ,fMinEntries(800)
145   ,fRebin(1)
146   ,fNumberFit(0)
147   ,fNumberFitSuccess(0)
148   ,fNumberEnt(0)
149   ,fStatisticMean(0.0)
150   ,fDebugStreamer(0x0)
151   ,fDebugLevel(0)
152   ,fFitVoir(0)
153   ,fMagneticField(0.5)
154   ,fCalibraMode(new AliTRDCalibraMode())
155   ,fCurrentCoefE(0.0)
156   ,fCurrentCoefE2(0.0)
157   ,fDect1(0)
158   ,fDect2(0)
159   ,fScaleFitFactor(0.0)
160   ,fEntriesCurrent(0)
161   ,fCountDet(0)
162   ,fCount(0)
163   ,fCalDet(0x0)
164   ,fCalROC(0x0)
165   ,fCalDet2(0x0)
166   ,fCalROC2(0x0)
167   ,fCurrentCoefDetector(0x0)
168   ,fCurrentCoefDetector2(0x0)
169   ,fVectorFit(0)
170   ,fVectorFit2(0)
171 {
172   //
173   // Default constructor
174   //
175
176   fGeo         = new AliTRDgeometry();  
177  
178   // Current variables initialised
179   for (Int_t k = 0; k < 2; k++) {
180     fCurrentCoef[k]  = 0.0;
181     fCurrentCoef2[k] = 0.0;
182   }
183   for (Int_t i = 0; i < 3; i++) {
184     fPhd[i]          = 0.0;
185     fDet[i]          = 0;
186   }
187  
188 }
189 //______________________________________________________________________________________
190 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
191 :TObject(c)
192 ,fGeo(0)
193 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
194 ,fMethod(c.fMethod)
195 ,fBeginFitCharge(c.fBeginFitCharge)
196 ,fFitPHPeriode(c.fFitPHPeriode)
197 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
198 ,fT0Shift(c.fT0Shift)
199 ,fRangeFitPRF(c.fRangeFitPRF)
200 ,fAccCDB(c.fAccCDB)
201 ,fMinEntries(c.fMinEntries)
202 ,fRebin(c.fRebin)
203 ,fNumberFit(c.fNumberFit)
204 ,fNumberFitSuccess(c.fNumberFitSuccess)
205 ,fNumberEnt(c.fNumberEnt)
206 ,fStatisticMean(c.fStatisticMean)
207 ,fDebugStreamer(0x0)
208 ,fDebugLevel(c.fDebugLevel)
209 ,fFitVoir(c.fFitVoir)
210 ,fMagneticField(c.fMagneticField)
211 ,fCalibraMode(0x0)
212 ,fCurrentCoefE(c.fCurrentCoefE)
213 ,fCurrentCoefE2(c.fCurrentCoefE2)
214 ,fDect1(c.fDect1)
215 ,fDect2(c.fDect2)
216 ,fScaleFitFactor(c.fScaleFitFactor)
217 ,fEntriesCurrent(c.fEntriesCurrent)
218 ,fCountDet(c.fCountDet)
219 ,fCount(c.fCount)
220 ,fCalDet(0x0)
221 ,fCalROC(0x0)
222 ,fCalDet2(0x0)
223 ,fCalROC2(0x0)
224 ,fCurrentCoefDetector(0x0)
225 ,fCurrentCoefDetector2(0x0)
226 ,fVectorFit(0)
227 ,fVectorFit2(0)
228 {
229   //
230   // Copy constructor
231   //
232
233   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
234  
235   //Current variables initialised
236   for (Int_t k = 0; k < 2; k++) {
237     fCurrentCoef[k]  = 0.0;
238     fCurrentCoef2[k] = 0.0;
239   }
240   for (Int_t i = 0; i < 3; i++) {
241     fPhd[i]          = 0.0;
242     fDet[i]          = 0;
243   }
244   if(c.fCalDet) fCalDet   = new AliTRDCalDet(*c.fCalDet);
245   if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
246
247   if(c.fCalROC) fCalROC   = new AliTRDCalROC(*c.fCalROC);
248   if(c.fCalROC2) fCalROC  = new AliTRDCalROC(*c.fCalROC2);
249
250   fVectorFit.SetName(c.fVectorFit.GetName());
251   for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
252     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
253     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
254     Int_t ntotal = 1;
255     if (GetChamber(detector) == 2) {
256       ntotal = 1728;
257     }
258     else {
259       ntotal = 2304;
260     }
261     Float_t *coef = new Float_t[ntotal];
262     for (Int_t i = 0; i < ntotal; i++) {
263       coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
264     }
265     fitInfo->SetCoef(coef);
266     fitInfo->SetDetector(detector);
267     fVectorFit.Add((TObject *) fitInfo);
268   }
269   fVectorFit.SetName(c.fVectorFit.GetName());
270   for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
271     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
272     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
273     Int_t ntotal = 1;
274     if (GetChamber(detector) == 2) {
275       ntotal = 1728;
276     }
277     else {
278       ntotal = 2304;
279     }
280     Float_t *coef = new Float_t[ntotal];
281     for (Int_t i = 0; i < ntotal; i++) {
282       coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
283     }
284     fitInfo->SetCoef(coef);
285     fitInfo->SetDetector(detector);
286     fVectorFit2.Add((TObject *) fitInfo);
287   }
288   if (fGeo) {
289     delete fGeo;
290   }
291   fGeo = new AliTRDgeometry();
292
293 }
294 //____________________________________________________________________________________
295 AliTRDCalibraFit::~AliTRDCalibraFit()
296 {
297   //
298   // AliTRDCalibraFit destructor
299   //
300   if ( fDebugStreamer ) delete fDebugStreamer;
301   if ( fCalDet )  delete fCalDet;
302   if ( fCalDet2 ) delete fCalDet2;
303   if ( fCalROC )  delete fCalROC;
304   if ( fCalROC2 ) delete fCalROC2; 
305   fVectorFit.Delete();
306   fVectorFit2.Delete();
307   if (fGeo) {
308     delete fGeo;
309   }
310
311 }
312 //_____________________________________________________________________________
313 void AliTRDCalibraFit::Destroy() 
314 {
315   //
316   // Delete instance 
317   //
318
319   if (fgInstance) {
320     delete fgInstance;
321     fgInstance = 0x0;
322   }
323
324 }
325 //____________Functions fit Online CH2d________________________________________
326 Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch)
327 {
328   //
329   // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
330   // calibration group normalized the resulted coefficients (to 1 normally)
331   //
332
333   // Set the calibration mode
334   const char *name = ch->GetTitle();
335   SetModeCalibration(name,0);
336
337   // Number of Ybins (detectors or groups of pads)
338   Int_t    nbins   = ch->GetNbinsX();// charge
339   Int_t    nybins  = ch->GetNbinsY();// groups number
340   if (!InitFit(nybins,0)) {
341     return kFALSE;
342   }
343   if (!InitFitCH()) {
344     return kFALSE;
345   }
346   fStatisticMean        = 0.0;
347   fNumberFit            = 0;
348   fNumberFitSuccess     = 0;
349   fNumberEnt            = 0;
350   // Init fCountDet and fCount
351   InitfCountDetAndfCount(0);
352   // Beginning of the loop betwwen dect1 and dect2
353   for (Int_t idect = fDect1; idect < fDect2; idect++) {
354     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
355     UpdatefCountDetAndfCount(idect,0);
356     ReconstructFitRowMinRowMax(idect, 0);
357     // Take the histo
358     TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
359     projch->SetDirectory(0);
360     // Number of entries for this calibration group
361     Double_t nentries = 0.0;
362     Double_t mean = 0.0;
363     for (Int_t k = 0; k < nbins; k++) {
364       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
365       nentries += ch->GetBinContent(binnb);
366       mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
367       projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
368     }
369     projch->SetEntries(nentries);
370     //printf("The number of entries for the group %d is %f\n",idect,nentries);
371     if (nentries > 0) {
372       fNumberEnt++;
373       mean /= nentries;
374     }
375     // Rebin and statistic stuff
376     if (fRebin > 1) {
377       projch = ReBin((TH1I *) projch);
378     }
379     // This detector has not enough statistics or was off
380     if (nentries <= fMinEntries) {
381       NotEnoughStatisticCH(idect);
382       if (fDebugLevel != 1) {
383         delete projch;
384       }
385       continue;
386     }
387     // Statistics of the group fitted
388     fStatisticMean += nentries;
389     fNumberFit++;
390     //Method choosen
391     switch(fMethod)
392       {
393       case 0: FitMeanW((TH1 *) projch, nentries); break;
394       case 1: FitMean((TH1 *) projch, nentries, mean); break;
395       case 2: FitCH((TH1 *) projch, mean); break;
396       case 3: FitBisCH((TH1 *) projch, mean); break;
397       default: return kFALSE;
398       }
399     // Fill Infos Fit
400     FillInfosFitCH(idect);
401     // Memory!!!
402     if (fDebugLevel != 1) {
403       delete projch;
404     }
405   } // Boucle object
406   // Normierungcharge
407   if (fDebugLevel != 1) {
408     NormierungCharge();
409   }
410   // Mean Statistic
411   if (fNumberFit > 0) {
412     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
413     fStatisticMean = fStatisticMean / fNumberFit;
414   }
415   else {
416     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
417   }
418   delete fDebugStreamer;
419   fDebugStreamer = 0x0;
420
421   return kTRUE;
422 }
423 //____________Functions fit Online CH2d________________________________________
424 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
425 {
426   //
427   // Reconstruct a 1D histo from the vectorCH for each calibration group,
428   // fit the histo, normalized the resulted coefficients (to 1 normally)
429   //
430
431   // Set the calibraMode
432   const char *name = calvect->GetNameCH();
433   SetModeCalibration(name,0);  
434
435   // Number of Xbins (detectors or groups of pads)
436   if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
437     return kFALSE;
438   }
439   if (!InitFitCH()) {
440     return kFALSE;
441   }
442   fStatisticMean        = 0.0;
443   fNumberFit            = 0;
444   fNumberFitSuccess     = 0;
445   fNumberEnt            = 0;
446   // Init fCountDet and fCount
447   InitfCountDetAndfCount(0);
448   // Beginning of the loop between dect1 and dect2
449   for (Int_t idect = fDect1; idect < fDect2; idect++) {
450     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
451     UpdatefCountDetAndfCount(idect,0);
452     ReconstructFitRowMinRowMax(idect,0);
453     // Take the histo
454     Double_t nentries = 0.0;
455     Double_t mean = 0.0;
456     TH1F *projch = 0x0;
457     Bool_t something = kTRUE;
458     if(!calvect->GetCHEntries(fCountDet)) something = kFALSE;
459     if(something){
460       TString name("CH");
461       name += idect;
462       projch  = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) name);
463       projch->SetDirectory(0);
464       for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
465         nentries += projch->GetBinContent(k+1);
466         mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
467       }
468       if (nentries > 0) {
469         fNumberEnt++;
470         mean /= nentries;
471       }
472       //printf("The number of entries for the group %d is %f\n",idect,nentries);
473       // Rebin
474       if (fRebin >  1) {
475         projch = ReBin((TH1F *) projch);
476       }
477     }
478     // This detector has not enough statistics or was not found in VectorCH
479     if (nentries <= fMinEntries) {
480       NotEnoughStatisticCH(idect);
481       if (fDebugLevel != 1) {
482         if(projch) delete projch;
483       }     
484       continue;
485     }
486     // Statistic of the histos fitted
487     fStatisticMean += nentries;
488     fNumberFit++;
489     //Method choosen
490     switch(fMethod)
491       {
492       case 0: FitMeanW((TH1 *) projch, nentries); break;
493       case 1: FitMean((TH1 *) projch, nentries, mean); break;
494       case 2: FitCH((TH1 *) projch, mean); break;
495       case 3: FitBisCH((TH1 *) projch, mean); break;
496       default: return kFALSE;
497       }
498     // Fill Infos Fit
499     FillInfosFitCH(idect); 
500     // Memory!!!
501     if (fDebugLevel != 1) {
502       delete projch;
503     }
504   } // Boucle object
505   // Normierungcharge
506   if (fDebugLevel != 1) {
507     NormierungCharge();
508   }
509   // Mean Statistics
510   if (fNumberFit > 0) {
511     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
512     fStatisticMean = fStatisticMean / fNumberFit;
513   }
514   else {
515     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
516   }
517   delete fDebugStreamer;
518   fDebugStreamer = 0x0;
519   return kTRUE;
520 }
521 //________________functions fit Online PH2d____________________________________
522 Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph)
523 {
524   //
525   // Take the 1D profiles (average pulse height), projections of the 2D PH
526   // on the Xaxis, for each calibration group
527   // Reconstruct a drift velocity
528   // A first calibration of T0 is also made  using the same method
529   //
530
531   // Set the calibration mode
532   const char *name = ph->GetTitle();
533   SetModeCalibration(name,1);
534   
535   // Number of Xbins (detectors or groups of pads)
536   Int_t    nbins   = ph->GetNbinsX();// time
537   Int_t    nybins  = ph->GetNbinsY();// calibration group
538   if (!InitFit(nybins,1)) {
539     return kFALSE;
540   }
541   if (!InitFitPH()) {
542     return kFALSE;
543   }
544   fStatisticMean        = 0.0;
545   fNumberFit            = 0;
546   fNumberFitSuccess     = 0;
547   fNumberEnt            = 0;
548   // Init fCountDet and fCount
549   InitfCountDetAndfCount(1);
550   // Beginning of the loop
551   for (Int_t idect = fDect1; idect < fDect2; idect++) {
552     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
553     UpdatefCountDetAndfCount(idect,1);
554     ReconstructFitRowMinRowMax(idect,1);
555     // Take the histo
556     TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
557     projph->SetDirectory(0); 
558     // Number of entries for this calibration group
559     Double_t nentries = 0;
560     for (Int_t k = 0; k < nbins; k++) {
561       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
562       nentries += ph->GetBinEntries(binnb);
563     }
564     if (nentries > 0) {
565       fNumberEnt++;
566     }  
567     //printf("The number of entries for the group %d is %f\n",idect,nentries);
568     // This detector has not enough statistics or was off
569     if (nentries  <= fMinEntries) {
570       //printf("Not enough statistic!\n");
571       NotEnoughStatisticPH(idect);     
572       if (fDebugLevel != 1) {
573         delete projph;
574       }
575       continue;
576     }
577     // Statistics of the histos fitted
578     fNumberFit++;
579     fStatisticMean += nentries;
580     // Calcul of "real" coef
581     CalculVdriftCoefMean();
582     CalculT0CoefMean();
583     //Method choosen
584     switch(fMethod)
585       {
586       case 0: FitLagrangePoly((TH1 *) projph); break;
587       case 1: FitPente((TH1 *) projph); break;
588       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
589       default: return kFALSE;
590       }
591     // Fill the tree if end of a detector or only the pointer to the branch!!!
592     FillInfosFitPH(idect);
593     // Memory!!!
594     if (fDebugLevel != 1) {
595       delete projph;
596     }
597   } // Boucle object
598   // Mean Statistic
599   if (fNumberFit > 0) {
600     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
601     fStatisticMean = fStatisticMean / fNumberFit;
602   }
603   else {
604     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
605   }
606   delete fDebugStreamer;
607   fDebugStreamer = 0x0;
608   return kTRUE;
609 }
610 //____________Functions fit Online PH2d________________________________________
611 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
612 {
613   //
614   // Reconstruct the average pulse height from the vectorPH for each
615   // calibration group
616   // Reconstruct a drift velocity
617   // A first calibration of T0 is also made  using the same method (slope method)
618   //
619
620   // Set the calibration mode
621   const char *name = calvect->GetNamePH();
622   SetModeCalibration(name,1);
623
624   // Number of Xbins (detectors or groups of pads)
625   if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
626     return kFALSE;
627   }
628   if (!InitFitPH()) {
629     return kFALSE;
630   }
631   fStatisticMean        = 0.0;
632   fNumberFit            = 0;
633   fNumberFitSuccess     = 0;
634   fNumberEnt            = 0;
635   // Init fCountDet and fCount
636   InitfCountDetAndfCount(1);
637   // Beginning of the loop
638   for (Int_t idect = fDect1; idect < fDect2; idect++) {
639     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
640     UpdatefCountDetAndfCount(idect,1);
641     ReconstructFitRowMinRowMax(idect,1);
642     // Take the histo
643     TH1F *projph = 0x0;
644     fEntriesCurrent = 0;
645     Bool_t something = kTRUE;
646     if(!calvect->GetPHEntries(fCountDet)) something = kFALSE;
647     if(something){
648       TString name("PH");
649       name += idect;
650       projph  = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
651       projph->SetDirectory(0);
652     }
653     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
654     // This detector has not enough statistics or was off
655     if (fEntriesCurrent <=  fMinEntries) {
656       //printf("Not enough stat!\n");
657       NotEnoughStatisticPH(idect);
658       if (fDebugLevel != 1) {
659         if(projph) delete projph;
660       }
661       continue;
662     }
663     // Statistic of the histos fitted
664     fNumberFit++;
665     fStatisticMean += fEntriesCurrent;
666     // Calcul of "real" coef
667     CalculVdriftCoefMean();
668     CalculT0CoefMean();
669     //Method choosen
670     switch(fMethod)
671       {
672       case 0: FitLagrangePoly((TH1 *) projph); break;
673       case 1: FitPente((TH1 *) projph); break;
674       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
675       default: return kFALSE;
676       }
677     // Fill the tree if end of a detector or only the pointer to the branch!!!
678     FillInfosFitPH(idect);
679     // Memory!!!
680     if (fDebugLevel != 1) {
681       delete projph;
682     }
683   } // Boucle object
684  
685   // Mean Statistic
686   if (fNumberFit > 0) {
687     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
688     fStatisticMean = fStatisticMean / fNumberFit;
689   }
690   else {
691     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
692   }
693   delete fDebugStreamer;
694   fDebugStreamer = 0x0;
695   return kTRUE;
696 }
697 //____________Functions fit Online PRF2d_______________________________________
698 Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf)
699 {
700   //
701   // Take the 1D profiles (pad response function), projections of the 2D PRF
702   // on the Xaxis, for each calibration group
703   // Fit with a gaussian to reconstruct the sigma of the pad response function
704   //
705
706   // Set the calibration mode
707   const char *name = prf->GetTitle();
708   SetModeCalibration(name,2);
709
710   // Number of Ybins (detectors or groups of pads)
711   Int_t    nybins  = prf->GetNbinsY();// calibration groups
712   Int_t    nbins   = prf->GetNbinsX();// bins
713   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
714   if((nbg > 0) || (nbg == -1)) return kFALSE;
715   if (!InitFit(nybins,2)) {
716     return kFALSE;
717   }
718   if (!InitFitPRF()) {
719     return kFALSE;
720   }
721   fStatisticMean        = 0.0;
722   fNumberFit            = 0;
723   fNumberFitSuccess     = 0;
724   fNumberEnt            = 0;
725   // Init fCountDet and fCount
726   InitfCountDetAndfCount(2);
727   // Beginning of the loop
728   for (Int_t idect = fDect1; idect < fDect2; idect++) {
729     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
730     UpdatefCountDetAndfCount(idect,2);
731     ReconstructFitRowMinRowMax(idect,2);
732     // Take the histo
733     TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
734     projprf->SetDirectory(0);
735     // Number of entries for this calibration group
736     Double_t nentries = 0;
737     for (Int_t k = 0; k < nbins; k++) {
738       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
739       nentries += prf->GetBinEntries(binnb);
740     }
741     if(nentries > 0) fNumberEnt++;
742     // This detector has not enough statistics or was off
743     if (nentries <= fMinEntries) {
744       NotEnoughStatisticPRF(idect);
745       if (fDebugLevel != 1) {
746         delete projprf;
747       }
748       continue;
749     }
750     // Statistics of the histos fitted
751     fNumberFit++;
752     fStatisticMean += nentries;
753     // Calcul of "real" coef
754     CalculPRFCoefMean();
755     //Method choosen
756     switch(fMethod)
757       {
758       case 0: FitPRF((TH1 *) projprf); break;
759       case 1: RmsPRF((TH1 *) projprf); break;
760       default: return kFALSE;
761       }
762     // Fill the tree if end of a detector or only the pointer to the branch!!!
763     FillInfosFitPRF(idect);
764     // Memory!!!
765     if (fDebugLevel != 1) {
766       delete projprf;
767     }
768   } // Boucle object
769   // Mean Statistic
770   if (fNumberFit > 0) {
771     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
772     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
773     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
774                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
775     fStatisticMean = fStatisticMean / fNumberFit;
776   }
777   else {
778     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
779   }
780   delete fDebugStreamer;
781   fDebugStreamer = 0x0;
782   return kTRUE;
783 }
784 //____________Functions fit Online PRF2d_______________________________________
785 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf)
786 {
787   //
788   // Take the 1D profiles (pad response function), projections of the 2D PRF
789   // on the Xaxis, for each calibration group
790   // Fit with a gaussian to reconstruct the sigma of the pad response function
791   //
792
793   // Set the calibration mode
794   const char *name = prf->GetTitle();
795   SetModeCalibration(name,2);
796
797   // Number of Ybins (detectors or groups of pads)
798   TAxis   *xprf    = prf->GetXaxis();
799   TAxis   *yprf    = prf->GetYaxis();
800   Int_t    nybins  = yprf->GetNbins();// calibration groups
801   Int_t    nbins   = xprf->GetNbins();// bins
802   Float_t  lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
803   Float_t  upedge  = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
804   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
805   if(nbg == -1) return kFALSE;
806   if(nbg > 0) fMethod = 1;
807   else fMethod = 0;
808   if (!InitFit(nybins,2)) {
809     return kFALSE;
810   }
811   if (!InitFitPRF()) {
812     return kFALSE;
813   }
814   fStatisticMean        = 0.0;
815   fNumberFit            = 0;
816   fNumberFitSuccess     = 0;
817   fNumberEnt            = 0;
818   // Init fCountDet and fCount
819   InitfCountDetAndfCount(2);
820   // Beginning of the loop
821   for (Int_t idect = fDect1; idect < fDect2; idect++) {
822     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
823     UpdatefCountDetAndfCount(idect,2);
824     ReconstructFitRowMinRowMax(idect,2);
825     // Build the array of entries and sum
826     TArrayD arraye  = TArrayD(nbins);
827     TArrayD arraym  = TArrayD(nbins);
828     TArrayD arrayme = TArrayD(nbins);
829     Double_t nentries = 0;
830     //printf("nbins %d\n",nbins);
831     for (Int_t k = 0; k < nbins; k++) {
832       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
833       Double_t entries = (Double_t)prf->GetBinEntries(binnb);
834       Double_t mean    = (Double_t)prf->GetBinContent(binnb);
835       Double_t error   = (Double_t)prf->GetBinError(binnb); 
836       //printf("for %d we have %f\n",k,entries);
837       nentries += entries;
838       arraye.AddAt(entries,k);
839       arraym.AddAt(mean,k);
840       arrayme.AddAt(error,k);
841     }
842     if(nentries > 0) fNumberEnt++;
843     //printf("The number of entries for the group %d is %f\n",idect,nentries);
844     // This detector has not enough statistics or was off
845     if (nentries <= fMinEntries) {
846       NotEnoughStatisticPRF(idect);
847       continue;
848     }
849     // Statistics of the histos fitted
850     fNumberFit++;
851     fStatisticMean += nentries;
852     // Calcul of "real" coef
853     CalculPRFCoefMean();
854     //Method choosen
855     switch(fMethod)
856       {
857       case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
858       case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
859       default: return kFALSE;
860       }
861     // Fill the tree if end of a detector or only the pointer to the branch!!!
862     FillInfosFitPRF(idect);
863   } // Boucle object
864   // Mean Statistic
865   if (fNumberFit > 0) {
866     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
867     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
868     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
869                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
870     fStatisticMean = fStatisticMean / fNumberFit;
871   }
872   else {
873     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
874   }
875   delete fDebugStreamer;
876   fDebugStreamer = 0x0;
877   return kTRUE;
878 }
879 //____________Functions fit Online PRF2d_______________________________________
880 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
881 {
882   //
883   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
884   // each calibration group
885   // Fit with a gaussian to reconstruct the sigma of the pad response function
886   //
887
888   // Set the calibra mode
889   const char *name = calvect->GetNamePRF();
890   SetModeCalibration(name,2);
891   //printf("test0 %s\n",name);
892
893   // Number of Xbins (detectors or groups of pads)
894   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
895     //printf("test1\n");
896     return kFALSE;
897   }
898   if (!InitFitPRF()) {
899     ///printf("test2\n");
900     return kFALSE;
901   }
902   fStatisticMean        = 0.0;
903   fNumberFit            = 0;
904   fNumberFitSuccess     = 0;
905   fNumberEnt            = 0;
906   // Init fCountDet and fCount
907   InitfCountDetAndfCount(2);
908   // Beginning of the loop
909   for (Int_t idect = fDect1; idect < fDect2; idect++) {
910     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
911     UpdatefCountDetAndfCount(idect,2);
912     ReconstructFitRowMinRowMax(idect,2);
913     // Take the histo
914     TH1F *projprf = 0x0;
915     fEntriesCurrent = 0;
916     Bool_t something = kTRUE;
917     if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
918     if(something){
919       TString name("PRF");
920       name += idect;
921       projprf  = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name)));
922       projprf->SetDirectory(0);
923     }
924     // This detector has not enough statistics or was off
925     if (fEntriesCurrent <= fMinEntries) {
926       NotEnoughStatisticPRF(idect);
927       if (fDebugLevel != 1) {
928         if(projprf) delete projprf;
929       }
930       continue;
931     }
932     // Statistic of the histos fitted
933     fNumberFit++;
934     fStatisticMean += fEntriesCurrent;
935     // Calcul of "real" coef
936     CalculPRFCoefMean();
937     //Method choosen
938     switch(fMethod)
939       {
940       case 1: FitPRF((TH1 *) projprf); break;
941       case 2: RmsPRF((TH1 *) projprf); break;
942       default: return kFALSE;
943       }    
944     // Fill the tree if end of a detector or only the pointer to the branch!!!
945     FillInfosFitPRF(idect);
946     // Memory!!!
947     if (fDebugLevel != 1) {
948       delete projprf;
949     }
950   } // Boucle object
951   // Mean Statistics
952   if (fNumberFit > 0) {
953     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
954   }
955   else {
956     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
957   }
958   delete fDebugStreamer;
959   fDebugStreamer = 0x0;
960   return kTRUE;
961 }
962 //____________Functions fit Online PRF2d_______________________________________
963 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
964 {
965   //
966   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
967   // each calibration group
968   // Fit with a gaussian to reconstruct the sigma of the pad response function
969   //
970
971   // Set the calibra mode
972   const char *name = calvect->GetNamePRF();
973   SetModeCalibration(name,2);
974   //printf("test0 %s\n",name);
975   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
976   printf("test1 %d\n",nbg);
977   if(nbg == -1) return kFALSE;
978   if(nbg > 0) fMethod = 1;
979   else fMethod = 0;
980   // Number of Xbins (detectors or groups of pads)
981   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
982     //printf("test2\n");
983     return kFALSE;
984   }
985   if (!InitFitPRF()) {
986     //printf("test3\n");
987     return kFALSE;
988   }
989   fStatisticMean        = 0.0;
990   fNumberFit            = 0;
991   fNumberFitSuccess     = 0;
992   fNumberEnt            = 0;
993   // Variables
994   Int_t nbins           = 0;
995   Double_t *arrayx       = 0;
996   Double_t *arraye       = 0;
997   Double_t *arraym       = 0;
998   Double_t *arrayme      = 0;
999   Float_t lowedge       = 0.0;
1000   Float_t upedge        = 0.0;
1001   // Init fCountDet and fCount
1002   InitfCountDetAndfCount(2);
1003   // Beginning of the loop
1004   for (Int_t idect = fDect1; idect < fDect2; idect++) {
1005     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1006     UpdatefCountDetAndfCount(idect,2);
1007     ReconstructFitRowMinRowMax(idect,2);
1008     // Take the histo
1009     TGraphErrors *projprftree = 0x0;
1010     fEntriesCurrent  = 0;
1011     Bool_t something = kTRUE;
1012     if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE;
1013     if(something){
1014       TString name("PRF");
1015       name += idect;
1016       projprftree  = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) name);
1017       nbins   = projprftree->GetN();
1018       arrayx  = (Double_t *)projprftree->GetX();
1019       arraye  = (Double_t *)projprftree->GetEX();
1020       arraym  = (Double_t *)projprftree->GetY();
1021       arrayme = (Double_t *)projprftree->GetEY();
1022       Float_t step = arrayx[1]-arrayx[0];
1023       lowedge = arrayx[0] - step/2.0;
1024       upedge  = arrayx[(nbins-1)] + step/2.0;
1025       //printf("nbins est %d\n",nbins);
1026       for(Int_t k = 0; k < nbins; k++){
1027         fEntriesCurrent += (Int_t)arraye[k];
1028         //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1029         if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1030       }
1031       if(fEntriesCurrent > 0) fNumberEnt++;
1032     }
1033     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1034     // This detector has not enough statistics or was off
1035     if (fEntriesCurrent <= fMinEntries) {
1036       NotEnoughStatisticPRF(idect);
1037       if(projprftree) delete projprftree;
1038       continue;
1039     }
1040     // Statistic of the histos fitted
1041     fNumberFit++;
1042     fStatisticMean += fEntriesCurrent;
1043     // Calcul of "real" coef
1044     CalculPRFCoefMean();
1045     //Method choosen
1046     switch(fMethod)
1047       {
1048       case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1049       case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1050       default: return kFALSE;
1051       }    
1052     // Fill the tree if end of a detector or only the pointer to the branch!!!
1053     FillInfosFitPRF(idect);
1054     // Memory!!!
1055     if (fDebugLevel != 1) {
1056       delete projprftree;
1057     }
1058   } // Boucle object
1059   // Mean Statistics
1060   if (fNumberFit > 0) {
1061     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1062   }
1063   else {
1064     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1065   }
1066   delete fDebugStreamer;
1067   fDebugStreamer = 0x0;
1068   return kTRUE;
1069 }
1070 //____________Functions fit Online CH2d________________________________________
1071 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1072 {
1073   //
1074   // The linear method
1075   //
1076
1077   fStatisticMean        = 0.0;
1078   fNumberFit            = 0;
1079   fNumberFitSuccess     = 0;
1080   fNumberEnt            = 0;
1081   if(!InitFitLinearFitter()) return kFALSE;
1082
1083   
1084   for(Int_t idet = 0; idet < 540; idet++){
1085
1086
1087     //printf("detector number %d\n",idet);
1088
1089     // Take the result
1090     TVectorD param(2);
1091     TVectorD error(3);
1092     fEntriesCurrent = 0;
1093     fCountDet       = idet;
1094     Bool_t here     = calivdli->GetParam(idet,&param);
1095     Bool_t heree    = calivdli->GetError(idet,&error);
1096     //printf("here %d and heree %d\n",here, heree);
1097     if(heree) {
1098       fEntriesCurrent = (Int_t) error[2];
1099       fNumberEnt++;
1100     }
1101     //printf("Number of entries %d\n",fEntriesCurrent);
1102     // Nothing found or not enough statistic
1103     if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1104       NotEnoughStatisticLinearFitter();
1105       continue;
1106     }
1107     //param.Print();
1108     //error.Print();
1109     //Statistics
1110     fNumberFit++;
1111     fStatisticMean += fEntriesCurrent;     
1112
1113     // Check the fit
1114     if((-(param[1])) <= 0.0) {
1115       NotEnoughStatisticLinearFitter();
1116       continue;
1117     }
1118
1119     // CalculDatabaseVdriftandTan
1120     CalculVdriftLorentzCoef();
1121
1122     // Statistics   
1123     fNumberFitSuccess ++;
1124
1125     // Put the fCurrentCoef
1126     fCurrentCoef[0]  = -param[1];
1127     // here the database must be the one of the reconstruction for the lorentz angle....
1128     fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1129     fCurrentCoefE    = error[1];
1130     fCurrentCoefE2   = error[0];
1131     if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1132       fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1133     }    
1134
1135     // Fill
1136     FillInfosFitLinearFitter();
1137
1138     
1139   }
1140   // Mean Statistics
1141   if (fNumberFit > 0) {
1142     AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
1143   }
1144   else {
1145     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1146   }
1147   delete fDebugStreamer;
1148   fDebugStreamer = 0x0;
1149   return kTRUE;
1150   
1151 }
1152 //____________Functions for seeing if the pad is really okey___________________
1153 //_____________________________________________________________________________
1154 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1155 {
1156   //
1157   // Get numberofgroupsprf
1158   //
1159   
1160   // Some patterns
1161   const Char_t *pattern0 = "Ngp0";
1162   const Char_t *pattern1 = "Ngp1";
1163   const Char_t *pattern2 = "Ngp2";
1164   const Char_t *pattern3 = "Ngp3";
1165   const Char_t *pattern4 = "Ngp4";
1166   const Char_t *pattern5 = "Ngp5";
1167   const Char_t *pattern6 = "Ngp6";
1168
1169   // Nrphi mode
1170   if (strstr(nametitle,pattern0)) {
1171     return 0;
1172   }
1173   if (strstr(nametitle,pattern1)) {
1174     return 1;
1175   }
1176   if (strstr(nametitle,pattern2)) {
1177     return 2;
1178   }
1179   if (strstr(nametitle,pattern3)) {
1180     return 3;
1181   }
1182   if (strstr(nametitle,pattern4)) {
1183     return 4;
1184   }
1185   if (strstr(nametitle,pattern5)) {
1186     return 5;
1187   }
1188   if (strstr(nametitle,pattern6)){
1189     return 6;
1190   }
1191   else return -1;
1192  
1193
1194 }
1195 //_____________________________________________________________________________
1196 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1197 {
1198   //
1199   // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1200   // corresponding to the given name
1201   //
1202
1203   if(!SetNzFromTObject(name,i)) return kFALSE;
1204   if(!SetNrphiFromTObject(name,i)) return kFALSE;
1205   
1206   return kTRUE; 
1207
1208 }
1209 //_____________________________________________________________________________
1210 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1211 {
1212   //
1213   // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1214   // corresponding to the given TObject
1215   //
1216   
1217   // Some patterns
1218   const Char_t *patternrphi0 = "Nrphi0";
1219   const Char_t *patternrphi1 = "Nrphi1";
1220   const Char_t *patternrphi2 = "Nrphi2";
1221   const Char_t *patternrphi3 = "Nrphi3";
1222   const Char_t *patternrphi4 = "Nrphi4";
1223   const Char_t *patternrphi5 = "Nrphi5";
1224   const Char_t *patternrphi6 = "Nrphi6";
1225
1226   // Nrphi mode
1227   if (strstr(name,patternrphi0)) {
1228     fCalibraMode->SetNrphi(i ,0);
1229     return kTRUE;
1230   }
1231   if (strstr(name,patternrphi1)) {
1232     fCalibraMode->SetNrphi(i, 1);
1233     return kTRUE;
1234   }
1235   if (strstr(name,patternrphi2)) {
1236     fCalibraMode->SetNrphi(i, 2);
1237     return kTRUE;
1238   }
1239   if (strstr(name,patternrphi3)) {
1240     fCalibraMode->SetNrphi(i, 3);
1241     return kTRUE;
1242   }
1243   if (strstr(name,patternrphi4)) {
1244     fCalibraMode->SetNrphi(i, 4);
1245     return kTRUE;
1246   }
1247   if (strstr(name,patternrphi5)) {
1248     fCalibraMode->SetNrphi(i, 5);
1249     return kTRUE;
1250   }
1251   if (strstr(name,patternrphi6)) {
1252     fCalibraMode->SetNrphi(i, 6);
1253     return kTRUE;
1254   }
1255   
1256   fCalibraMode->SetNrphi(i ,0);
1257   return kFALSE;
1258     
1259 }
1260 //_____________________________________________________________________________
1261 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1262 {
1263   //
1264   // Set fNz[i] of the AliTRDCalibraFit::Instance()
1265   // corresponding to the given TObject
1266   //
1267
1268   // Some patterns
1269   const Char_t *patternz0    = "Nz0";
1270   const Char_t *patternz1    = "Nz1";
1271   const Char_t *patternz2    = "Nz2";
1272   const Char_t *patternz3    = "Nz3";
1273   const Char_t *patternz4    = "Nz4";
1274   
1275   if (strstr(name,patternz0)) {
1276     fCalibraMode->SetNz(i, 0);
1277     return kTRUE;
1278   }
1279   if (strstr(name,patternz1)) {
1280     fCalibraMode->SetNz(i ,1);
1281     return kTRUE;
1282   }
1283   if (strstr(name,patternz2)) {
1284     fCalibraMode->SetNz(i ,2);
1285     return kTRUE;
1286   }
1287   if (strstr(name,patternz3)) {
1288     fCalibraMode->SetNz(i ,3);
1289     return kTRUE;  
1290   }
1291   if (strstr(name,patternz4)) {
1292     fCalibraMode->SetNz(i ,4);
1293     return kTRUE;
1294   }
1295
1296   fCalibraMode->SetNz(i ,0);
1297   return kFALSE;
1298 }
1299 //_____________________________________________________________________________
1300 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool_t perdetector)
1301 {
1302   //
1303   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1304   // It takes the mean value of the coefficients per detector 
1305   // This object has to be written in the database
1306   //
1307   
1308   // Create the DetObject
1309   AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1310
1311   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1312   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1313   Int_t detector = -1;
1314   Float_t value  = 0.0;
1315
1316   for (Int_t k = 0; k < loop; k++) {
1317     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1318     Float_t mean  = 0.0;
1319     if(perdetector){
1320       mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1321     }
1322     else {
1323       Int_t   count = 0;
1324       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1325       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1326       for (Int_t row = 0; row < rowMax; row++) {
1327         for (Int_t col = 0; col < colMax; col++) {
1328           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1329           mean += TMath::Abs(value);
1330           count++;       
1331         } // Col
1332       } // Row
1333       if(count > 0) mean = mean/count;
1334     }
1335     object->SetValue(detector,mean);
1336   }
1337   
1338   return object;
1339 }
1340 //_____________________________________________________________________________
1341 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector)
1342 {
1343   //
1344   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1345   // It takes the mean value of the coefficients per detector 
1346   // This object has to be written in the database
1347   //
1348   
1349   // Create the DetObject
1350   AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1351   
1352  
1353   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1354   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1355   Int_t detector = -1;
1356   Float_t value  = 0.0;
1357
1358   for (Int_t k = 0; k < loop; k++) {
1359     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();  
1360     Float_t mean  = 0.0;
1361     if(perdetector){
1362       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1363       if(value > 0) value = value*scaleFitFactor;
1364       mean = TMath::Abs(value);
1365     }
1366     else{
1367       Int_t   count = 0;
1368       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1369       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1370       for (Int_t row = 0; row < rowMax; row++) {
1371         for (Int_t col = 0; col < colMax; col++) {
1372           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1373           if(value > 0) value = value*scaleFitFactor;
1374           mean += TMath::Abs(value);
1375           count++;       
1376         } // Col
1377       } // Row
1378       if(count > 0) mean = mean/count;
1379     }
1380     object->SetValue(detector,mean);
1381   }
1382  
1383   return object;
1384 }
1385 //_____________________________________________________________________________
1386 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector)
1387 {
1388   //
1389   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1390   // It takes the min value of the coefficients per detector 
1391   // This object has to be written in the database
1392   //
1393   
1394   // Create the DetObject
1395   AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1396   
1397   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1398   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1399   Int_t detector = -1;
1400   Float_t value  = 0.0;
1401
1402   for (Int_t k = 0; k < loop; k++) {
1403     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();   
1404     Float_t min  = 100.0;
1405     if(perdetector){
1406       min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1407     }
1408     else{
1409       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1410       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1411       for (Int_t row = 0; row < rowMax; row++) {
1412         for (Int_t col = 0; col < colMax; col++) {
1413           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1414           if(min > value) min = value;
1415         } // Col
1416       } // Row
1417     }
1418     object->SetValue(detector,min);
1419   }
1420
1421   return object;
1422
1423 }
1424 //_____________________________________________________________________________
1425 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit)
1426 {
1427   //
1428   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1429   // It takes the min value of the coefficients per detector 
1430   // This object has to be written in the database
1431   //
1432   
1433   // Create the DetObject
1434   AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1435   
1436   
1437   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1438   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1439   Int_t detector = -1;
1440   Float_t value  = 0.0;
1441
1442   for (Int_t k = 0; k < loop; k++) {
1443     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1444     /*
1445       Int_t rowMax    = fGeo->GetRowMax(GetPlane(detector),GetChamber(detector),GetSector(detector));
1446       Int_t colMax    = fGeo->GetColMax(GetPlane(detector));
1447       Float_t min  = 100.0;
1448       for (Int_t row = 0; row < rowMax; row++) {
1449       for (Int_t col = 0; col < colMax; col++) {
1450       value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1451       mean += -TMath::Abs(value);
1452       count++;       
1453       } // Col
1454       } // Row
1455       if(count > 0) mean = mean/count;
1456     */
1457     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1458     object->SetValue(detector,-TMath::Abs(value));
1459   }
1460
1461   return object;
1462   
1463 }
1464 //_____________________________________________________________________________
1465 TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject)
1466 {
1467   //
1468   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1469   // You need first to create the object for the detectors,
1470   // where the mean value is put.
1471   // This object has to be written in the database
1472   //
1473   
1474   // Create the DetObject
1475   AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1476   
1477   if(!vectorFit){
1478     for(Int_t k = 0; k < 540; k++){
1479       AliTRDCalROC *calROC = object->GetCalROC(k);
1480       Int_t nchannels = calROC->GetNchannels();
1481       for(Int_t ch = 0; ch < nchannels; ch++){
1482         calROC->SetValue(ch,1.0);
1483       }
1484     }
1485   }
1486   else{
1487
1488     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1489     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1490     Int_t detector = -1;
1491     Float_t value  = 0.0;
1492     
1493     for (Int_t k = 0; k < loop; k++) {
1494       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1495       AliTRDCalROC *calROC = object->GetCalROC(detector);
1496       Float_t mean         = detobject->GetValue(detector);
1497       if(mean == 0) continue;
1498       Int_t rowMax    = calROC->GetNrows();
1499       Int_t colMax    = calROC->GetNcols();
1500       for (Int_t row = 0; row < rowMax; row++) {
1501         for (Int_t col = 0; col < colMax; col++) {
1502           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1503           if(value > 0) value = value*scaleFitFactor;
1504           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1505         } // Col
1506       } // Row
1507     } 
1508   }
1509
1510   return object;  
1511 }
1512 //_____________________________________________________________________________
1513 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject)
1514 {
1515   //
1516   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1517   // You need first to create the object for the detectors,
1518   // where the mean value is put.
1519   // This object has to be written in the database
1520   //
1521
1522   // Create the DetObject
1523   AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1524
1525   if(!vectorFit){
1526     for(Int_t k = 0; k < 540; k++){
1527       AliTRDCalROC *calROC = object->GetCalROC(k);
1528       Int_t nchannels = calROC->GetNchannels();
1529       for(Int_t ch = 0; ch < nchannels; ch++){
1530         calROC->SetValue(ch,1.0);
1531       }
1532     }
1533   }
1534   else {
1535     
1536     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1537     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1538     Int_t detector = -1;
1539     Float_t value  = 0.0;
1540     
1541     for (Int_t k = 0; k < loop; k++) {
1542       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1543       AliTRDCalROC *calROC = object->GetCalROC(detector);
1544       Float_t mean         = detobject->GetValue(detector);
1545       if(mean == 0) continue;
1546       Int_t rowMax    = calROC->GetNrows();
1547       Int_t colMax    = calROC->GetNcols();
1548       for (Int_t row = 0; row < rowMax; row++) {
1549         for (Int_t col = 0; col < colMax; col++) {
1550           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1551           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1552         } // Col
1553       } // Row
1554     } 
1555   }
1556   return object;    
1557
1558 }
1559 //_____________________________________________________________________________
1560 TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject)
1561 {
1562   //
1563   // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1564   // You need first to create the object for the detectors,
1565   // where the mean value is put.
1566   // This object has to be written in the database
1567   //
1568   
1569   // Create the DetObject
1570   AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1571
1572   if(!vectorFit){
1573     for(Int_t k = 0; k < 540; k++){
1574       AliTRDCalROC *calROC = object->GetCalROC(k);
1575       Int_t nchannels = calROC->GetNchannels();
1576       for(Int_t ch = 0; ch < nchannels; ch++){
1577         calROC->SetValue(ch,0.0);
1578       }
1579     }
1580   }
1581   else {
1582     
1583     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1584     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1585     Int_t detector = -1;
1586     Float_t value  = 0.0;
1587     
1588     for (Int_t k = 0; k < loop; k++) {
1589       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1590       AliTRDCalROC *calROC = object->GetCalROC(detector);
1591       Float_t min          = detobject->GetValue(detector);
1592       Int_t rowMax    = calROC->GetNrows();
1593       Int_t colMax    = calROC->GetNcols();
1594       for (Int_t row = 0; row < rowMax; row++) {
1595         for (Int_t col = 0; col < colMax; col++) {
1596           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1597           calROC->SetValue(col,row,value-min);
1598         } // Col
1599       } // Row
1600     } 
1601   }
1602   return object;    
1603
1604 }
1605 //_____________________________________________________________________________
1606 TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit)
1607 {
1608   //
1609   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1610   // This object has to be written in the database
1611   //
1612   
1613   // Create the DetObject
1614   AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1615
1616   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1617   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1618   Int_t detector = -1;
1619   Float_t value  = 0.0;
1620
1621   for (Int_t k = 0; k < loop; k++) {
1622     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1623     AliTRDCalROC *calROC = object->GetCalROC(detector);
1624     Int_t rowMax    = calROC->GetNrows();
1625     Int_t colMax    = calROC->GetNcols();
1626     for (Int_t row = 0; row < rowMax; row++) {
1627       for (Int_t col = 0; col < colMax; col++) {
1628         value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1629         calROC->SetValue(col,row,TMath::Abs(value));
1630       } // Col
1631     } // Row
1632   } 
1633
1634   return object;  
1635
1636 }
1637 //_____________________________________________________________________________
1638 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean)
1639 {
1640   //
1641   // It Creates the AliTRDCalDet object from AliTRDFitInfo
1642   // 0 successful fit 1 not successful fit
1643   // mean is the mean value over the successful fit
1644   // do not use it for t0: no meaning
1645   //
1646   
1647   // Create the CalObject
1648   AliTRDCalDet *object = new AliTRDCalDet(name,name);
1649   mean = 0.0;
1650   Int_t count = 0;
1651   
1652   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1653   if(loop != 540) {
1654     AliInfo("The Vector Fit is not complete! We initialise all outliers");
1655     for(Int_t k = 0; k < 540; k++){
1656       object->SetValue(k,1.0);
1657     }
1658   }
1659   Int_t detector = -1;
1660   Float_t value  = 0.0;
1661   
1662   for (Int_t k = 0; k < loop; k++) {
1663     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1664     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1665     if(value <= 0) object->SetValue(detector,1.0);
1666     else {
1667       object->SetValue(detector,0.0);
1668       mean += value;
1669       count++;
1670     }
1671   }
1672   if(count > 0) mean /= count;
1673   return object;  
1674 }
1675 //_____________________________________________________________________________
1676 TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean)
1677 {
1678   //
1679   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1680   // 0 not successful fit 1 successful fit
1681   // mean mean value over the successful fit
1682   //
1683   
1684   // Create the CalObject
1685   AliTRDCalPad *object = new AliTRDCalPad(name,name);
1686   mean = 0.0;
1687   Int_t count = 0;
1688   
1689   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1690   if(loop != 540) {
1691     AliInfo("The Vector Fit is not complete! We initialise all outliers");
1692     for(Int_t k = 0; k < 540; k++){
1693       AliTRDCalROC *calROC = object->GetCalROC(k);
1694       Int_t nchannels = calROC->GetNchannels();
1695       for(Int_t ch = 0; ch < nchannels; ch++){
1696         calROC->SetValue(ch,1.0);
1697       }
1698     }
1699   }
1700   Int_t detector = -1;
1701   Float_t value  = 0.0;
1702   
1703   for (Int_t k = 0; k < loop; k++) {
1704     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1705     AliTRDCalROC *calROC = object->GetCalROC(detector);
1706     Int_t nchannels    = calROC->GetNchannels();
1707     for (Int_t ch = 0; ch < nchannels; ch++) {
1708       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
1709       if(value <= 0) calROC->SetValue(ch,1.0);
1710       else {
1711         calROC->SetValue(ch,0.0);
1712         mean += value;
1713         count++;
1714       }
1715     } // channels
1716   }
1717   if(count > 0) mean /= count;
1718   return object;  
1719 }
1720 //_____________________________________________________________________________
1721 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
1722
1723   //
1724   // Set FitPH if 1 then each detector will be fitted
1725   //
1726
1727   if (periodeFitPH > 0) {
1728     fFitPHPeriode   = periodeFitPH; 
1729   }
1730   else {
1731     AliInfo("periodeFitPH must be higher than 0!");
1732   }
1733
1734 }
1735 //_____________________________________________________________________________
1736 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
1737
1738   //
1739   // The fit of the deposited charge distribution begins at
1740   // histo->Mean()/beginFitCharge
1741   // You can here set beginFitCharge
1742   //
1743
1744   if (beginFitCharge > 0) {
1745     fBeginFitCharge = beginFitCharge; 
1746   }
1747   else {
1748     AliInfo("beginFitCharge must be strict positif!");
1749   }
1750
1751 }
1752
1753 //_____________________________________________________________________________
1754 void AliTRDCalibraFit::SetT0Shift(Float_t t0Shift) 
1755
1756   //
1757   // The t0 calculated with the maximum positif slope is shift from t0Shift
1758   // You can here set t0Shift
1759   //
1760
1761   if (t0Shift > 0) {
1762     fT0Shift = t0Shift; 
1763   } 
1764   else {
1765     AliInfo("t0Shift must be strict positif!");
1766   }
1767
1768 }
1769
1770 //_____________________________________________________________________________
1771 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
1772
1773   //
1774   // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
1775   // You can here set rangeFitPRF
1776   //
1777
1778   if ((rangeFitPRF >    0) && 
1779       (rangeFitPRF <= 1.5)) {
1780     fRangeFitPRF = rangeFitPRF;
1781   } 
1782   else {
1783     AliInfo("rangeFitPRF must be between 0 and 1.0");
1784   }
1785
1786 }
1787
1788 //_____________________________________________________________________________
1789 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
1790
1791   //
1792   // Minimum entries for fitting
1793   //
1794
1795   if (minEntries >    0) {
1796     fMinEntries = minEntries;
1797   } 
1798   else {
1799     AliInfo("fMinEntries must be >= 0.");
1800   }
1801
1802 }
1803
1804 //_____________________________________________________________________________
1805 void AliTRDCalibraFit::SetRebin(Short_t rebin)
1806
1807   //
1808   // Rebin with rebin time less bins the Ch histo
1809   // You can set here rebin that should divide the number of bins of CH histo
1810   //
1811
1812   if (rebin > 0) {
1813     fRebin = rebin; 
1814     AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
1815   } 
1816   else {
1817     AliInfo("You have to choose a positiv value!");
1818   }
1819
1820 }
1821 //_____________________________________________________________________________
1822 Bool_t AliTRDCalibraFit::FillVectorFit()
1823 {
1824   //
1825   // For the Fit functions fill the vector Fit
1826   //
1827
1828   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1829
1830   Int_t ntotal = 1;
1831   if (GetChamber(fCountDet) == 2) {
1832     ntotal = 1728;
1833   }
1834   else {
1835     ntotal = 2304;
1836   }
1837
1838   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1839   Float_t *coef = new Float_t[ntotal];
1840   for (Int_t i = 0; i < ntotal; i++) {
1841     coef[i] = fCurrentCoefDetector[i];
1842   }
1843   
1844   Int_t detector = fCountDet;
1845   // Set
1846   fitInfo->SetCoef(coef);
1847   fitInfo->SetDetector(detector);
1848   fVectorFit.Add((TObject *) fitInfo);
1849
1850   return kTRUE;
1851
1852 }
1853 //_____________________________________________________________________________
1854 Bool_t AliTRDCalibraFit::FillVectorFit2()
1855 {
1856   //
1857   // For the Fit functions fill the vector Fit
1858   //
1859
1860   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
1861
1862   Int_t ntotal = 1;
1863   if (GetChamber(fCountDet) == 2) {
1864     ntotal = 1728;
1865   }
1866   else {
1867     ntotal = 2304;
1868   }
1869
1870   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
1871   Float_t *coef = new Float_t[ntotal];
1872   for (Int_t i = 0; i < ntotal; i++) {
1873     coef[i] = fCurrentCoefDetector2[i];
1874   }
1875   
1876   Int_t detector = fCountDet;
1877   // Set
1878   fitInfo->SetCoef(coef);
1879   fitInfo->SetDetector(detector);
1880   fVectorFit2.Add((TObject *) fitInfo);
1881
1882   return kTRUE;
1883
1884 }
1885 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1886 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
1887 {
1888   //
1889   // Init the number of expected bins and fDect1[i] fDect2[i] 
1890   //
1891
1892   gStyle->SetPalette(1);
1893   gStyle->SetOptStat(1111);
1894   gStyle->SetPadBorderMode(0);
1895   gStyle->SetCanvasColor(10);
1896   gStyle->SetPadLeftMargin(0.13);
1897   gStyle->SetPadRightMargin(0.01);
1898   
1899   // Mode groups of pads: the total number of bins!
1900   CalculNumberOfBinsExpected(i);
1901   
1902   // Quick verification that we have the good pad calibration mode!
1903   if (fNumberOfBinsExpected != nbins) {
1904     AliInfo("It doesn't correspond to the mode of pad group calibration!");
1905     return kFALSE;
1906   }
1907   
1908   // Security for fDebug 3 and 4
1909   if ((fDebugLevel >= 3) && 
1910       ((fDet[0] >  5) || 
1911        (fDet[1] >  4) || 
1912        (fDet[2] > 17))) {
1913     AliInfo("This detector doesn't exit!");
1914     return kFALSE;
1915   }
1916
1917   // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
1918   CalculDect1Dect2(i);
1919
1920  
1921   return kTRUE;
1922 }
1923 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1924 Bool_t AliTRDCalibraFit::InitFitCH()
1925 {
1926   //
1927   // Init the fVectorFitCH for normalisation
1928   // Init the histo for debugging 
1929   //
1930
1931   gDirectory = gROOT;
1932  
1933   fScaleFitFactor = 0.0;
1934   fCurrentCoefDetector   = new Float_t[2304];
1935   for (Int_t k = 0; k < 2304; k++) {
1936     fCurrentCoefDetector[k] = 0.0;    
1937   }
1938   fVectorFit.SetName("gainfactorscoefficients");
1939
1940   // fDebug == 0 nothing
1941   // fDebug == 1 and fFitVoir no histo
1942   if (fDebugLevel == 1) {
1943     if(!CheckFitVoir()) return kFALSE;
1944   }
1945   //Get the CalDet object
1946   if(fAccCDB){
1947     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
1948     if (!cal) {
1949       AliInfo("Could not get calibDB");
1950       return kFALSE;
1951     }
1952     if(fCalDet) delete fCalDet;
1953     fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
1954   }
1955   else{
1956     Float_t devalue = 1.0;
1957     if(fCalDet) delete fCalDet;
1958     fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1959     for(Int_t k = 0; k < 540; k++){
1960       fCalDet->SetValue(k,devalue);
1961     }
1962   }
1963   return kTRUE;
1964   
1965 }
1966 //____________Functions for initialising the AliTRDCalibraFit in the code_________
1967 Bool_t AliTRDCalibraFit::InitFitPH()
1968 {
1969   //
1970   // Init the arrays of results 
1971   // Init the histos for debugging 
1972   //
1973
1974   gDirectory = gROOT;
1975   fVectorFit.SetName("driftvelocitycoefficients");
1976   fVectorFit2.SetName("t0coefficients");
1977
1978   fCurrentCoefDetector   = new Float_t[2304];
1979   for (Int_t k = 0; k < 2304; k++) {
1980     fCurrentCoefDetector[k] = 0.0;    
1981   }
1982
1983   fCurrentCoefDetector2   = new Float_t[2304];
1984   for (Int_t k = 0; k < 2304; k++) {
1985     fCurrentCoefDetector2[k] = 0.0;    
1986   }
1987  
1988   //fDebug == 0 nothing
1989   // fDebug == 1 and fFitVoir no histo
1990   if (fDebugLevel == 1) {
1991     if(!CheckFitVoir()) return kFALSE;
1992   }
1993   //Get the CalDet object
1994   if(fAccCDB){
1995     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
1996     if (!cal) {
1997       AliInfo("Could not get calibDB");
1998       return kFALSE;
1999     }
2000     if(fCalDet) delete fCalDet;
2001     if(fCalDet2) delete fCalDet2;
2002     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2003     fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det())); 
2004   }
2005   else{
2006     Float_t devalue  = 1.5;
2007     Float_t devalue2 = 0.0; 
2008     if(fCalDet) delete fCalDet;
2009     if(fCalDet2) delete fCalDet2;
2010     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2011     fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2012     for(Int_t k = 0; k < 540; k++){
2013       fCalDet->SetValue(k,devalue);
2014       fCalDet2->SetValue(k,devalue2);
2015     }
2016   }
2017   return kTRUE;
2018 }
2019 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2020 Bool_t AliTRDCalibraFit::InitFitPRF()
2021 {
2022   //
2023   // Init the calibration mode (Nz, Nrphi), the histograms for
2024   // debugging the fit methods if fDebug > 0, 
2025   //
2026   
2027   gDirectory = gROOT;
2028   fVectorFit.SetName("prfwidthcoefficients");
2029  
2030   fCurrentCoefDetector   = new Float_t[2304];
2031   for (Int_t k = 0; k < 2304; k++) {
2032     fCurrentCoefDetector[k] = 0.0;    
2033   }
2034   
2035   // fDebug == 0 nothing
2036   // fDebug == 1 and fFitVoir no histo
2037   if (fDebugLevel == 1) {
2038     if(!CheckFitVoir()) return kFALSE;
2039   }
2040   return kTRUE;
2041 }
2042 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2043 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2044 {
2045   //
2046   // Init the fCalDet, fVectorFit fCurrentCoefDetector 
2047   //
2048   
2049   gDirectory = gROOT;
2050  
2051   fCurrentCoefDetector   = new Float_t[2304];
2052   fCurrentCoefDetector2  = new Float_t[2304];
2053   for (Int_t k = 0; k < 2304; k++) {
2054     fCurrentCoefDetector[k]  = 0.0;
2055     fCurrentCoefDetector2[k] = 0.0;    
2056   }
2057
2058   //printf("test0\n");
2059   
2060   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2061   if (!cal) {
2062     AliInfo("Could not get calibDB");
2063     return kFALSE;
2064   }
2065   
2066   //Get the CalDet object
2067   if(fAccCDB){
2068     if(fCalDet) delete fCalDet;
2069     if(fCalDet2) delete fCalDet2;
2070     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2071     //printf("test1\n");
2072     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2073     //printf("test2\n");
2074     for(Int_t k = 0; k < 540; k++){
2075       fCalDet2->SetValue(k,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField));
2076     }
2077     //printf("test3\n");
2078   }
2079   else{
2080     Float_t devalue  = 1.5;
2081     Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField); 
2082     if(fCalDet) delete fCalDet;
2083     if(fCalDet2) delete fCalDet2;
2084     //printf("test1\n");
2085     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2086     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2087     //printf("test2\n");
2088     for(Int_t k = 0; k < 540; k++){
2089       fCalDet->SetValue(k,devalue);
2090       fCalDet2->SetValue(k,devalue2);
2091     }
2092     //printf("test3\n");
2093   }
2094   return kTRUE;
2095 }
2096
2097 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2098 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2099 {
2100   //
2101   // Init the current detector where we are fCountDet and the
2102   // next fCount for the functions Fit... 
2103   //
2104
2105   // Loop on the Xbins of ch!!
2106   fCountDet = -1; // Current detector
2107   fCount    =  0; // To find the next detector
2108   
2109   // If fDebug >= 3
2110   if (fDebugLevel >= 3) {
2111     // Set countdet to the detector
2112     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2113     // Set counter to write at the end of the detector
2114     fCount = fDect2;
2115     // Get the right calib objects
2116     SetCalROC(i);
2117   }
2118   if(fDebugLevel == 1) {
2119     fCountDet = 0;
2120     fCalibraMode->CalculXBins(fCountDet,i);
2121     while(fCalibraMode->GetXbins(i) <=fFitVoir){
2122       fCountDet++;
2123       fCalibraMode->CalculXBins(fCountDet,i);
2124     }      
2125     fCount    = fCalibraMode->GetXbins(i);
2126     fCountDet--;
2127     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2128     fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2129     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2130                                        ,(Int_t) GetChamber(fCountDet)
2131                                        ,(Int_t) GetSector(fCountDet),i);
2132   }
2133 }
2134 //_______________________________________________________________________________
2135 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2136 {
2137   //
2138   // Calculate the number of bins expected (calibration groups)
2139   //
2140   
2141   fNumberOfBinsExpected = 0;
2142   fCalibraMode->ModePadCalibration(2,i);
2143   fCalibraMode->ModePadFragmentation(0,2,0,i);
2144   fCalibraMode->SetDetChamb2(i);
2145   if (fDebugLevel > 1) {
2146     AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2147   }
2148   fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2149   fCalibraMode->ModePadCalibration(0,i);
2150   fCalibraMode->ModePadFragmentation(0,0,0,i);
2151   fCalibraMode->SetDetChamb0(i);
2152   if (fDebugLevel > 1) {
2153     AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2154   }
2155   fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2156  
2157 }
2158 //_______________________________________________________________________________
2159 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2160 {
2161   //
2162   // Calculate the range of fits
2163   //
2164   
2165   fDect1 = -1;
2166   fDect2 = -1;
2167   if (fDebugLevel == 1) {
2168     fDect1 = fFitVoir;
2169     fDect2 = fDect1 +1;
2170   }
2171   if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2172     fDect1 = 0;
2173     fDect2 = fNumberOfBinsExpected;
2174   }
2175   if (fDebugLevel >= 3) {
2176     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2177     fCalibraMode->CalculXBins(fCountDet,i);
2178     fDect1 = fCalibraMode->GetXbins(i);
2179     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2180     fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2181     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2182                                        ,(Int_t) GetChamber(fCountDet)
2183                                        ,(Int_t) GetSector(fCountDet),i);
2184     // Set for the next detector
2185     fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2186   }
2187 }
2188 //_______________________________________________________________________________
2189 Bool_t AliTRDCalibraFit::CheckFitVoir()
2190 {
2191   //
2192   // Check if fFitVoir is in the range
2193   //
2194   
2195   if (fFitVoir < fNumberOfBinsExpected) {
2196     AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2197   }
2198   else {
2199     AliInfo("fFitVoir is out of range of the histo!");
2200     return kFALSE;
2201   }
2202   return kTRUE;
2203 }
2204 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2205 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2206 {
2207   //
2208   // See if we are in a new detector and update the
2209   // variables fNfragZ and fNfragRphi if yes 
2210   // Will never happen for only one detector (3 and 4)
2211   // Doesn't matter for 2
2212   //
2213   if (fCount == idect) {
2214      // On en est au detector
2215      fCountDet += 1;
2216      // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2217      fCalibraMode->ModePadCalibration((Int_t) GetChamber(fCountDet),i);
2218      fCalibraMode->ModePadFragmentation((Int_t) GetPlane(fCountDet)
2219                                 ,(Int_t) GetChamber(fCountDet)
2220                           ,(Int_t) GetSector(fCountDet),i);
2221      // Set for the next detector
2222      fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2223      // calib objects
2224      SetCalROC(i);
2225     }
2226 }
2227 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2228 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2229 {
2230   //
2231   // Reconstruct the min pad row, max pad row, min pad col and
2232   // max pad col of the calibration group for the Fit functions
2233   //
2234   if (fDebugLevel !=  1) {
2235     fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2236   }
2237 }
2238 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2239 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2240 {
2241   //
2242   // For the case where there are not enough entries in the histograms
2243   // of the calibration group, the value present in the choosen database
2244   // will be put. A negativ sign enables to know that a fit was not possible.
2245   //
2246   
2247   if (fDebugLevel == 1) {
2248     AliInfo("The element has not enough statistic to be fitted");
2249   }
2250   
2251   else {
2252
2253     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2254                  ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2255     
2256     // Calcul the coef from the database choosen
2257     CalculChargeCoefMean(kFALSE);
2258
2259     //chamber 2, not chamber 2
2260     Int_t factor = 0;
2261     if(GetChamber(fCountDet) == 2) factor = 12;
2262     else factor = 16;
2263     
2264     // Fill the fCurrentCoefDetector with negative value to say: not fitted
2265     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2266       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2267         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2268       }
2269     }
2270     
2271     //Put default value negative
2272     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2273     fCurrentCoefE   = 0.0;
2274    
2275     FillFillCH(idect);
2276   }
2277   
2278   return kTRUE;
2279 }
2280
2281
2282 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2283 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect)
2284 {
2285   //
2286   // For the case where there are not enough entries in the histograms
2287   // of the calibration group, the value present in the choosen database
2288   // will be put. A negativ sign enables to know that a fit was not possible.
2289   //
2290   if (fDebugLevel == 1) {
2291     AliInfo("The element has not enough statistic to be fitted");
2292   }
2293   else {
2294
2295     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2296                  ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2297
2298     CalculVdriftCoefMean();
2299     CalculT0CoefMean();
2300   
2301     //chamber 2 and not chamber 2
2302     Int_t factor = 0;
2303     if(GetChamber(fCountDet) == 2) factor = 12;
2304     else factor = 16;
2305
2306
2307     // Fill the fCurrentCoefDetector 2
2308     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2309       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2310         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2311         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1];
2312       }
2313     }
2314
2315     // Put the default value
2316     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2317     fCurrentCoefE    = 0.0;
2318     fCurrentCoef2[0] = fCurrentCoef2[1];
2319     fCurrentCoefE2   = 0.0;
2320      
2321     FillFillPH(idect);
2322     
2323   }
2324   
2325   return kTRUE;
2326
2327 }
2328
2329
2330 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2331 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2332 {
2333   //
2334   // For the case where there are not enough entries in the histograms
2335   // of the calibration group, the value present in the choosen database
2336   // will be put. A negativ sign enables to know that a fit was not possible.
2337   //
2338   
2339   if (fDebugLevel == 1) {
2340     AliInfo("The element has not enough statistic to be fitted");
2341   }
2342   else {
2343     
2344     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2345                  ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2346     
2347     CalculPRFCoefMean();
2348     
2349     // chamber 2 and not chamber 2
2350     Int_t factor = 0;
2351     if(GetChamber(fCountDet) == 2) factor = 12;
2352     else factor = 16;
2353
2354     
2355     // Fill the fCurrentCoefDetector
2356     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2357       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2358         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1];
2359       }
2360     }
2361
2362     // Put the default value
2363     fCurrentCoef[0] = -fCurrentCoef[1];
2364     fCurrentCoefE   = 0.0;
2365     
2366     FillFillPRF(idect);
2367   }
2368   
2369   return kTRUE;
2370
2371 }
2372 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2373 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
2374 {
2375   //
2376   // For the case where there are not enough entries in the histograms
2377   // of the calibration group, the value present in the choosen database
2378   // will be put. A negativ sign enables to know that a fit was not possible.
2379   //
2380   
2381   // Calcul the coef from the database choosen
2382   CalculVdriftLorentzCoef();
2383
2384   Int_t factor = 0;
2385   if(GetChamber(fCountDet) == 2) factor = 1728;
2386   else factor = 2304;
2387     
2388     
2389   // Fill the fCurrentCoefDetector
2390   for (Int_t k = 0; k < factor; k++) {
2391     fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
2392     // should be negative
2393     fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
2394   }
2395    
2396   
2397   //Put default opposite sign
2398   fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2399   fCurrentCoefE    = 0.0;
2400   fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
2401   fCurrentCoefE2 = 0.0; 
2402   
2403   FillFillLinearFitter();
2404     
2405   return kTRUE;
2406 }
2407
2408 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2409 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
2410 {
2411   //
2412   // Fill the coefficients found with the fits or other
2413   // methods from the Fit functions
2414   //
2415
2416   if (fDebugLevel != 1) {
2417     
2418     Int_t factor = 0;
2419     if(GetChamber(fCountDet) == 2) factor = 12;
2420     else factor = 16; 
2421     
2422     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2423       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2424         fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2425       }
2426     }
2427     
2428     FillFillCH(idect);
2429     
2430   }
2431
2432   return kTRUE;
2433
2434 }
2435 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2436 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect)
2437 {
2438   //
2439   // Fill the coefficients found with the fits or other
2440   // methods from the Fit functions
2441   //
2442
2443   if (fDebugLevel != 1) {
2444
2445     Int_t factor = 0;
2446     if(GetChamber(fCountDet) == 2) factor = 12;
2447     else factor = 16; 
2448     
2449     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2450       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2451         fCurrentCoefDetector[(Int_t)(j*factor+k)]  = fCurrentCoef[0];
2452         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
2453       }
2454     }                
2455     FillFillPH(idect);
2456   }
2457   return kTRUE;
2458 }
2459 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2460 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
2461 {
2462   //
2463   // Fill the coefficients found with the fits or other
2464   // methods from the Fit functions
2465   //
2466   
2467   if (fDebugLevel != 1) {
2468
2469     Int_t factor = 0;
2470     if(GetChamber(fCountDet) == 2) factor = 12;
2471     else factor = 16; 
2472     
2473     // Pointer to the branch
2474     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2475       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2476         fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
2477       }
2478     }
2479     FillFillPRF(idect);   
2480   }
2481
2482   return kTRUE;
2483
2484 }
2485 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2486 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
2487 {
2488   //
2489   // Fill the coefficients found with the fits or other
2490   // methods from the Fit functions
2491   //
2492   
2493   Int_t factor = 0;
2494   if(GetChamber(fCountDet) == 2) factor = 1728;
2495   else factor = 2304; 
2496   
2497   // Pointer to the branch
2498   for (Int_t k = 0; k < factor; k++) {
2499     fCurrentCoefDetector[k]  = fCurrentCoef[0];
2500     fCurrentCoefDetector2[k] = fCurrentCoef2[0];
2501   }
2502   
2503   FillFillLinearFitter();
2504   
2505   return kTRUE;
2506
2507 }
2508 //________________________________________________________________________________
2509 void AliTRDCalibraFit::FillFillCH(Int_t idect)
2510 {
2511   //
2512   // DebugStream and fVectorFit
2513   //
2514
2515   // End of one detector
2516   if ((idect == (fCount-1))) {
2517     FillVectorFit();
2518     // Reset
2519     for (Int_t k = 0; k < 2304; k++) {
2520       fCurrentCoefDetector[k] = 0.0;
2521     }
2522   }
2523
2524   if(fDebugLevel > 1){ 
2525
2526     if ( !fDebugStreamer ) {
2527       //debug stream
2528       TDirectory *backup = gDirectory;
2529       fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2530       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2531     } 
2532     
2533     Int_t   detector   = fCountDet;
2534     Int_t   caligroup  = idect;
2535     Short_t rowmin     = fCalibraMode->GetRowMin(0);
2536     Short_t rowmax     = fCalibraMode->GetRowMax(0);
2537     Short_t colmin     = fCalibraMode->GetColMin(0);
2538     Short_t colmax     = fCalibraMode->GetColMax(0);
2539     Float_t gf         = fCurrentCoef[0]; 
2540     Float_t gfs        = fCurrentCoef[1]; 
2541     Float_t gfE        = fCurrentCoefE;
2542     
2543     (*fDebugStreamer) << "GAIN" <<
2544       "detector=" << detector <<
2545       "caligroup=" << caligroup <<
2546       "rowmin=" << rowmin <<
2547       "rowmax=" << rowmax <<
2548       "colmin=" << colmin <<
2549       "colmax=" << colmax <<
2550       "gf=" << gf <<
2551       "gfs=" << gfs <<
2552       "gfE=" << gfE <<
2553       "\n"; 
2554     
2555   }
2556 }
2557 //________________________________________________________________________________
2558 void AliTRDCalibraFit::FillFillPH(Int_t idect)
2559 {
2560   //
2561   // DebugStream and fVectorFit and fVectorFit2
2562   //
2563   
2564   // End of one detector
2565     if ((idect == (fCount-1))) {
2566       FillVectorFit();
2567       FillVectorFit2();
2568       // Reset
2569       for (Int_t k = 0; k < 2304; k++) {
2570         fCurrentCoefDetector[k] = 0.0;
2571         fCurrentCoefDetector2[k] = 0.0;
2572       }
2573     }
2574
2575     if(fDebugLevel > 1){ 
2576
2577       if ( !fDebugStreamer ) {
2578         //debug stream
2579         TDirectory *backup = gDirectory;
2580         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2581         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2582       } 
2583       
2584        
2585       Int_t   detector     = fCountDet;
2586       Int_t   caligroup    = idect;
2587       Short_t rowmin       = fCalibraMode->GetRowMin(1);
2588       Short_t rowmax       = fCalibraMode->GetRowMax(1);
2589       Short_t colmin       = fCalibraMode->GetColMin(1);
2590       Short_t colmax       = fCalibraMode->GetColMax(1);
2591       Float_t vf           = fCurrentCoef[0]; 
2592       Float_t vs           = fCurrentCoef[1]; 
2593       Float_t vfE          = fCurrentCoefE;
2594       Float_t t0f          = fCurrentCoef2[0]; 
2595       Float_t t0s          = fCurrentCoef2[1]; 
2596       Float_t t0E          = fCurrentCoefE2;
2597    
2598
2599
2600       (* fDebugStreamer) << "PH"<<
2601         "detector="<<detector<<
2602         "caligroup="<<caligroup<<
2603         "rowmin="<<rowmin<<
2604         "rowmax="<<rowmax<<
2605         "colmin="<<colmin<<
2606         "colmax="<<colmax<<
2607         "vf="<<vf<<
2608         "vs="<<vs<<
2609         "vfE="<<vfE<<
2610         "t0f="<<t0f<<
2611         "t0s="<<t0s<<
2612         "t0E="<<t0E<<
2613         "\n";  
2614     }
2615
2616 }
2617 //________________________________________________________________________________
2618 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
2619 {
2620   //
2621   // DebugStream and fVectorFit
2622   //
2623
2624     // End of one detector
2625     if ((idect == (fCount-1))) {
2626       FillVectorFit();
2627       // Reset
2628       for (Int_t k = 0; k < 2304; k++) {
2629         fCurrentCoefDetector[k] = 0.0;
2630       }
2631     }
2632
2633     
2634     if(fDebugLevel > 1){
2635
2636       if ( !fDebugStreamer ) {
2637         //debug stream
2638         TDirectory *backup = gDirectory;
2639         fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2640         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2641       } 
2642       
2643       Int_t   detector     = fCountDet;
2644       Int_t   plane        = GetPlane(fCountDet);
2645       Int_t   caligroup    = idect;
2646       Short_t rowmin       = fCalibraMode->GetRowMin(2);
2647       Short_t rowmax       = fCalibraMode->GetRowMax(2);
2648       Short_t colmin       = fCalibraMode->GetColMin(2);
2649       Short_t colmax       = fCalibraMode->GetColMax(2);
2650       Float_t widf         = fCurrentCoef[0]; 
2651       Float_t wids         = fCurrentCoef[1]; 
2652       Float_t widfE        = fCurrentCoefE;
2653
2654       (* fDebugStreamer) << "PRF"<<
2655         "detector="<<detector<<
2656         "plane="<<plane<<
2657         "caligroup="<<caligroup<<
2658         "rowmin="<<rowmin<<
2659         "rowmax="<<rowmax<<
2660         "colmin="<<colmin<<
2661         "colmax="<<colmax<<
2662         "widf="<<widf<<
2663         "wids="<<wids<<
2664         "widfE="<<widfE<<
2665         "\n";  
2666     }
2667
2668 }
2669 //________________________________________________________________________________
2670 void AliTRDCalibraFit::FillFillLinearFitter()
2671 {
2672   //
2673   // DebugStream and fVectorFit
2674   //
2675
2676   // End of one detector
2677   FillVectorFit();
2678   FillVectorFit2();
2679   
2680   
2681   // Reset
2682   for (Int_t k = 0; k < 2304; k++) {
2683   fCurrentCoefDetector[k]  = 0.0;
2684   fCurrentCoefDetector2[k] = 0.0;
2685   }
2686   
2687
2688   if(fDebugLevel > 1){
2689
2690     if ( !fDebugStreamer ) {
2691       //debug stream
2692       TDirectory *backup = gDirectory;
2693       fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
2694       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2695     } 
2696     
2697     //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
2698     AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetPlane(fCountDet),GetChamber(fCountDet));
2699     Float_t rowmd            = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
2700     Float_t r                = AliTRDgeometry::GetTime0(GetPlane(fCountDet)); 
2701     Float_t tiltangle        = padplane->GetTiltingAngle();
2702     Int_t   detector         = fCountDet;
2703     Int_t   chamber          = GetChamber(fCountDet);
2704     Int_t   plane            = GetChamber(fCountDet);
2705     Float_t vf               = fCurrentCoef[0]; 
2706     Float_t vs               = fCurrentCoef[1]; 
2707     Float_t vfE              = fCurrentCoefE;
2708     Float_t lorentzangler    = fCurrentCoef2[0];
2709     Float_t Elorentzangler   = fCurrentCoefE2;
2710     Float_t lorentzangles    = fCurrentCoef2[1];
2711    
2712     (* fDebugStreamer) << "LinearFitter"<<
2713       "detector="<<detector<<
2714       "chamber="<<chamber<<
2715       "plane="<<plane<<
2716       "rowmd="<<rowmd<<
2717       "r="<<r<<
2718       "tiltangle="<<tiltangle<<
2719       "vf="<<vf<<
2720       "vs="<<vs<<
2721       "vfE="<<vfE<<
2722       "lorentzangler="<<lorentzangler<<
2723       "Elorentzangler="<<Elorentzangler<<
2724       "lorentzangles="<<lorentzangles<<
2725       "\n";  
2726   }
2727   
2728 }
2729 //
2730 //____________Calcul Coef Mean_________________________________________________
2731 //
2732 //_____________________________________________________________________________
2733 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
2734 {
2735   //
2736   // For the detector Dect calcul the mean time 0
2737   // for the calibration group idect from the choosen database
2738   //
2739
2740   fCurrentCoef2[1] = 0.0;
2741   if(fDebugLevel != 1){
2742     if((fCalibraMode->GetNz(1) > 0) ||
2743        (fCalibraMode->GetNrphi(1) > 0)) {
2744       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2745         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2746           fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2747         }
2748       }
2749       fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2750     }
2751     else {
2752       if(!fAccCDB){
2753         fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
2754       }
2755       else{
2756         for(Int_t row = 0; row < fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)); row++){
2757           for(Int_t col = 0; col < fGeo->GetColMax(GetPlane(fCountDet)); col++){
2758             fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
2759           }
2760         }
2761         fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetPlane(fCountDet),GetChamber(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetPlane(fCountDet))));
2762       }
2763     }
2764   }
2765   return kTRUE;
2766 }
2767
2768 //_____________________________________________________________________________
2769 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
2770 {
2771   //
2772   // For the detector Dect calcul the mean gain factor
2773   // for the calibration group idect from the choosen database
2774   //
2775
2776   fCurrentCoef[1] = 0.0;
2777   if(fDebugLevel != 1){
2778     if ((fCalibraMode->GetNz(0)    > 0) || 
2779         (fCalibraMode->GetNrphi(0) > 0)) {
2780       for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
2781         for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
2782           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2783           if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2784         }
2785       }
2786       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
2787     }
2788     else {
2789       //Per detectors
2790       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2791       if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
2792     }    
2793   }
2794   return kTRUE;
2795 }
2796 //_____________________________________________________________________________
2797 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
2798 {
2799   //
2800   // For the detector Dect calcul the mean sigma of pad response
2801   // function for the calibration group idect from the choosen database
2802   //
2803   
2804   fCurrentCoef[1] = 0.0;
2805   if(fDebugLevel != 1){
2806     for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
2807       for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
2808         fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
2809       }
2810     }
2811     fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
2812   }
2813   return kTRUE;
2814 }
2815 //_____________________________________________________________________________
2816 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
2817 {
2818   //
2819   // For the detector dect calcul the mean drift velocity for the
2820   // calibration group idect from the choosen database
2821   //
2822
2823   fCurrentCoef[1] = 0.0;
2824   if(fDebugLevel != 1){
2825     if ((fCalibraMode->GetNz(1)    > 0) || 
2826         (fCalibraMode->GetNrphi(1) > 0)) {
2827       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
2828         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
2829           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
2830         }
2831       }
2832       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
2833     }
2834     else {
2835       //per detectors
2836       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
2837     }  
2838   }
2839   return kTRUE;
2840 }
2841 //_____________________________________________________________________________
2842 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
2843 {
2844   //
2845   // For the detector fCountDet, mean drift velocity and tan lorentzangle
2846   //
2847
2848   fCurrentCoef[1]  = fCalDet->GetValue(fCountDet);
2849   fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet); 
2850
2851   return kTRUE;
2852 }
2853 //_____________________________________________________________________________
2854 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t plane) const
2855 {
2856   //
2857   // Default width of the PRF if there is no database as reference
2858   //
2859   switch(plane)
2860     {
2861       // default database
2862       //case 0:  return 0.515;
2863       //case 1:  return 0.502;
2864       //case 2:  return 0.491;
2865       //case 3:  return 0.481;
2866       //case 4:  return 0.471;
2867       //case 5:  return 0.463;
2868       //default: return 0.0;
2869
2870       // fit database
2871     case 0:  return 0.538429;
2872     case 1:  return 0.524302;
2873     case 2:  return 0.511591;
2874     case 3:  return 0.500140;
2875     case 4:  return 0.489821;
2876     case 5:  return 0.480524;
2877     default: return 0.0;
2878   }
2879 }
2880 //________________________________________________________________________________
2881 void AliTRDCalibraFit::SetCalROC(Int_t i)
2882 {
2883   //
2884   // Set the calib object for fCountDet
2885   //
2886
2887   Float_t value = 0.0;
2888   
2889   //Get the CalDet object
2890   if(fAccCDB){
2891     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2892     if (!cal) {
2893       AliInfo("Could not get calibDB");
2894       return;
2895     }
2896     switch (i)
2897       {
2898       case 0: 
2899         if(fCalROC) delete fCalROC;
2900         fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet))); 
2901         break;
2902       case 1:
2903         if(fCalROC)  delete fCalROC;
2904         if(fCalROC2) delete fCalROC2;
2905         fCalROC  = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet))); 
2906         fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet))); 
2907         break;
2908       case 2:
2909         if(fCalROC) delete fCalROC; 
2910         fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break; 
2911       default: return;
2912       }
2913   }
2914   else{
2915     switch (i)
2916       {
2917       case 0:
2918         if(fCalROC) delete fCalROC;
2919         fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet)); 
2920         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2921           fCalROC->SetValue(k,1.0);
2922         }
2923         break;
2924       case 1:
2925         if(fCalROC)  delete fCalROC;
2926         if(fCalROC2) delete fCalROC2;
2927         fCalROC  = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2928         fCalROC2 = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet));
2929         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2930           fCalROC->SetValue(k,1.0);
2931           fCalROC2->SetValue(k,0.0);
2932         }
2933         break;
2934       case 2:
2935         if(fCalROC) delete fCalROC;
2936         value = GetPRFDefault(GetPlane(fCountDet));
2937         fCalROC = new AliTRDCalROC(GetPlane(fCountDet),GetChamber(fCountDet)); 
2938         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
2939           fCalROC->SetValue(k,value);
2940         }
2941         break;
2942       default: return; 
2943       }
2944   }
2945   
2946 }
2947 //____________Fit Methods______________________________________________________
2948
2949 //_____________________________________________________________________________
2950 void AliTRDCalibraFit::FitPente(TH1* projPH)
2951 {
2952   //
2953   // Slope methode for the drift velocity
2954   //
2955   
2956   // Constants
2957   const Float_t kDrWidth = AliTRDgeometry::DrThick();
2958   Int_t binmax           = 0;
2959   Int_t binmin           = 0;
2960   fPhd[0]                = 0.0;
2961   fPhd[1]                = 0.0;
2962   fPhd[2]                = 0.0;
2963   Int_t ju               = 0;
2964   fCurrentCoefE          = 0.0;
2965   fCurrentCoefE2         = 0.0;
2966   fCurrentCoef[0]        = 0.0;
2967   fCurrentCoef2[0]       = 0.0;
2968   TLine *line            = new TLine();
2969
2970   // Some variables
2971   TAxis   *xpph    = projPH->GetXaxis();
2972   Int_t    nbins   = xpph->GetNbins();
2973   Double_t lowedge = xpph->GetBinLowEdge(1);
2974   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
2975   Double_t widbins = (upedge - lowedge) / nbins;
2976   Double_t limit   = upedge + 0.5 * widbins; 
2977   Bool_t put = kTRUE;
2978
2979   // Beginning of the signal
2980   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
2981   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
2982     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
2983   }
2984   binmax = (Int_t) pentea->GetMaximumBin();
2985   if (binmax <= 1) {
2986     binmax = 2;
2987     AliInfo("Put the binmax from 1 to 2 to enable the fit");
2988   }
2989   if (binmax >= nbins) {
2990     binmax = nbins-1;
2991     put = kFALSE;
2992     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
2993   }
2994   pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
2995   Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
2996   Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
2997   Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
2998   Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
2999   if (l3P2am != 0) {
3000     fPhd[0] = -(l3P1am / (2 * l3P2am));
3001   }
3002   if(!fTakeTheMaxPH){
3003     if((l3P1am != 0.0) && (l3P2am != 0.0)){
3004       fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3005     }
3006   }
3007   // Amplification region
3008   binmax = 0;
3009   ju     = 0;
3010   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3011     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3012       binmax = kbin;
3013       ju     = 1;
3014     }
3015   }
3016   if (binmax <= 1) {
3017     binmax = 2;
3018     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3019   }
3020   if (binmax >= nbins) {
3021     binmax = nbins-1;
3022     put = kFALSE;
3023     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3024   }
3025   projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3026   Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3027   Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3028   Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3029   Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3030   if (l3P2amf != 0) {
3031     fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3032   }
3033   if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3034     fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3035   }
3036   if(fTakeTheMaxPH){
3037     fCurrentCoefE2 = fCurrentCoefE;
3038   }
3039   // Drift region
3040   TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3041   for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
3042     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3043   }
3044   binmin = 0;
3045   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3046   if (binmin <= 1) {
3047     binmin = 2;
3048     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3049   }
3050   if (binmin >= nbins) {
3051     binmin = nbins-1;
3052     put = kFALSE;
3053     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3054   }
3055   pente->Fit("pol2"
3056             ,"0MR"
3057             ,""
3058             ,TMath::Max(pente->GetBinCenter(binmin-1),             0.0)
3059             ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
3060   Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
3061   Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
3062   Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
3063   Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
3064   if (l3P2dr != 0) {
3065     fPhd[2] = -(l3P1dr / (2 * l3P2dr));
3066   }
3067   if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
3068     fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2]; 
3069   }
3070   Float_t fPhdt0 = 0.0;
3071   if(fTakeTheMaxPH) fPhdt0 = fPhd[1];
3072   else fPhdt0 = fPhd[0];
3073
3074   if ((fPhd[2] > fPhd[0]) && 
3075       (fPhd[2] > fPhd[1]) && 
3076       (fPhd[1] > fPhd[0]) &&
3077       (put)) {
3078     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
3079     fNumberFitSuccess++;
3080
3081     if (fPhdt0 >= 0.0) {
3082       fCurrentCoef2[0] = (fPhdt0 - fT0Shift) / widbins;
3083       if (fCurrentCoef2[0] < -1.0) {
3084         fCurrentCoef2[0] = fCurrentCoef2[1];
3085       }
3086     }
3087     else {
3088       fCurrentCoef2[0] = fCurrentCoef2[1];
3089     }
3090
3091   }
3092   else {
3093     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3094     fCurrentCoef2[0] = fCurrentCoef2[1];
3095   }
3096
3097   if (fDebugLevel == 1) {
3098     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
3099     cpentei->cd();
3100     projPH->Draw();
3101     line->SetLineColor(2);
3102     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
3103     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
3104     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
3105     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
3106     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
3107     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
3108     AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
3109     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
3110     cpentei2->cd();
3111     pentea->Draw();
3112     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
3113     cpentei3->cd();
3114     pente->Draw();
3115   }
3116   else {
3117     delete pentea;
3118     delete pente;
3119   }
3120 }
3121 //_____________________________________________________________________________
3122 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
3123 {
3124   //
3125   // Slope methode but with polynomes de Lagrange
3126   //
3127   
3128   // Constants
3129   const Float_t kDrWidth = AliTRDgeometry::DrThick();
3130   Int_t binmax      = 0;
3131   Int_t binmin      = 0;
3132   Double_t    *x    = new Double_t[5];
3133   Double_t    *y    = new Double_t[5];
3134   x[0]              = 0.0;
3135   x[1]              = 0.0;
3136   x[2]              = 0.0;
3137   x[3]              = 0.0;
3138   x[4]              = 0.0;
3139   y[0]              = 0.0;
3140   y[1]              = 0.0;
3141   y[2]              = 0.0;
3142   y[3]              = 0.0;
3143   y[4]              = 0.0;
3144   fPhd[0]           = 0.0;
3145   fPhd[1]           = 0.0;
3146   fPhd[2]           = 0.0;
3147   Int_t ju          = 0;
3148   fCurrentCoefE     = 0.0;
3149   fCurrentCoefE2    = 1.0;
3150   fCurrentCoef[0]   = 0.0;
3151   fCurrentCoef2[0]  = 0.0;
3152   TLine *line = new TLine();
3153   TF1 * polynome = 0x0;
3154   TF1 * polynomea = 0x0;
3155   TF1 * polynomeb = 0x0;
3156   Double_t *c = 0x0;
3157   
3158   // Some variables
3159   TAxis   *xpph    = projPH->GetXaxis();
3160   Int_t    nbins   = xpph->GetNbins();
3161   Double_t lowedge = xpph->GetBinLowEdge(1);
3162   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
3163   Double_t widbins = (upedge - lowedge) / nbins;
3164   Double_t limit   = upedge + 0.5 * widbins;
3165
3166   
3167   Bool_t put = kTRUE;
3168
3169   // Beginning of the signal
3170   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3171   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
3172     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3173   }
3174
3175   binmax = (Int_t) pentea->GetMaximumBin();
3176
3177   Double_t minnn = 0.0;
3178   Double_t maxxx = 0.0;
3179
3180   Int_t kase = nbins-binmax;
3181   
3182   switch(kase)
3183     {
3184     case 0:
3185       put = kFALSE;
3186       break;
3187     case 1:
3188       minnn = pentea->GetBinCenter(binmax-2);
3189       maxxx = pentea->GetBinCenter(binmax);
3190       x[0] = pentea->GetBinCenter(binmax-2);
3191       x[1] = pentea->GetBinCenter(binmax-1);
3192       x[2] = pentea->GetBinCenter(binmax);
3193       y[0] = pentea->GetBinContent(binmax-2);
3194       y[1] = pentea->GetBinContent(binmax-1);
3195       y[2] = pentea->GetBinContent(binmax);
3196       c = CalculPolynomeLagrange2(x,y);
3197       AliInfo("At the limit for beginning!");
3198       break;  
3199     case 2:
3200       minnn = pentea->GetBinCenter(binmax-2);
3201       maxxx = pentea->GetBinCenter(binmax+1);
3202       x[0] = pentea->GetBinCenter(binmax-2);
3203       x[1] = pentea->GetBinCenter(binmax-1);
3204       x[2] = pentea->GetBinCenter(binmax);
3205       x[3] = pentea->GetBinCenter(binmax+1);
3206       y[0] = pentea->GetBinContent(binmax-2);
3207       y[1] = pentea->GetBinContent(binmax-1);
3208       y[2] = pentea->GetBinContent(binmax);
3209       y[3] = pentea->GetBinContent(binmax+1);
3210       c = CalculPolynomeLagrange3(x,y);
3211       break;
3212     default:
3213       switch(binmax){
3214       case 0:
3215         put = kFALSE;
3216         break;
3217       case 1:
3218         minnn = pentea->GetBinCenter(binmax);
3219         maxxx = pentea->GetBinCenter(binmax+2);
3220         x[0] = pentea->GetBinCenter(binmax);
3221         x[1] = pentea->GetBinCenter(binmax+1);
3222         x[2] = pentea->GetBinCenter(binmax+2);
3223         y[0] = pentea->GetBinContent(binmax);
3224         y[1] = pentea->GetBinContent(binmax+1);
3225         y[2] = pentea->GetBinContent(binmax+2);
3226         c = CalculPolynomeLagrange2(x,y);
3227         break;
3228       case 2:
3229         minnn = pentea->GetBinCenter(binmax-1);
3230         maxxx = pentea->GetBinCenter(binmax+2);
3231         x[0] = pentea->GetBinCenter(binmax-1);
3232         x[1] = pentea->GetBinCenter(binmax);
3233         x[2] = pentea->GetBinCenter(binmax+1);
3234         x[3] = pentea->GetBinCenter(binmax+2);
3235         y[0] = pentea->GetBinContent(binmax-1);
3236         y[1] = pentea->GetBinContent(binmax);
3237         y[2] = pentea->GetBinContent(binmax+1);
3238         y[3] = pentea->GetBinContent(binmax+2);
3239         c = CalculPolynomeLagrange3(x,y);
3240         break;
3241       default:
3242         minnn = pentea->GetBinCenter(binmax-2);
3243         maxxx = pentea->GetBinCenter(binmax+2);
3244         x[0] = pentea->GetBinCenter(binmax-2);
3245         x[1] = pentea->GetBinCenter(binmax-1);
3246         x[2] = pentea->GetBinCenter(binmax);
3247         x[3] = pentea->GetBinCenter(binmax+1);
3248         x[4] = pentea->GetBinCenter(binmax+2);
3249         y[0] = pentea->GetBinContent(binmax-2);
3250         y[1] = pentea->GetBinContent(binmax-1);
3251         y[2] = pentea->GetBinContent(binmax);
3252         y[3] = pentea->GetBinContent(binmax+1);
3253         y[4] = pentea->GetBinContent(binmax+2);
3254         c = CalculPolynomeLagrange4(x,y);
3255         break;
3256       }
3257       break;
3258     }
3259   
3260   
3261   if(put) {
3262     polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
3263     polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3264       
3265     Double_t step = (maxxx-minnn)/10000;
3266     Double_t l = minnn;
3267     Double_t maxvalue = 0.0;
3268     Double_t placemaximum = minnn;
3269     for(Int_t o = 0; o < 10000; o++){
3270       if(o == 0) maxvalue = polynomeb->Eval(l);
3271       if(maxvalue < (polynomeb->Eval(l))){
3272         maxvalue = polynomeb->Eval(l);
3273         placemaximum = l;
3274       }
3275       l += step;
3276     }
3277     fPhd[0] = placemaximum;
3278   }
3279   
3280   // Amplification region
3281   binmax = 0;
3282   ju     = 0;
3283   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3284     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3285       binmax = kbin;
3286       ju     = 1;
3287     }
3288   }
3289    
3290   Double_t minn = 0.0;
3291   Double_t maxx = 0.0;
3292   x[0] = 0.0;
3293   x[1] = 0.0;
3294   x[2] = 0.0;
3295   x[3] = 0.0;
3296   x[4] = 0.0;
3297   y[0] = 0.0;
3298   y[1] = 0.0;
3299   y[2] = 0.0;
3300   y[3] = 0.0;
3301   y[4] = 0.0;
3302
3303   Int_t    kase1 = nbins - binmax;
3304
3305   //Determination of minn and maxx
3306   //case binmax = nbins
3307   //pol2
3308   switch(kase1)
3309     {
3310     case 0:
3311       minn = projPH->GetBinCenter(binmax-2);
3312       maxx = projPH->GetBinCenter(binmax);
3313       x[0] = projPH->GetBinCenter(binmax-2);
3314       x[1] = projPH->GetBinCenter(binmax-1);
3315       x[2] = projPH->GetBinCenter(binmax);
3316       y[0] = projPH->GetBinContent(binmax-2);
3317       y[1] = projPH->GetBinContent(binmax-1);
3318       y[2] = projPH->GetBinContent(binmax);
3319       c = CalculPolynomeLagrange2(x,y);
3320       //AliInfo("At the limit for the drift!");
3321       break;
3322     case 1:
3323       minn = projPH->GetBinCenter(binmax-2);
3324       maxx = projPH->GetBinCenter(binmax+1);
3325       x[0] = projPH->GetBinCenter(binmax-2);
3326       x[1] = projPH->GetBinCenter(binmax-1);
3327       x[2] = projPH->GetBinCenter(binmax);
3328       x[3] = projPH->GetBinCenter(binmax+1);
3329       y[0] = projPH->GetBinContent(binmax-2);
3330       y[1] = projPH->GetBinContent(binmax-1);
3331       y[2] = projPH->GetBinContent(binmax);
3332       y[3] = projPH->GetBinContent(binmax+1);
3333       c = CalculPolynomeLagrange3(x,y);
3334       break;
3335     default:
3336       switch(binmax)
3337         {
3338         case 0:
3339           put = kFALSE;
3340           break;
3341         case 1:
3342           minn = projPH->GetBinCenter(binmax);
3343           maxx = projPH->GetBinCenter(binmax+2);
3344           x[0] = projPH->GetBinCenter(binmax);
3345           x[1] = projPH->GetBinCenter(binmax+1);
3346           x[2] = projPH->GetBinCenter(binmax+2);
3347           y[0] = projPH->GetBinContent(binmax);
3348           y[1] = projPH->GetBinContent(binmax+1);
3349           y[2] = projPH->GetBinContent(binmax+2);
3350           c = CalculPolynomeLagrange2(x,y);
3351           break;
3352         case 2:
3353           minn = projPH->GetBinCenter(binmax-1);
3354           maxx = projPH->GetBinCenter(binmax+2);
3355           x[0] = projPH->GetBinCenter(binmax-1);
3356           x[1] = projPH->GetBinCenter(binmax);
3357           x[2] = projPH->GetBinCenter(binmax+1);
3358           x[3] = projPH->GetBinCenter(binmax+2);
3359           y[0] = projPH->GetBinContent(binmax-1);
3360           y[1] = projPH->GetBinContent(binmax);
3361           y[2] = projPH->GetBinContent(binmax+1);
3362           y[3] = projPH->GetBinContent(binmax+2);
3363           c = CalculPolynomeLagrange3(x,y);
3364           break;
3365         default:
3366           minn = projPH->GetBinCenter(binmax-2);
3367           maxx = projPH->GetBinCenter(binmax+2);
3368           x[0] = projPH->GetBinCenter(binmax-2);
3369           x[1] = projPH->GetBinCenter(binmax-1);
3370           x[2] = projPH->GetBinCenter(binmax);
3371           x[3] = projPH->GetBinCenter(binmax+1);
3372           x[4] = projPH->GetBinCenter(binmax+2);
3373           y[0] = projPH->GetBinContent(binmax-2);
3374           y[1] = projPH->GetBinContent(binmax-1);
3375           y[2] = projPH->GetBinContent(binmax);
3376           y[3] = projPH->GetBinContent(binmax+1);
3377           y[4] = projPH->GetBinContent(binmax+2);
3378           c = CalculPolynomeLagrange4(x,y);
3379           break;
3380         }
3381       break;
3382     }
3383   
3384   if(put) {
3385     polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
3386     polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3387        
3388     Double_t step = (maxx-minn)/1000;
3389     Double_t l = minn;
3390     Double_t maxvalue = 0.0;
3391     Double_t placemaximum = minn;
3392     for(Int_t o = 0; o < 1000; o++){
3393       if(o == 0) maxvalue = polynomea->Eval(l);
3394       if(maxvalue < (polynomea->Eval(l))){
3395         maxvalue = polynomea->Eval(l);
3396         placemaximum = l;
3397       }
3398       l += step;
3399     }
3400     fPhd[1] = placemaximum;
3401   }
3402   
3403   // Drift region
3404   TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
3405   for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
3406     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3407   }
3408   binmin = 0;
3409   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3410
3411   //should not happen
3412   if (binmin <= 1) {
3413     binmin = 2;
3414     put = 1;
3415     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3416   }
3417   
3418   //check
3419   if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) put = kFALSE;
3420   if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE;
3421   
3422   x[0] = 0.0;
3423   x[1] = 0.0;
3424   x[2] = 0.0;
3425   x[3] = 0.0;
3426   x[4] = 0.0;
3427   y[0] = 0.0;
3428   y[1] = 0.0;
3429   y[2] = 0.0;
3430   y[3] = 0.0;
3431   y[4] = 0.0;
3432   Double_t min = 0.0;
3433   Double_t max = 0.0;
3434   Bool_t case1 = kFALSE;
3435   Bool_t case2 = kFALSE;
3436   Bool_t case4 = kFALSE;
3437
3438   //Determination of min and max
3439   //case binmin <= nbins-3
3440   //pol4 case 3
3441   if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3442     min = pente->GetBinCenter(binmin-2);
3443     max = pente->GetBinCenter(binmin+2);
3444     x[0] = pente->GetBinCenter(binmin-2);
3445     x[1] = pente->GetBinCenter(binmin-1);
3446     x[2] = pente->GetBinCenter(binmin);
3447     x[3] = pente->GetBinCenter(binmin+1);
3448     x[4] = pente->GetBinCenter(binmin+2);
3449     y[0] = pente->GetBinContent(binmin-2);
3450     y[1] = pente->GetBinContent(binmin-1);
3451     y[2] = pente->GetBinContent(binmin);
3452     y[3] = pente->GetBinContent(binmin+1);
3453     y[4] = pente->GetBinContent(binmin+2);
3454     //Calcul the polynome de Lagrange
3455     c = CalculPolynomeLagrange4(x,y);
3456     //richtung +/-
3457     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3458        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3459     if(((binmin+3) <= (nbins-1)) &&
3460        (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
3461        ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
3462        (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) put = kFALSE;
3463     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
3464        (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE;
3465     if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
3466        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE;
3467   }
3468   //case binmin = nbins-2
3469   //pol3 case 1
3470   if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3471      (case1)){
3472     min = pente->GetBinCenter(binmin-2);
3473     max = pente->GetBinCenter(binmin+1);
3474     x[0] = pente->GetBinCenter(binmin-2);
3475     x[1] = pente->GetBinCenter(binmin-1);
3476     x[2] = pente->GetBinCenter(binmin);
3477     x[3] = pente->GetBinCenter(binmin+1);
3478     y[0] = pente->GetBinContent(binmin-2);
3479     y[1] = pente->GetBinContent(binmin-1);
3480     y[2] = pente->GetBinContent(binmin);
3481     y[3] = pente->GetBinContent(binmin+1);
3482     //Calcul the polynome de Lagrange
3483     c = CalculPolynomeLagrange3(x,y);
3484     //richtung +: nothing
3485     //richtung -
3486     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case2 = kTRUE;
3487   }
3488   //pol3 case 4
3489   if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3490      (case4)){
3491     min = pente->GetBinCenter(binmin-1);
3492     max = pente->GetBinCenter(binmin+2);
3493     x[0] = pente->GetBinCenter(binmin-1);
3494     x[1] = pente->GetBinCenter(binmin);
3495     x[2] = pente->GetBinCenter(binmin+1);
3496     x[3] = pente->GetBinCenter(binmin+2);
3497     y[0] = pente->GetBinContent(binmin-1);
3498     y[1] = pente->GetBinContent(binmin);
3499     y[2] = pente->GetBinContent(binmin+1);
3500     y[3] = pente->GetBinContent(binmin+2);
3501     //Calcul the polynome de Lagrange
3502     c = CalculPolynomeLagrange3(x,y);
3503     //richtung +
3504     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) case2 = kTRUE;
3505   }
3506   //pol2 case 5
3507   if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
3508     min = pente->GetBinCenter(binmin);
3509     max = pente->GetBinCenter(binmin+2);
3510     x[0] = pente->GetBinCenter(binmin);
3511     x[1] = pente->GetBinCenter(binmin+1);
3512     x[2] = pente->GetBinCenter(binmin+2);
3513     y[0] = pente->GetBinContent(binmin);
3514     y[1] = pente->GetBinContent(binmin+1);
3515     y[2] = pente->GetBinContent(binmin+2);
3516     //Calcul the polynome de Lagrange
3517     c = CalculPolynomeLagrange2(x,y);
3518     //richtung +
3519     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) put = kFALSE;
3520   }
3521   //pol2 case 2
3522   if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
3523      (case2)){
3524     min = pente->GetBinCenter(binmin-1);
3525     max = pente->GetBinCenter(binmin+1);
3526     x[0] = pente->GetBinCenter(binmin-1);
3527     x[1] = pente->GetBinCenter(binmin);
3528     x[2] = pente->GetBinCenter(binmin+1);
3529     y[0] = pente->GetBinContent(binmin-1);
3530     y[1] = pente->GetBinContent(binmin);
3531     y[2] = pente->GetBinContent(binmin+1);
3532     //Calcul the polynome de Lagrange
3533     c = CalculPolynomeLagrange2(x,y);
3534     //richtung +: nothing
3535     //richtung -: nothing
3536   }
3537   //case binmin = nbins-1
3538   //pol2 case 0
3539   if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
3540     min = pente->GetBinCenter(binmin-2);
3541     max = pente->GetBinCenter(binmin);
3542     x[0] = pente->GetBinCenter(binmin-2);
3543     x[1] = pente->GetBinCenter(binmin-1);
3544     x[2] = pente->GetBinCenter(binmin);
3545     y[0] = pente->GetBinContent(binmin-2);
3546     y[1] = pente->GetBinContent(binmin-1);
3547     y[2] = pente->GetBinContent(binmin);
3548     //Calcul the polynome de Lagrange
3549     c = CalculPolynomeLagrange2(x,y);
3550     //AliInfo("At the limit for the drift!");
3551     //fluctuation too big!
3552     //richtung +: nothing
3553     //richtung -
3554     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) put = kFALSE;
3555   }
3556   if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
3557     put = kFALSE;
3558     AliInfo("At the limit for the drift and not usable!");
3559   }
3560
3561   //pass
3562   if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3563     put = kFALSE;
3564     AliInfo("For the drift...problem!");
3565   }
3566   //pass but should not happen
3567   if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){
3568     put = kFALSE;
3569     AliInfo("For the drift...problem!");
3570   }
3571   
3572   if(put) {
3573     polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
3574     polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
3575     //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
3576     Double_t step = (max-min)/1000;
3577     Double_t l = min;
3578     Double_t minvalue = 0.0;
3579     Double_t placeminimum = min;
3580     for(Int_t o = 0; o < 1000; o++){
3581       if(o == 0) minvalue = polynome->Eval(l);
3582       if(minvalue > (polynome->Eval(l))){
3583         minvalue = polynome->Eval(l);
3584         placeminimum = l;
3585       }
3586       l += step;
3587     }
3588     fPhd[2] = placeminimum;
3589   }