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: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
  32 // _fDebug = 2: a comparaison of the coefficients found and the default values
  33 //              in the choosen database.
  34 //              fCoef , histogram of the coefs as function of the calibration group number
  35 //              fDelta , histogram of the relative difference of the coef with the default
  36 //                        value in the database as function of the calibration group number
  37 //              fError , dirstribution of this relative difference
  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 <TLine.h>
  50 #include <TH1I.h>
  51 #include <TStyle.h>
  52 #include <TProfile2D.h>
  53 #include <TCanvas.h>
  54 #include <TGraphErrors.h>
  55 #include <TObjArray.h>
  56 #include <TH1.h>
  57 #include <TH1F.h>
  58 #include <TF1.h>
  59 #include <TAxis.h>
  60 #include <TMath.h>
  61 #include <TDirectory.h>
  62 #include <TROOT.h>
  63 #include <TTreeStream.h>
  64 #include <TLinearFitter.h>
  65 #include <TVectorD.h>
  66 #include <TArrayF.h>
  67
  68 #include "AliLog.h"
  69 #include "AliMathBase.h"
  70
  71 #include "AliTRDCalibraFit.h"
  72 #include "AliTRDCalibraMode.h"
  73 #include "AliTRDCalibraVector.h"
  74 #include "AliTRDCalibraVdriftLinearFit.h"
  75 #include "AliTRDcalibDB.h"
  76 #include "AliTRDgeometry.h"
  77 #include "AliTRDpadPlane.h"
  78 #include "AliTRDgeometry.h"
  79 #include "AliTRDCommonParam.h"
  80 #include "./Cal/AliTRDCalROC.h"
  81 #include "./Cal/AliTRDCalPad.h"
  82 #include "./Cal/AliTRDCalDet.h"
  83
  84
  85 ClassImp(AliTRDCalibraFit)
  86
  87 AliTRDCalibraFit* AliTRDCalibraFit::fgInstance = 0;
  88 Bool_t AliTRDCalibraFit::fgTerminated = kFALSE;
  89
  90 //_____________singleton implementation_________________________________________________
  91 AliTRDCalibraFit *AliTRDCalibraFit::Instance()
  92 {
  93   //
  94   // Singleton implementation
  95   //
  96
  97   if (fgTerminated != kFALSE) {
  98     return 0;
  99   }
 100
 101   if (fgInstance == 0) {
 102     fgInstance = new AliTRDCalibraFit();
 103   }
 104
 105   return fgInstance;
 106
 107 }
 108 //______________________________________________________________________________________
 109 void AliTRDCalibraFit::Terminate()
 110 {
 111   //
 112   // Singleton implementation
 113   // Deletes the instance of this class
 114   //
 115
 116   fgTerminated = kTRUE;
 117
 118   if (fgInstance != 0) {
 119     delete fgInstance;
 120     fgInstance = 0;
 121   }
 122
 123 }
 124 //______________________________________________________________________________________
 125 AliTRDCalibraFit::AliTRDCalibraFit()
 126   :TObject()
 127   ,fGeo(0)
 128   ,fNumberOfBinsExpected(0)
 129   ,fMethod(0)
 130   ,fBeginFitCharge(3.5)
 131   ,fFitPHPeriode(1)
 132   ,fTakeTheMaxPH(kTRUE)
 133   ,fT0Shift0(0.124797)
 134   ,fT0Shift1(0.267451)
 135   ,fRangeFitPRF(1.0)
 136   ,fAccCDB(kFALSE)
 137   ,fMinEntries(800)
 138   ,fRebin(1)
 139   ,fNumberFit(0)
 140   ,fNumberFitSuccess(0)
 141   ,fNumberEnt(0)
 142   ,fStatisticMean(0.0)
 143   ,fDebugStreamer(0x0)
 144   ,fDebugLevel(0)
 145   ,fFitVoir(0)
 146   ,fMagneticField(0.5)
 147   ,fCalibraMode(new AliTRDCalibraMode())
 148   ,fCurrentCoefE(0.0)
 149   ,fCurrentCoefE2(0.0)
 150   ,fDect1(0)
 151   ,fDect2(0)
 152   ,fScaleFitFactor(0.0)
 153   ,fEntriesCurrent(0)
 154   ,fCountDet(0)
 155   ,fCount(0)
 156   ,fNbDet(0)
 157   ,fCalDet(0x0)
 158   ,fCalROC(0x0)
 159   ,fCalDet2(0x0)
 160   ,fCalROC2(0x0)
 161   ,fCurrentCoefDetector(0x0)
 162   ,fCurrentCoefDetector2(0x0)
 163   ,fVectorFit(0)
 164   ,fVectorFit2(0)
 165 {
 166   //
 167   // Default constructor
 168   //
 169
 170   fGeo         = new AliTRDgeometry();
 171
 172   // Current variables initialised
 173   for (Int_t k = 0; k < 2; k++) {
 174     fCurrentCoef[k]  = 0.0;
 175     fCurrentCoef2[k] = 0.0;
 176   }
 177   for (Int_t i = 0; i < 3; i++) {
 178     fPhd[i]          = 0.0;
 179     fDet[i]          = 0;
 180   }
 181
 182 }
 183 //______________________________________________________________________________________
 184 AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c)
 185 :TObject(c)
 186 ,fGeo(0)
 187 ,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
 188 ,fMethod(c.fMethod)
 189 ,fBeginFitCharge(c.fBeginFitCharge)
 190 ,fFitPHPeriode(c.fFitPHPeriode)
 191 ,fTakeTheMaxPH(c.fTakeTheMaxPH)
 192 ,fT0Shift0(c.fT0Shift0)
 193 ,fT0Shift1(c.fT0Shift1)
 194 ,fRangeFitPRF(c.fRangeFitPRF)
 195 ,fAccCDB(c.fAccCDB)
 196 ,fMinEntries(c.fMinEntries)
 197 ,fRebin(c.fRebin)
 198 ,fNumberFit(c.fNumberFit)
 199 ,fNumberFitSuccess(c.fNumberFitSuccess)
 200 ,fNumberEnt(c.fNumberEnt)
 201 ,fStatisticMean(c.fStatisticMean)
 202 ,fDebugStreamer(0x0)
 203 ,fDebugLevel(c.fDebugLevel)
 204 ,fFitVoir(c.fFitVoir)
 205 ,fMagneticField(c.fMagneticField)
 206 ,fCalibraMode(0x0)
 207 ,fCurrentCoefE(c.fCurrentCoefE)
 208 ,fCurrentCoefE2(c.fCurrentCoefE2)
 209 ,fDect1(c.fDect1)
 210 ,fDect2(c.fDect2)
 211 ,fScaleFitFactor(c.fScaleFitFactor)
 212 ,fEntriesCurrent(c.fEntriesCurrent)
 213 ,fCountDet(c.fCountDet)
 214 ,fCount(c.fCount)
 215 ,fNbDet(c.fNbDet)
 216 ,fCalDet(0x0)
 217 ,fCalROC(0x0)
 218 ,fCalDet2(0x0)
 219 ,fCalROC2(0x0)
 220 ,fCurrentCoefDetector(0x0)
 221 ,fCurrentCoefDetector2(0x0)
 222 ,fVectorFit(0)
 223 ,fVectorFit2(0)
 224 {
 225   //
 226   // Copy constructor
 227   //
 228
 229   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
 230
 231   //Current variables initialised
 232   for (Int_t k = 0; k < 2; k++) {
 233     fCurrentCoef[k]  = 0.0;
 234     fCurrentCoef2[k] = 0.0;
 235   }
 236   for (Int_t i = 0; i < 3; i++) {
 237     fPhd[i]          = 0.0;
 238     fDet[i]          = 0;
 239   }
 240   if(c.fCalDet) fCalDet   = new AliTRDCalDet(*c.fCalDet);
 241   if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2);
 242
 243   if(c.fCalROC) fCalROC   = new AliTRDCalROC(*c.fCalROC);
 244   if(c.fCalROC2) fCalROC  = new AliTRDCalROC(*c.fCalROC2);
 245
 246   fVectorFit.SetName(c.fVectorFit.GetName());
 247   for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
 248     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
 249     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetDetector();
 250     Int_t ntotal = 1;
 251     if (GetStack(detector) == 2) {
 252       ntotal = 1728;
 253     }
 254     else {
 255       ntotal = 2304;
 256     }
 257     Float_t *coef = new Float_t[ntotal];
 258     for (Int_t i = 0; i < ntotal; i++) {
 259       coef[i] = ((AliTRDFitInfo *)c.fVectorFit.UncheckedAt(k))->GetCoef()[i];
 260     }
 261     fitInfo->SetCoef(coef);
 262     fitInfo->SetDetector(detector);
 263     fVectorFit.Add((TObject *) fitInfo);
 264   }
 265   fVectorFit.SetName(c.fVectorFit.GetName());
 266   for(Int_t k = 0; k < c.fVectorFit2.GetEntriesFast(); k++){
 267     AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
 268     Int_t detector         = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetDetector();
 269     Int_t ntotal = 1;
 270     if (GetStack(detector) == 2) {
 271       ntotal = 1728;
 272     }
 273     else {
 274       ntotal = 2304;
 275     }
 276     Float_t *coef = new Float_t[ntotal];
 277     for (Int_t i = 0; i < ntotal; i++) {
 278       coef[i] = ((AliTRDFitInfo *)c.fVectorFit2.UncheckedAt(k))->GetCoef()[i];
 279     }
 280     fitInfo->SetCoef(coef);
 281     fitInfo->SetDetector(detector);
 282     fVectorFit2.Add((TObject *) fitInfo);
 283   }
 284   if (fGeo) {
 285     delete fGeo;
 286   }
 287   fGeo = new AliTRDgeometry();
 288
 289 }
 290 //____________________________________________________________________________________
 291 AliTRDCalibraFit::~AliTRDCalibraFit()
 292 {
 293   //
 294   // AliTRDCalibraFit destructor
 295   //
 296   if ( fDebugStreamer ) delete fDebugStreamer;
 297   if ( fCalDet )  delete fCalDet;
 298   if ( fCalDet2 ) delete fCalDet2;
 299   if ( fCalROC )  delete fCalROC;
 300   if ( fCalROC2 ) delete fCalROC2;
 301   if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
 302   if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
 303   fVectorFit.Delete();
 304   fVectorFit2.Delete();
 305   if (fGeo) {
 306     delete fGeo;
 307   }
 308
 309 }
 310 //_____________________________________________________________________________
 311 void AliTRDCalibraFit::Destroy()
 312 {
 313   //
 314   // Delete instance
 315   //
 316
 317   if (fgInstance) {
 318     delete fgInstance;
 319     fgInstance = 0x0;
 320   }
 321
 322 }
 323 //_____________________________________________________________________________
 324 void AliTRDCalibraFit::DestroyDebugStreamer()
 325 {
 326   //
 327   // Delete DebugStreamer
 328   //
 329
 330   if ( fDebugStreamer ) delete fDebugStreamer;
 331   fDebugStreamer = 0x0;
 332
 333 }
 334 //__________________________________________________________________________________
 335 void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
 336 {
 337   //
 338   // From the drift velocity and t0
 339   // return the position of the peak and maximum negative slope
 340   //
 341
 342   const Float_t kDrWidth = AliTRDgeometry::DrThick();    // drift region
 343   Double_t widbins = 0.1;                                // 0.1 mus
 344
 345   //peak and maxnegslope in mus
 346   Double_t begind = t0*widbins + fT0Shift0;
 347   Double_t peakd  = t0*widbins + fT0Shift1;
 348   Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
 349
 350   // peak and maxnegslope in timebin
 351   begin = TMath::Nint(begind*widbins);
 352   peak  = TMath::Nint(peakd*widbins);
 353   end   = TMath::Nint(maxnegslope*widbins);
 354
 355 }
 356 //____________Functions fit Online CH2d________________________________________
 357 Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
 358 {
 359   //
 360   // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
 361   // calibration group normalized the resulted coefficients (to 1 normally)
 362   //
 363
 364   // Set the calibration mode
 365   const char *name = ch->GetTitle();
 366   if(!SetModeCalibration(name,0)) return kFALSE;
 367
 368   // Number of Ybins (detectors or groups of pads)
 369   Int_t    nbins   = ch->GetNbinsX();// charge
 370   Int_t    nybins  = ch->GetNbinsY();// groups number
 371   if (!InitFit(nybins,0)) {
 372     return kFALSE;
 373   }
 374   if (!InitFitCH()) {
 375     return kFALSE;
 376   }
 377   fStatisticMean        = 0.0;
 378   fNumberFit            = 0;
 379   fNumberFitSuccess     = 0;
 380   fNumberEnt            = 0;
 381   // Init fCountDet and fCount
 382   InitfCountDetAndfCount(0);
 383   // Beginning of the loop betwwen dect1 and dect2
 384   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 385     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...
 386     UpdatefCountDetAndfCount(idect,0);
 387     ReconstructFitRowMinRowMax(idect, 0);
 388     // Take the histo
 389     TH1I *projch = (TH1I *) ch->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
 390     projch->SetDirectory(0);
 391     // Number of entries for this calibration group
 392     Double_t nentries = 0.0;
 393     Double_t mean = 0.0;
 394     for (Int_t k = 0; k < nbins; k++) {
 395       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 396       nentries += ch->GetBinContent(binnb);
 397       mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
 398       projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
 399     }
 400     projch->SetEntries(nentries);
 401     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 402     if (nentries > 0) {
 403       fNumberEnt++;
 404       mean /= nentries;
 405     }
 406     // Rebin and statistic stuff
 407     if (fRebin > 1) {
 408       projch = ReBin((TH1I *) projch);
 409     }
 410     // This detector has not enough statistics or was off
 411     if (nentries <= fMinEntries) {
 412       NotEnoughStatisticCH(idect);
 413       if (fDebugLevel != 1) {
 414         delete projch;
 415       }
 416       continue;
 417     }
 418     // Statistics of the group fitted
 419     fStatisticMean += nentries;
 420     fNumberFit++;
 421     //Method choosen
 422     switch(fMethod)
 423       {
 424       case 0: FitMeanW((TH1 *) projch, nentries); break;
 425       case 1: FitMean((TH1 *) projch, nentries, mean); break;
 426       case 2: FitCH((TH1 *) projch, mean); break;
 427       case 3: FitBisCH((TH1 *) projch, mean); break;
 428       default: return kFALSE;
 429       }
 430     // Fill Infos Fit
 431     FillInfosFitCH(idect);
 432     // Memory!!!
 433     if (fDebugLevel != 1) {
 434       delete projch;
 435     }
 436   } // Boucle object
 437   // Normierungcharge
 438   if (fDebugLevel != 1) {
 439     NormierungCharge();
 440   }
 441   // Mean Statistic
 442   if (fNumberFit > 0) {
 443     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));
 444     fStatisticMean = fStatisticMean / fNumberFit;
 445   }
 446   else {
 447     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 448   }
 449   delete fDebugStreamer;
 450   fDebugStreamer = 0x0;
 451
 452   return kTRUE;
 453 }
 454 //____________Functions fit Online CH2d________________________________________
 455 Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect)
 456 {
 457   //
 458   // Reconstruct a 1D histo from the vectorCH for each calibration group,
 459   // fit the histo, normalized the resulted coefficients (to 1 normally)
 460   //
 461
 462   // Set the calibraMode
 463   const char *name = calvect->GetNameCH();
 464   if(!SetModeCalibration(name,0)) return kFALSE;
 465
 466   // Number of Xbins (detectors or groups of pads)
 467   if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) {
 468     return kFALSE;
 469   }
 470   if (!InitFitCH()) {
 471     return kFALSE;
 472   }
 473   fStatisticMean        = 0.0;
 474   fNumberFit            = 0;
 475   fNumberFitSuccess     = 0;
 476   fNumberEnt            = 0;
 477   // Init fCountDet and fCount
 478   InitfCountDetAndfCount(0);
 479   // Beginning of the loop between dect1 and dect2
 480   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 481     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
 482     UpdatefCountDetAndfCount(idect,0);
 483     ReconstructFitRowMinRowMax(idect,0);
 484     // Take the histo
 485     Double_t nentries = 0.0;
 486     Double_t mean = 0.0;
 487     if(!calvect->GetCHEntries(fCountDet)) {
 488       NotEnoughStatisticCH(idect);
 489       continue;
 490     }
 491
 492     TString tname("CH");
 493     tname += idect;
 494     TH1F *projch  = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname);
 495     projch->SetDirectory(0);
 496     for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) {
 497       nentries += projch->GetBinContent(k+1);
 498       mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
 499     }
 500     if (nentries > 0) {
 501       fNumberEnt++;
 502       mean /= nentries;
 503     }
 504     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 505     // Rebin
 506     if (fRebin >  1) {
 507       projch = ReBin((TH1F *) projch);
 508     }
 509     // This detector has not enough statistics or was not found in VectorCH
 510     if (nentries <= fMinEntries) {
 511       NotEnoughStatisticCH(idect);
 512       continue;
 513     }
 514     // Statistic of the histos fitted
 515     fStatisticMean += nentries;
 516     fNumberFit++;
 517     //Method choosen
 518     switch(fMethod)
 519       {
 520       case 0: FitMeanW((TH1 *) projch, nentries); break;
 521       case 1: FitMean((TH1 *) projch, nentries, mean); break;
 522       case 2: FitCH((TH1 *) projch, mean); break;
 523       case 3: FitBisCH((TH1 *) projch, mean); break;
 524       default: return kFALSE;
 525       }
 526     // Fill Infos Fit
 527     FillInfosFitCH(idect);
 528   } // Boucle object
 529   // Normierungcharge
 530   if (fDebugLevel != 1) {
 531     NormierungCharge();
 532   }
 533   // Mean Statistics
 534   if (fNumberFit > 0) {
 535     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));
 536     fStatisticMean = fStatisticMean / fNumberFit;
 537   }
 538   else {
 539     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 540   }
 541   delete fDebugStreamer;
 542   fDebugStreamer = 0x0;
 543   return kTRUE;
 544 }
 545 //________________functions fit Online PH2d____________________________________
 546 Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
 547 {
 548   //
 549   // Take the 1D profiles (average pulse height), projections of the 2D PH
 550   // on the Xaxis, for each calibration group
 551   // Reconstruct a drift velocity
 552   // A first calibration of T0 is also made  using the same method
 553   //
 554
 555   // Set the calibration mode
 556   const char *name = ph->GetTitle();
 557   if(!SetModeCalibration(name,1)) return kFALSE;
 558
 559   // Number of Xbins (detectors or groups of pads)
 560   Int_t    nbins   = ph->GetNbinsX();// time
 561   Int_t    nybins  = ph->GetNbinsY();// calibration group
 562   if (!InitFit(nybins,1)) {
 563     return kFALSE;
 564   }
 565   if (!InitFitPH()) {
 566     return kFALSE;
 567   }
 568   fStatisticMean        = 0.0;
 569   fNumberFit            = 0;
 570   fNumberFitSuccess     = 0;
 571   fNumberEnt            = 0;
 572   // Init fCountDet and fCount
 573   InitfCountDetAndfCount(1);
 574   // Beginning of the loop
 575   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 576     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
 577     UpdatefCountDetAndfCount(idect,1);
 578     ReconstructFitRowMinRowMax(idect,1);
 579     // Take the histo
 580     TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
 581     projph->SetDirectory(0);
 582     // Number of entries for this calibration group
 583     Double_t nentries = 0;
 584     for (Int_t k = 0; k < nbins; k++) {
 585       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 586       nentries += ph->GetBinEntries(binnb);
 587     }
 588     if (nentries > 0) {
 589       fNumberEnt++;
 590     }
 591     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 592     // This detector has not enough statistics or was off
 593     if (nentries  <= fMinEntries) {
 594       //printf("Not enough statistic!\n");
 595       NotEnoughStatisticPH(idect,nentries);
 596       if (fDebugLevel != 1) {
 597         delete projph;
 598       }
 599       continue;
 600     }
 601     // Statistics of the histos fitted
 602     fNumberFit++;
 603     fStatisticMean += nentries;
 604     // Calcul of "real" coef
 605     CalculVdriftCoefMean();
 606     CalculT0CoefMean();
 607     //Method choosen
 608     switch(fMethod)
 609       {
 610       case 0: FitLagrangePoly((TH1 *) projph); break;
 611       case 1: FitPente((TH1 *) projph); break;
 612       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
 613       default: return kFALSE;
 614       }
 615     // Fill the tree if end of a detector or only the pointer to the branch!!!
 616     FillInfosFitPH(idect,nentries);
 617     // Memory!!!
 618     if (fDebugLevel != 1) {
 619       delete projph;
 620     }
 621   } // Boucle object
 622   // Mean Statistic
 623   if (fNumberFit > 0) {
 624     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));
 625     fStatisticMean = fStatisticMean / fNumberFit;
 626   }
 627   else {
 628     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 629   }
 630   delete fDebugStreamer;
 631   fDebugStreamer = 0x0;
 632   return kTRUE;
 633 }
 634 //____________Functions fit Online PH2d________________________________________
 635 Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
 636 {
 637   //
 638   // Reconstruct the average pulse height from the vectorPH for each
 639   // calibration group
 640   // Reconstruct a drift velocity
 641   // A first calibration of T0 is also made  using the same method (slope method)
 642   //
 643
 644   // Set the calibration mode
 645   const char *name = calvect->GetNamePH();
 646   if(!SetModeCalibration(name,1)) return kFALSE;
 647
 648   // Number of Xbins (detectors or groups of pads)
 649   if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
 650     return kFALSE;
 651   }
 652   if (!InitFitPH()) {
 653     return kFALSE;
 654   }
 655   fStatisticMean        = 0.0;
 656   fNumberFit            = 0;
 657   fNumberFitSuccess     = 0;
 658   fNumberEnt            = 0;
 659   // Init fCountDet and fCount
 660   InitfCountDetAndfCount(1);
 661   // Beginning of the loop
 662   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 663     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
 664     UpdatefCountDetAndfCount(idect,1);
 665     ReconstructFitRowMinRowMax(idect,1);
 666     // Take the histo
 667     fEntriesCurrent = 0;
 668     if(!calvect->GetPHEntries(fCountDet)) {
 669       NotEnoughStatisticPH(idect,fEntriesCurrent);
 670       continue;
 671     }
 672     TString tname("PH");
 673     tname += idect;
 674     TH1F *projph  = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
 675     projph->SetDirectory(0);
 676     if(fEntriesCurrent > 0) fNumberEnt++;
 677     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
 678     // This detector has not enough statistics or was off
 679     if (fEntriesCurrent <=  fMinEntries) {
 680       //printf("Not enough stat!\n");
 681       NotEnoughStatisticPH(idect,fEntriesCurrent);
 682       continue;
 683     }
 684     // Statistic of the histos fitted
 685     fNumberFit++;
 686     fStatisticMean += fEntriesCurrent;
 687     // Calcul of "real" coef
 688     CalculVdriftCoefMean();
 689     CalculT0CoefMean();
 690     //Method choosen
 691     switch(fMethod)
 692       {
 693       case 0: FitLagrangePoly((TH1 *) projph); break;
 694       case 1: FitPente((TH1 *) projph); break;
 695       case 2: FitPH((TH1 *) projph,(Int_t) (idect - fDect1)); break;
 696       default: return kFALSE;
 697       }
 698     // Fill the tree if end of a detector or only the pointer to the branch!!!
 699     FillInfosFitPH(idect,fEntriesCurrent);
 700   } // Boucle object
 701
 702   // Mean Statistic
 703   if (fNumberFit > 0) {
 704     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));
 705     fStatisticMean = fStatisticMean / fNumberFit;
 706   }
 707   else {
 708     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 709   }
 710   delete fDebugStreamer;
 711   fDebugStreamer = 0x0;
 712   return kTRUE;
 713 }
 714 //____________Functions fit Online PRF2d_______________________________________
 715 Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf)
 716 {
 717   //
 718   // Take the 1D profiles (pad response function), projections of the 2D PRF
 719   // on the Xaxis, for each calibration group
 720   // Fit with a gaussian to reconstruct the sigma of the pad response function
 721   //
 722
 723   // Set the calibration mode
 724   const char *name = prf->GetTitle();
 725   if(!SetModeCalibration(name,2)) return kFALSE;
 726
 727   // Number of Ybins (detectors or groups of pads)
 728   Int_t    nybins  = prf->GetNbinsY();// calibration groups
 729   Int_t    nbins   = prf->GetNbinsX();// bins
 730   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)prf->GetTitle());
 731   if((nbg > 0) || (nbg == -1)) return kFALSE;
 732   if (!InitFit(nybins,2)) {
 733     return kFALSE;
 734   }
 735   if (!InitFitPRF()) {
 736     return kFALSE;
 737   }
 738   fStatisticMean        = 0.0;
 739   fNumberFit            = 0;
 740   fNumberFitSuccess     = 0;
 741   fNumberEnt            = 0;
 742   // Init fCountDet and fCount
 743   InitfCountDetAndfCount(2);
 744   // Beginning of the loop
 745   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 746     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
 747     UpdatefCountDetAndfCount(idect,2);
 748     ReconstructFitRowMinRowMax(idect,2);
 749     // Take the histo
 750     TH1D *projprf = (TH1D *) prf->ProjectionX("projprf",idect+1,idect+1,(Option_t *) "e");
 751     projprf->SetDirectory(0);
 752     // Number of entries for this calibration group
 753     Double_t nentries = 0;
 754     for (Int_t k = 0; k < nbins; k++) {
 755       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 756       nentries += prf->GetBinEntries(binnb);
 757     }
 758     if(nentries > 0) fNumberEnt++;
 759     // This detector has not enough statistics or was off
 760     if (nentries <= fMinEntries) {
 761       NotEnoughStatisticPRF(idect);
 762       if (fDebugLevel != 1) {
 763         delete projprf;
 764       }
 765       continue;
 766     }
 767     // Statistics of the histos fitted
 768     fNumberFit++;
 769     fStatisticMean += nentries;
 770     // Calcul of "real" coef
 771     CalculPRFCoefMean();
 772     //Method choosen
 773     switch(fMethod)
 774       {
 775       case 0: FitPRF((TH1 *) projprf); break;
 776       case 1: RmsPRF((TH1 *) projprf); break;
 777       default: return kFALSE;
 778       }
 779     // Fill the tree if end of a detector or only the pointer to the branch!!!
 780     FillInfosFitPRF(idect);
 781     // Memory!!!
 782     if (fDebugLevel != 1) {
 783       delete projprf;
 784     }
 785   } // Boucle object
 786   // Mean Statistic
 787   if (fNumberFit > 0) {
 788     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
 789     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
 790     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
 791                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 792     fStatisticMean = fStatisticMean / fNumberFit;
 793   }
 794   else {
 795     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 796   }
 797   delete fDebugStreamer;
 798   fDebugStreamer = 0x0;
 799   return kTRUE;
 800 }
 801 //____________Functions fit Online PRF2d_______________________________________
 802 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf)
 803 {
 804   //
 805   // Take the 1D profiles (pad response function), projections of the 2D PRF
 806   // on the Xaxis, for each calibration group
 807   // Fit with a gaussian to reconstruct the sigma of the pad response function
 808   //
 809
 810   // Set the calibration mode
 811   const char *name = prf->GetTitle();
 812   if(!SetModeCalibration(name,2)) return kFALSE;
 813
 814   // Number of Ybins (detectors or groups of pads)
 815   TAxis   *xprf    = prf->GetXaxis();
 816   TAxis   *yprf    = prf->GetYaxis();
 817   Int_t    nybins  = yprf->GetNbins();// calibration groups
 818   Int_t    nbins   = xprf->GetNbins();// bins
 819   Float_t  lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
 820   Float_t  upedge  = (Float_t) xprf->GetBinUpEdge(nbins);//upedge in bins
 821   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
 822   if(nbg == -1) return kFALSE;
 823   if(nbg > 0) fMethod = 1;
 824   else fMethod = 0;
 825   if (!InitFit(nybins,2)) {
 826     return kFALSE;
 827   }
 828   if (!InitFitPRF()) {
 829     return kFALSE;
 830   }
 831   fStatisticMean        = 0.0;
 832   fNumberFit            = 0;
 833   fNumberFitSuccess     = 0;
 834   fNumberEnt            = 0;
 835   // Init fCountDet and fCount
 836   InitfCountDetAndfCount(2);
 837   // Beginning of the loop
 838   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 839     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
 840     UpdatefCountDetAndfCount(idect,2);
 841     ReconstructFitRowMinRowMax(idect,2);
 842     // Build the array of entries and sum
 843     TArrayD arraye  = TArrayD(nbins);
 844     TArrayD arraym  = TArrayD(nbins);
 845     TArrayD arrayme = TArrayD(nbins);
 846     Double_t nentries = 0;
 847     //printf("nbins %d\n",nbins);
 848     for (Int_t k = 0; k < nbins; k++) {
 849       Int_t binnb = (nbins+2)*(idect+1)+(k+1);
 850       Double_t entries = (Double_t)prf->GetBinEntries(binnb);
 851       Double_t mean    = (Double_t)prf->GetBinContent(binnb);
 852       Double_t error   = (Double_t)prf->GetBinError(binnb);
 853       //printf("for %d we have %f\n",k,entries);
 854       nentries += entries;
 855       arraye.AddAt(entries,k);
 856       arraym.AddAt(mean,k);
 857       arrayme.AddAt(error,k);
 858     }
 859     if(nentries > 0) fNumberEnt++;
 860     //printf("The number of entries for the group %d is %f\n",idect,nentries);
 861     // This detector has not enough statistics or was off
 862     if (nentries <= fMinEntries) {
 863       NotEnoughStatisticPRF(idect);
 864       continue;
 865     }
 866     // Statistics of the histos fitted
 867     fNumberFit++;
 868     fStatisticMean += nentries;
 869     // Calcul of "real" coef
 870     CalculPRFCoefMean();
 871     //Method choosen
 872     switch(fMethod)
 873       {
 874       case 0: FitPRFGausMI( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbins, lowedge, upedge); break;
 875       case 1: FitTnpRange( arraye.GetArray(), arraym.GetArray(), arrayme.GetArray(), nbg, nbins); break;
 876       default: return kFALSE;
 877       }
 878     // Fill the tree if end of a detector or only the pointer to the branch!!!
 879     FillInfosFitPRF(idect);
 880   } // Boucle object
 881   // Mean Statistic
 882   if (fNumberFit > 0) {
 883     AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
 884     AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
 885     AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
 886                 ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
 887     fStatisticMean = fStatisticMean / fNumberFit;
 888   }
 889   else {
 890     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 891   }
 892   delete fDebugStreamer;
 893   fDebugStreamer = 0x0;
 894   return kTRUE;
 895 }
 896 //____________Functions fit Online PRF2d_______________________________________
 897 Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect)
 898 {
 899   //
 900   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
 901   // each calibration group
 902   // Fit with a gaussian to reconstruct the sigma of the pad response function
 903   //
 904
 905   // Set the calibra mode
 906   const char *name = calvect->GetNamePRF();
 907   if(!SetModeCalibration(name,2)) return kFALSE;
 908   //printf("test0 %s\n",name);
 909
 910   // Number of Xbins (detectors or groups of pads)
 911   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
 912     //printf("test1\n");
 913     return kFALSE;
 914   }
 915   if (!InitFitPRF()) {
 916     ///printf("test2\n");
 917     return kFALSE;
 918   }
 919   fStatisticMean        = 0.0;
 920   fNumberFit            = 0;
 921   fNumberFitSuccess     = 0;
 922   fNumberEnt            = 0;
 923   // Init fCountDet and fCount
 924   InitfCountDetAndfCount(2);
 925   // Beginning of the loop
 926   for (Int_t idect = fDect1; idect < fDect2; idect++) {
 927     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi........
 928     UpdatefCountDetAndfCount(idect,2);
 929     ReconstructFitRowMinRowMax(idect,2);
 930     // Take the histo
 931     fEntriesCurrent = 0;
 932     if(!calvect->GetPRFEntries(fCountDet)) {
 933       NotEnoughStatisticPRF(idect);
 934       continue;
 935     }
 936     TString tname("PRF");
 937     tname += idect;
 938     TH1F *projprf  = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
 939     projprf->SetDirectory(0);
 940     if(fEntriesCurrent > 0) fNumberEnt++;
 941     // This detector has not enough statistics or was off
 942     if (fEntriesCurrent <= fMinEntries) {
 943       NotEnoughStatisticPRF(idect);
 944       continue;
 945     }
 946     // Statistic of the histos fitted
 947     fNumberFit++;
 948     fStatisticMean += fEntriesCurrent;
 949     // Calcul of "real" coef
 950     CalculPRFCoefMean();
 951     //Method choosen
 952     switch(fMethod)
 953       {
 954       case 1: FitPRF((TH1 *) projprf); break;
 955       case 2: RmsPRF((TH1 *) projprf); break;
 956       default: return kFALSE;
 957       }
 958     // Fill the tree if end of a detector or only the pointer to the branch!!!
 959     FillInfosFitPRF(idect);
 960   } // Boucle object
 961   // Mean Statistics
 962   if (fNumberFit > 0) {
 963     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));
 964   }
 965   else {
 966     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
 967   }
 968   delete fDebugStreamer;
 969   fDebugStreamer = 0x0;
 970   return kTRUE;
 971 }
 972 //____________Functions fit Online PRF2d_______________________________________
 973 Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect)
 974 {
 975   //
 976   // Reconstruct the 1D histo (pad response function) from the vectorPRD for
 977   // each calibration group
 978   // Fit with a gaussian to reconstruct the sigma of the pad response function
 979   //
 980
 981   // Set the calibra mode
 982   const char *name = calvect->GetNamePRF();
 983   if(!SetModeCalibration(name,2)) return kFALSE;
 984   //printf("test0 %s\n",name);
 985   Int_t    nbg     = GetNumberOfGroupsPRF((const char *)name);
 986   //printf("test1 %d\n",nbg);
 987   if(nbg == -1) return kFALSE;
 988   if(nbg > 0) fMethod = 1;
 989   else fMethod = 0;
 990   // Number of Xbins (detectors or groups of pads)
 991   if (!InitFit((432*calvect->GetDetCha0(2)+108*calvect->GetDetCha2(2)),2)) {
 992     //printf("test2\n");
 993     return kFALSE;
 994   }
 995   if (!InitFitPRF()) {
 996     //printf("test3\n");
 997     return kFALSE;
 998   }
 999   fStatisticMean        = 0.0;
1000   fNumberFit            = 0;
1001   fNumberFitSuccess     = 0;
1002   fNumberEnt            = 0;
1003   // Variables
1004   Int_t nbins           = 0;
1005   Double_t *arrayx       = 0;
1006   Double_t *arraye       = 0;
1007   Double_t *arraym       = 0;
1008   Double_t *arrayme      = 0;
1009   Float_t lowedge       = 0.0;
1010   Float_t upedge        = 0.0;
1011   // Init fCountDet and fCount
1012   InitfCountDetAndfCount(2);
1013   // Beginning of the loop
1014   for (Int_t idect = fDect1; idect < fDect2; idect++) {
1015     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi......
1016     UpdatefCountDetAndfCount(idect,2);
1017     ReconstructFitRowMinRowMax(idect,2);
1018     // Take the histo
1019     fEntriesCurrent  = 0;
1020     if(!calvect->GetPRFEntries(fCountDet)) {
1021       NotEnoughStatisticPRF(idect);
1022       continue;
1023     }
1024     TString tname("PRF");
1025     tname += idect;
1026     TGraphErrors *projprftree  = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname);
1027     nbins   = projprftree->GetN();
1028     arrayx  = (Double_t *)projprftree->GetX();
1029     arraye  = (Double_t *)projprftree->GetEX();
1030     arraym  = (Double_t *)projprftree->GetY();
1031     arrayme = (Double_t *)projprftree->GetEY();
1032     Float_t step = arrayx[1]-arrayx[0];
1033     lowedge = arrayx[0] - step/2.0;
1034     upedge  = arrayx[(nbins-1)] + step/2.0;
1035     //printf("nbins est %d\n",nbins);
1036     for(Int_t k = 0; k < nbins; k++){
1037       fEntriesCurrent += (Int_t)arraye[k];
1038       //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k]));
1039       if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]);
1040     }
1041     if(fEntriesCurrent > 0) fNumberEnt++;
1042     //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
1043     // This detector has not enough statistics or was off
1044     if (fEntriesCurrent <= fMinEntries) {
1045       NotEnoughStatisticPRF(idect);
1046       continue;
1047     }
1048     // Statistic of the histos fitted
1049     fNumberFit++;
1050     fStatisticMean += fEntriesCurrent;
1051     // Calcul of "real" coef
1052     CalculPRFCoefMean();
1053     //Method choosen
1054     switch(fMethod)
1055       {
1056       case 0: FitPRFGausMI(arraye,arraym,arrayme,nbins,lowedge,upedge); break;
1057       case 1: FitTnpRange(arraye,arraym,arrayme,nbg,nbins); break;
1058       default: return kFALSE;
1059       }
1060     // Fill the tree if end of a detector or only the pointer to the branch!!!
1061     FillInfosFitPRF(idect);
1062   } // Boucle object
1063   // Mean Statistics
1064   if (fNumberFit > 0) {
1065     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));
1066   }
1067   else {
1068     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1069   }
1070   delete fDebugStreamer;
1071   fDebugStreamer = 0x0;
1072   return kTRUE;
1073 }
1074 //____________Functions fit Online CH2d________________________________________
1075 Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *calivdli)
1076 {
1077   //
1078   // The linear method
1079   //
1080
1081   fStatisticMean        = 0.0;
1082   fNumberFit            = 0;
1083   fNumberFitSuccess     = 0;
1084   fNumberEnt            = 0;
1085   if(!InitFitLinearFitter()) return kFALSE;
1086
1087
1088   for(Int_t idet = 0; idet < 540; idet++){
1089
1090
1091     //printf("detector number %d\n",idet);
1092
1093     // Take the result
1094     TVectorD param(2);
1095     TVectorD error(3);
1096     fEntriesCurrent = 0;
1097     fCountDet       = idet;
1098     Bool_t here     = calivdli->GetParam(idet,&param);
1099     Bool_t heree    = calivdli->GetError(idet,&error);
1100     //printf("here %d and heree %d\n",here, heree);
1101     if(heree) {
1102       fEntriesCurrent = (Int_t) error[2];
1103       fNumberEnt++;
1104     }
1105     //printf("Number of entries %d\n",fEntriesCurrent);
1106     // Nothing found or not enough statistic
1107     if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
1108       NotEnoughStatisticLinearFitter();
1109       continue;
1110     }
1111     //param.Print();
1112     //error.Print();
1113     //Statistics
1114     fNumberFit++;
1115     fStatisticMean += fEntriesCurrent;
1116
1117     // Check the fit
1118     if((-(param[1])) <= 0.0) {
1119       NotEnoughStatisticLinearFitter();
1120       continue;
1121     }
1122
1123     // CalculDatabaseVdriftandTan
1124     CalculVdriftLorentzCoef();
1125
1126     // Statistics
1127     fNumberFitSuccess ++;
1128
1129     // Put the fCurrentCoef
1130     fCurrentCoef[0]  = -param[1];
1131     // here the database must be the one of the reconstruction for the lorentz angle....
1132     fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
1133     fCurrentCoefE    = error[1];
1134     fCurrentCoefE2   = error[0];
1135     if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
1136       fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
1137     }
1138
1139     // Fill
1140     FillInfosFitLinearFitter();
1141
1142
1143   }
1144   // Mean Statistics
1145   if (fNumberFit > 0) {
1146     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));
1147   }
1148   else {
1149     AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
1150   }
1151   delete fDebugStreamer;
1152   fDebugStreamer = 0x0;
1153   return kTRUE;
1154
1155 }
1156 //____________Functions for seeing if the pad is really okey___________________
1157 //_____________________________________________________________________________
1158 Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
1159 {
1160   //
1161   // Get numberofgroupsprf
1162   //
1163
1164   // Some patterns
1165   const Char_t *pattern0 = "Ngp0";
1166   const Char_t *pattern1 = "Ngp1";
1167   const Char_t *pattern2 = "Ngp2";
1168   const Char_t *pattern3 = "Ngp3";
1169   const Char_t *pattern4 = "Ngp4";
1170   const Char_t *pattern5 = "Ngp5";
1171   const Char_t *pattern6 = "Ngp6";
1172
1173   // Nrphi mode
1174   if (strstr(nametitle,pattern0)) {
1175     return 0;
1176   }
1177   if (strstr(nametitle,pattern1)) {
1178     return 1;
1179   }
1180   if (strstr(nametitle,pattern2)) {
1181     return 2;
1182   }
1183   if (strstr(nametitle,pattern3)) {
1184     return 3;
1185   }
1186   if (strstr(nametitle,pattern4)) {
1187     return 4;
1188   }
1189   if (strstr(nametitle,pattern5)) {
1190     return 5;
1191   }
1192   if (strstr(nametitle,pattern6)){
1193     return 6;
1194   }
1195   else return -1;
1196
1197
1198 }
1199 //_____________________________________________________________________________
1200 Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
1201 {
1202   //
1203   // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
1204   // corresponding to the given name
1205   //
1206
1207   if(!SetNzFromTObject(name,i)) return kFALSE;
1208   if(!SetNrphiFromTObject(name,i)) return kFALSE;
1209
1210   return kTRUE;
1211
1212 }
1213 //_____________________________________________________________________________
1214 Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
1215 {
1216   //
1217   // Set fNrphi[i] of the AliTRDCalibraFit::Instance()
1218   // corresponding to the given TObject
1219   //
1220
1221   // Some patterns
1222   const Char_t *patternrphi0 = "Nrphi0";
1223   const Char_t *patternrphi1 = "Nrphi1";
1224   const Char_t *patternrphi2 = "Nrphi2";
1225   const Char_t *patternrphi3 = "Nrphi3";
1226   const Char_t *patternrphi4 = "Nrphi4";
1227   const Char_t *patternrphi5 = "Nrphi5";
1228   const Char_t *patternrphi6 = "Nrphi6";
1229
1230
1231   const Char_t *patternrphi10 = "Nrphi10";
1232   const Char_t *patternrphi100 = "Nrphi100";
1233   const Char_t *patternz10 = "Nz10";
1234   const Char_t *patternz100 = "Nz100";
1235
1236   // Nrphi mode
1237   if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1238     fCalibraMode->SetAllTogether(i);
1239     fNbDet = 540;
1240     if (fDebugLevel > 1) {
1241       AliInfo(Form("fNbDet %d and 100",fNbDet));
1242     }
1243     return kTRUE;
1244   }
1245   if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1246     fCalibraMode->SetPerSuperModule(i);
1247     fNbDet = 30;
1248     if (fDebugLevel > 1) {
1249       AliInfo(Form("fNDet %d and 100",fNbDet));
1250     }
1251     return kTRUE;
1252   }
1253
1254   if (strstr(name,patternrphi0)) {
1255     fCalibraMode->SetNrphi(i ,0);
1256     if (fDebugLevel > 1) {
1257       AliInfo(Form("fNbDet %d and 0",fNbDet));
1258     }
1259     return kTRUE;
1260   }
1261   if (strstr(name,patternrphi1)) {
1262     fCalibraMode->SetNrphi(i, 1);
1263     if (fDebugLevel > 1) {
1264       AliInfo(Form("fNbDet %d and 1",fNbDet));
1265     }
1266     return kTRUE;
1267   }
1268   if (strstr(name,patternrphi2)) {
1269     fCalibraMode->SetNrphi(i, 2);
1270     if (fDebugLevel > 1) {
1271       AliInfo(Form("fNbDet %d and 2",fNbDet));
1272     }
1273     return kTRUE;
1274   }
1275   if (strstr(name,patternrphi3)) {
1276     fCalibraMode->SetNrphi(i, 3);
1277     if (fDebugLevel > 1) {
1278       AliInfo(Form("fNbDet %d and 3",fNbDet));
1279     }
1280     return kTRUE;
1281   }
1282   if (strstr(name,patternrphi4)) {
1283     fCalibraMode->SetNrphi(i, 4);
1284     if (fDebugLevel > 1) {
1285       AliInfo(Form("fNbDet %d and 4",fNbDet));
1286     }
1287     return kTRUE;
1288   }
1289   if (strstr(name,patternrphi5)) {
1290     fCalibraMode->SetNrphi(i, 5);
1291     if (fDebugLevel > 1) {
1292       AliInfo(Form("fNbDet %d and 5",fNbDet));
1293     }
1294     return kTRUE;
1295   }
1296   if (strstr(name,patternrphi6)) {
1297     fCalibraMode->SetNrphi(i, 6);
1298     if (fDebugLevel > 1) {
1299       AliInfo(Form("fNbDet %d and 6",fNbDet));
1300     }
1301     return kTRUE;
1302   }
1303
1304   if (fDebugLevel > 1) {
1305     AliInfo(Form("fNbDet %d and rest",fNbDet));
1306   }
1307   fCalibraMode->SetNrphi(i ,0);
1308   return kFALSE;
1309
1310 }
1311 //_____________________________________________________________________________
1312 Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
1313 {
1314   //
1315   // Set fNz[i] of the AliTRDCalibraFit::Instance()
1316   // corresponding to the given TObject
1317   //
1318
1319   // Some patterns
1320   const Char_t *patternz0    = "Nz0";
1321   const Char_t *patternz1    = "Nz1";
1322   const Char_t *patternz2    = "Nz2";
1323   const Char_t *patternz3    = "Nz3";
1324   const Char_t *patternz4    = "Nz4";
1325
1326   const Char_t *patternrphi10 = "Nrphi10";
1327   const Char_t *patternrphi100 = "Nrphi100";
1328   const Char_t *patternz10 = "Nz10";
1329   const Char_t *patternz100 = "Nz100";
1330
1331   if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1332     fCalibraMode->SetAllTogether(i);
1333     fNbDet = 540;
1334     if (fDebugLevel > 1) {
1335       AliInfo(Form("fNbDet %d and 100",fNbDet));
1336     }
1337     return kTRUE;
1338   }
1339   if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1340     fCalibraMode->SetPerSuperModule(i);
1341     fNbDet = 30;
1342     if (fDebugLevel > 1) {
1343       AliInfo(Form("fNbDet %d and 10",fNbDet));
1344     }
1345     return kTRUE;
1346   }
1347   if (strstr(name,patternz0)) {
1348     fCalibraMode->SetNz(i, 0);
1349     if (fDebugLevel > 1) {
1350       AliInfo(Form("fNbDet %d and 0",fNbDet));
1351     }
1352     return kTRUE;
1353   }
1354   if (strstr(name,patternz1)) {
1355     fCalibraMode->SetNz(i ,1);
1356     if (fDebugLevel > 1) {
1357       AliInfo(Form("fNbDet %d and 1",fNbDet));
1358     }
1359     return kTRUE;
1360   }
1361   if (strstr(name,patternz2)) {
1362     fCalibraMode->SetNz(i ,2);
1363     if (fDebugLevel > 1) {
1364       AliInfo(Form("fNbDet %d and 2",fNbDet));
1365     }
1366     return kTRUE;
1367   }
1368   if (strstr(name,patternz3)) {
1369     fCalibraMode->SetNz(i ,3);
1370     if (fDebugLevel > 1) {
1371       AliInfo(Form("fNbDet %d and 3",fNbDet));
1372     }
1373     return kTRUE;
1374   }
1375   if (strstr(name,patternz4)) {
1376     fCalibraMode->SetNz(i ,4);
1377     if (fDebugLevel > 1) {
1378       AliInfo(Form("fNbDet %d and 4",fNbDet));
1379     }
1380     return kTRUE;
1381   }
1382
1383   if (fDebugLevel > 1) {
1384     AliInfo(Form("fNbDet %d and rest",fNbDet));
1385   }
1386   fCalibraMode->SetNz(i ,0);
1387   return kFALSE;
1388 }
1389 //______________________________________________________________________
1390 void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
1391   //
1392   // ofwhat is equaled to 0: mean value of all passing detectors
1393   // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1394   //
1395
1396   Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
1397   if(loop != 540) {
1398     AliInfo("The Vector Fit is not complete!");
1399     return;
1400   }
1401   Int_t detector = -1;
1402   Int_t sector = -1;
1403   Float_t value  = 0.0;
1404
1405   /////////////////////////////////
1406   // Calculate the mean values
1407   ////////////////////////////////
1408   // Initialisation
1409   ////////////////////////
1410   Double_t meanAll = 0.0;
1411   Double_t meanSupermodule[18];
1412   Double_t meanDetector[540];
1413   Int_t countAll = 0;
1414   Int_t countSupermodule[18];
1415   Int_t countDetector[540];
1416   for(Int_t sm = 0; sm < 18; sm++){
1417     meanSupermodule[sm] = 0.0;
1418     countSupermodule[sm] = 0;
1419   }
1420   for(Int_t det = 0; det < 540; det++){
1421     meanDetector[det] = 0.0;
1422     countDetector[det] = 0;
1423   }
1424   // compute
1425   ////////////
1426   for (Int_t k = 0; k < loop; k++) {
1427     detector  = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1428     sector = GetSector(detector);
1429     if(perdetector){
1430       value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
1431       if(value > 0.0) {
1432         meanDetector[detector] += value;
1433         countDetector[detector]++;
1434         meanSupermodule[sector] += value;
1435         countSupermodule[sector]++;
1436         meanAll += value;
1437         countAll++;
1438       }
1439     }
1440     else {
1441       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1442       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1443       for (Int_t row = 0; row < rowMax; row++) {
1444         for (Int_t col = 0; col < colMax; col++) {
1445           value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1446           if(value > 0.0) {
1447             meanDetector[detector] += value;
1448             countDetector[detector]++;
1449             meanSupermodule[sector] += value;
1450             countSupermodule[sector]++;
1451             meanAll += value;
1452             countAll++;
1453           }
1454
1455         } // Col
1456       } // Row
1457     }
1458   }
1459   if(countAll > 0) meanAll = meanAll/countAll;
1460   for(Int_t sm = 0; sm < 18; sm++){
1461     if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1462   }
1463   for(Int_t det = 0; det < 540; det++){
1464     if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1465   }
1466   // Put the mean value for the no-fitted
1467   /////////////////////////////////////////////
1468   for (Int_t k = 0; k < loop; k++) {
1469     detector  = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
1470     sector = GetSector(detector);
1471     Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1472     Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1473     Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
1474
1475     for (Int_t row = 0; row < rowMax; row++) {
1476       for (Int_t col = 0; col < colMax; col++) {
1477         value = coef[(Int_t)(col*rowMax+row)];
1478         if(value < 0.0) {
1479           if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1480           if(ofwhat == 1){
1481             if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
1482             else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
1483             else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
1484           }
1485         }
1486         // Debug
1487         if(fDebugLevel > 1){
1488
1489           if ( !fDebugStreamer ) {
1490             //debug stream
1491             TDirectory *backup = gDirectory;
1492             fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1493             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1494           }
1495
1496           Float_t coefnow      = coef[(Int_t)(col*rowMax+row)];
1497
1498           (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<<
1499             "detector="<<detector<<
1500             "sector="<<sector<<
1501             "row="<<row<<
1502             "col="<<col<<
1503             "before="<<value<<
1504             "after="<<coefnow<<
1505             "\n";
1506         }
1507       } // Col
1508     } // Row
1509   }
1510
1511 }
1512 //______________________________________________________________________
1513 void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
1514   //
1515   // ofwhat is equaled to 0: mean value of all passing detectors
1516   // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all
1517   //
1518
1519   Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
1520   if(loop != 540) {
1521     AliInfo("The Vector Fit is not complete!");
1522     return;
1523   }
1524   Int_t detector = -1;
1525   Int_t sector = -1;
1526   Float_t value  = 0.0;
1527
1528   /////////////////////////////////
1529   // Calculate the mean values
1530   ////////////////////////////////
1531   // Initialisation
1532   ////////////////////////
1533   Double_t meanAll = 0.0;
1534   Double_t meanSupermodule[18];
1535   Double_t meanDetector[540];
1536   Int_t countAll = 0;
1537   Int_t countSupermodule[18];
1538   Int_t countDetector[540];
1539   for(Int_t sm = 0; sm < 18; sm++){
1540     meanSupermodule[sm] = 0.0;
1541     countSupermodule[sm] = 0;
1542   }
1543   for(Int_t det = 0; det < 540; det++){
1544     meanDetector[det] = 0.0;
1545     countDetector[det] = 0;
1546   }
1547   // compute
1548   ////////////
1549   for (Int_t k = 0; k < loop; k++) {
1550     detector  = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1551     sector = GetSector(detector);
1552     if(perdetector){
1553       value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
1554       if(value < 70.0) {
1555         meanDetector[detector] += value;
1556         countDetector[detector]++;
1557         meanSupermodule[sector] += value;
1558         countSupermodule[sector]++;
1559         meanAll += value;
1560         countAll++;
1561       }
1562     }
1563     else {
1564       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1565       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1566       for (Int_t row = 0; row < rowMax; row++) {
1567         for (Int_t col = 0; col < colMax; col++) {
1568           value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1569           if(value < 70.0) {
1570             meanDetector[detector] += value;
1571             countDetector[detector]++;
1572             meanSupermodule[sector] += value;
1573             countSupermodule[sector]++;
1574             meanAll += value;
1575             countAll++;
1576           }
1577
1578         } // Col
1579       } // Row
1580     }
1581   }
1582   if(countAll > 0) meanAll = meanAll/countAll;
1583   for(Int_t sm = 0; sm < 18; sm++){
1584     if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
1585   }
1586   for(Int_t det = 0; det < 540; det++){
1587     if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
1588   }
1589   // Put the mean value for the no-fitted
1590   /////////////////////////////////////////////
1591   for (Int_t k = 0; k < loop; k++) {
1592     detector  = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
1593     sector = GetSector(detector);
1594     Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1595     Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1596     Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
1597
1598     for (Int_t row = 0; row < rowMax; row++) {
1599       for (Int_t col = 0; col < colMax; col++) {
1600         value = coef[(Int_t)(col*rowMax+row)];
1601         if(value > 70.0) {
1602           if((ofwhat == 0) && (meanAll > -1.5)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1603           if(ofwhat == 1){
1604             if(meanDetector[detector] > -1.5) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
1605             else if(meanSupermodule[sector] > -1.5) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
1606             else if(meanAll > -1.5) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
1607           }
1608         }
1609         // Debug
1610         if(fDebugLevel > 1){
1611
1612           if ( !fDebugStreamer ) {
1613             //debug stream
1614             TDirectory *backup = gDirectory;
1615             fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root");
1616             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1617           }
1618
1619           Float_t coefnow      = coef[(Int_t)(col*rowMax+row)];
1620
1621           (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<<
1622             "detector="<<detector<<
1623             "sector="<<sector<<
1624             "row="<<row<<
1625             "col="<<col<<
1626             "before="<<value<<
1627             "after="<<coefnow<<
1628             "\n";
1629         }
1630       } // Col
1631     } // Row
1632   }
1633
1634 }
1635 //_____________________________________________________________________________
1636 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(const TObjArray *vectorFit, Bool_t perdetector)
1637 {
1638   //
1639   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1640   // It takes the mean value of the coefficients per detector
1641   // This object has to be written in the database
1642   //
1643
1644   // Create the DetObject
1645   AliTRDCalDet *object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
1646
1647   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1648   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1649   Int_t detector = -1;
1650   Float_t value  = 0.0;
1651
1652   //
1653   for (Int_t k = 0; k < loop; k++) {
1654     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1655     Float_t mean  = 0.0;
1656     if(perdetector){
1657       mean = TMath::Abs(((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]);
1658     }
1659     else {
1660       Int_t   count = 0;
1661       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1662       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1663       for (Int_t row = 0; row < rowMax; row++) {
1664         for (Int_t col = 0; col < colMax; col++) {
1665           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1666           mean += TMath::Abs(value);
1667           count++;
1668         } // Col
1669       } // Row
1670       if(count > 0) mean = mean/count;
1671     }
1672     object->SetValue(detector,mean);
1673   }
1674
1675   return object;
1676 }
1677 //_____________________________________________________________________________
1678 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector)
1679 {
1680   //
1681   // It creates the AliTRDCalDet object from the AliTRDFitInfo
1682   // It takes the mean value of the coefficients per detector
1683   // This object has to be written in the database
1684   //
1685
1686   // Create the DetObject
1687   AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
1688
1689
1690   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1691   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1692   Int_t detector = -1;
1693   Float_t value  = 0.0;
1694
1695   for (Int_t k = 0; k < loop; k++) {
1696     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1697     Float_t mean  = 0.0;
1698     if(perdetector){
1699       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1700       if(!meanOtherBefore){
1701         if(value > 0) value = value*scaleFitFactor;
1702       }
1703       else value = value*scaleFitFactor;
1704       mean = TMath::Abs(value);
1705     }
1706     else{
1707       Int_t   count = 0;
1708       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1709       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1710       for (Int_t row = 0; row < rowMax; row++) {
1711         for (Int_t col = 0; col < colMax; col++) {
1712           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1713           if(!meanOtherBefore) {
1714             if(value > 0) value = value*scaleFitFactor;
1715           }
1716           else value = value*scaleFitFactor;
1717           mean += TMath::Abs(value);
1718           count++;
1719         } // Col
1720       } // Row
1721       if(count > 0) mean = mean/count;
1722     }
1723     object->SetValue(detector,mean);
1724   }
1725
1726   return object;
1727 }
1728 //_____________________________________________________________________________
1729 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1730 {
1731   //
1732   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1733   // It takes the min value of the coefficients per detector
1734   // This object has to be written in the database
1735   //
1736
1737   // Create the DetObject
1738   AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1739
1740   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1741   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1742   Int_t detector = -1;
1743   Float_t value  = 0.0;
1744
1745   for (Int_t k = 0; k < loop; k++) {
1746     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1747     Float_t min  = 100.0;
1748     if(perdetector){
1749       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1750       // check successful
1751       if(value > 70.0) value = value-100.0;
1752       //
1753       min = value;
1754     }
1755     else{
1756       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1757       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1758       for (Int_t row = 0; row < rowMax; row++) {
1759         for (Int_t col = 0; col < colMax; col++) {
1760           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1761           // check successful
1762           if(value > 70.0) value = value-100.0;
1763           //
1764           if(min > value) min = value;
1765         } // Col
1766       } // Row
1767     }
1768     object->SetValue(detector,min);
1769   }
1770
1771   return object;
1772
1773 }
1774 //_____________________________________________________________________________
1775 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
1776 {
1777   //
1778   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1779   // It takes the min value of the coefficients per detector
1780   // This object has to be written in the database
1781   //
1782
1783   // Create the DetObject
1784   AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
1785
1786
1787   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1788   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1789   Int_t detector = -1;
1790   Float_t value  = 0.0;
1791
1792   for (Int_t k = 0; k < loop; k++) {
1793     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1794     /*
1795       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1796       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1797       Float_t min  = 100.0;
1798       for (Int_t row = 0; row < rowMax; row++) {
1799       for (Int_t col = 0; col < colMax; col++) {
1800       value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1801       mean += -TMath::Abs(value);
1802       count++;
1803       } // Col
1804       } // Row
1805       if(count > 0) mean = mean/count;
1806     */
1807     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1808     object->SetValue(detector,-TMath::Abs(value));
1809   }
1810
1811   return object;
1812
1813 }
1814 //_____________________________________________________________________________
1815 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
1816 {
1817   //
1818   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1819   // You need first to create the object for the detectors,
1820   // where the mean value is put.
1821   // This object has to be written in the database
1822   //
1823
1824   // Create the DetObject
1825   AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
1826
1827   if(!vectorFit){
1828     for(Int_t k = 0; k < 540; k++){
1829       AliTRDCalROC *calROC = object->GetCalROC(k);
1830       Int_t nchannels = calROC->GetNchannels();
1831       for(Int_t ch = 0; ch < nchannels; ch++){
1832         calROC->SetValue(ch,1.0);
1833       }
1834     }
1835   }
1836   else{
1837
1838     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1839     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1840     Int_t detector = -1;
1841     Float_t value  = 0.0;
1842
1843     for (Int_t k = 0; k < loop; k++) {
1844       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1845       AliTRDCalROC *calROC = object->GetCalROC(detector);
1846       Float_t mean         = detobject->GetValue(detector);
1847       if(mean == 0) continue;
1848       Int_t rowMax    = calROC->GetNrows();
1849       Int_t colMax    = calROC->GetNcols();
1850       for (Int_t row = 0; row < rowMax; row++) {
1851         for (Int_t col = 0; col < colMax; col++) {
1852           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1853           if(value > 0) value = value*scaleFitFactor;
1854           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1855         } // Col
1856       } // Row
1857     }
1858   }
1859
1860   return object;
1861 }
1862 //_____________________________________________________________________________
1863 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1864 {
1865   //
1866   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1867   // You need first to create the object for the detectors,
1868   // where the mean value is put.
1869   // This object has to be written in the database
1870   //
1871
1872   // Create the DetObject
1873   AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
1874
1875   if(!vectorFit){
1876     for(Int_t k = 0; k < 540; k++){
1877       AliTRDCalROC *calROC = object->GetCalROC(k);
1878       Int_t nchannels = calROC->GetNchannels();
1879       for(Int_t ch = 0; ch < nchannels; ch++){
1880         calROC->SetValue(ch,1.0);
1881       }
1882     }
1883   }
1884   else {
1885
1886     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1887     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1888     Int_t detector = -1;
1889     Float_t value  = 0.0;
1890
1891     for (Int_t k = 0; k < loop; k++) {
1892       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1893       AliTRDCalROC *calROC = object->GetCalROC(detector);
1894       Float_t mean         = detobject->GetValue(detector);
1895       if(mean == 0) continue;
1896       Int_t rowMax    = calROC->GetNrows();
1897       Int_t colMax    = calROC->GetNcols();
1898       for (Int_t row = 0; row < rowMax; row++) {
1899         for (Int_t col = 0; col < colMax; col++) {
1900           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1901           calROC->SetValue(col,row,TMath::Abs(value)/mean);
1902         } // Col
1903       } // Row
1904     }
1905   }
1906   return object;
1907
1908 }
1909 //_____________________________________________________________________________
1910 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
1911 {
1912   //
1913   // It Creates the AliTRDCalPad object from AliTRDFitInfo2
1914   // You need first to create the object for the detectors,
1915   // where the mean value is put.
1916   // This object has to be written in the database
1917   //
1918
1919   // Create the DetObject
1920   AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
1921
1922   if(!vectorFit){
1923     for(Int_t k = 0; k < 540; k++){
1924       AliTRDCalROC *calROC = object->GetCalROC(k);
1925       Int_t nchannels = calROC->GetNchannels();
1926       for(Int_t ch = 0; ch < nchannels; ch++){
1927         calROC->SetValue(ch,0.0);
1928       }
1929     }
1930   }
1931   else {
1932
1933     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1934     if(loop != 540) AliInfo("The Vector Fit is not complete!");
1935     Int_t detector = -1;
1936     Float_t value  = 0.0;
1937
1938     for (Int_t k = 0; k < loop; k++) {
1939       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1940       AliTRDCalROC *calROC = object->GetCalROC(detector);
1941       Float_t min          = detobject->GetValue(detector);
1942       Int_t rowMax    = calROC->GetNrows();
1943       Int_t colMax    = calROC->GetNcols();
1944       for (Int_t row = 0; row < rowMax; row++) {
1945         for (Int_t col = 0; col < colMax; col++) {
1946           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1947           // check successful
1948           if(value > 70.0) value = value - 100.0;
1949           //
1950           calROC->SetValue(col,row,value-min);
1951         } // Col
1952       } // Row
1953     }
1954   }
1955   return object;
1956
1957 }
1958 //_____________________________________________________________________________
1959 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
1960 {
1961   //
1962   // It Creates the AliTRDCalPad object from AliTRDFitInfo
1963   // This object has to be written in the database
1964   //
1965
1966   // Create the DetObject
1967   AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
1968
1969   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1970   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1971   Int_t detector = -1;
1972   Float_t value  = 0.0;
1973
1974   for (Int_t k = 0; k < loop; k++) {
1975     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1976     AliTRDCalROC *calROC = object->GetCalROC(detector);
1977     Int_t rowMax    = calROC->GetNrows();
1978     Int_t colMax    = calROC->GetNcols();
1979     for (Int_t row = 0; row < rowMax; row++) {
1980       for (Int_t col = 0; col < colMax; col++) {
1981         value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1982         calROC->SetValue(col,row,TMath::Abs(value));
1983       } // Col
1984     } // Row
1985   }
1986
1987   return object;
1988
1989 }
1990 //_____________________________________________________________________________
1991 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
1992 {
1993   //
1994   // It Creates the AliTRDCalDet object from AliTRDFitInfo
1995   // 0 successful fit 1 not successful fit
1996   // mean is the mean value over the successful fit
1997   // do not use it for t0: no meaning
1998   //
1999
2000   // Create the CalObject
2001   AliTRDCalDet *object = new AliTRDCalDet(name,name);
2002   mean = 0.0;
2003   Int_t count = 0;
2004
2005   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2006   if(loop != 540) {
2007     AliInfo("The Vector Fit is not complete! We initialise all outliers");
2008     for(Int_t k = 0; k < 540; k++){
2009       object->SetValue(k,1.0);
2010     }
2011   }
2012   Int_t detector = -1;
2013   Float_t value  = 0.0;
2014
2015   for (Int_t k = 0; k < loop; k++) {
2016     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2017     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2018     if(value <= 0) object->SetValue(detector,1.0);
2019     else {
2020       object->SetValue(detector,0.0);
2021       mean += value;
2022       count++;
2023     }
2024   }
2025   if(count > 0) mean /= count;
2026   return object;
2027 }
2028 //_____________________________________________________________________________
2029 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2030 {
2031   //
2032   // It Creates the AliTRDCalPad object from AliTRDFitInfo
2033   // 0 not successful fit 1 successful fit
2034   // mean mean value over the successful fit
2035   //
2036
2037   // Create the CalObject
2038   AliTRDCalPad *object = new AliTRDCalPad(name,name);
2039   mean = 0.0;
2040   Int_t count = 0;
2041
2042   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2043   if(loop != 540) {
2044     AliInfo("The Vector Fit is not complete! We initialise all outliers");
2045     for(Int_t k = 0; k < 540; k++){
2046       AliTRDCalROC *calROC = object->GetCalROC(k);
2047       Int_t nchannels = calROC->GetNchannels();
2048       for(Int_t ch = 0; ch < nchannels; ch++){
2049         calROC->SetValue(ch,1.0);
2050       }
2051     }
2052   }
2053   Int_t detector = -1;
2054   Float_t value  = 0.0;
2055
2056   for (Int_t k = 0; k < loop; k++) {
2057     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2058     AliTRDCalROC *calROC = object->GetCalROC(detector);
2059     Int_t nchannels    = calROC->GetNchannels();
2060     for (Int_t ch = 0; ch < nchannels; ch++) {
2061       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2062       if(value <= 0) calROC->SetValue(ch,1.0);
2063       else {
2064         calROC->SetValue(ch,0.0);
2065         mean += value;
2066         count++;
2067       }
2068     } // channels
2069   }
2070   if(count > 0) mean /= count;
2071   return object;
2072 }
2073 //_____________________________________________________________________________
2074 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2075 {
2076   //
2077   // Set FitPH if 1 then each detector will be fitted
2078   //
2079
2080   if (periodeFitPH > 0) {
2081     fFitPHPeriode   = periodeFitPH;
2082   }
2083   else {
2084     AliInfo("periodeFitPH must be higher than 0!");
2085   }
2086
2087 }
2088 //_____________________________________________________________________________
2089 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2090 {
2091   //
2092   // The fit of the deposited charge distribution begins at
2093   // histo->Mean()/beginFitCharge
2094   // You can here set beginFitCharge
2095   //
2096
2097   if (beginFitCharge > 0) {
2098     fBeginFitCharge = beginFitCharge;
2099   }
2100   else {
2101     AliInfo("beginFitCharge must be strict positif!");
2102   }
2103
2104 }
2105
2106 //_____________________________________________________________________________
2107 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2108 {
2109   //
2110   // The t0 calculated with the maximum positif slope is shift from t0Shift0
2111   // You can here set t0Shift0
2112   //
2113
2114   if (t0Shift > 0) {
2115     fT0Shift0 = t0Shift;
2116   }
2117   else {
2118     AliInfo("t0Shift0 must be strict positif!");
2119   }
2120
2121 }
2122
2123 //_____________________________________________________________________________
2124 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2125 {
2126   //
2127   // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2128   // You can here set t0Shift1
2129   //
2130
2131   if (t0Shift > 0) {
2132     fT0Shift1 = t0Shift;
2133   }
2134   else {
2135     AliInfo("t0Shift must be strict positif!");
2136   }
2137
2138 }
2139
2140 //_____________________________________________________________________________
2141 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2142 {
2143   //
2144   // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2145   // You can here set rangeFitPRF
2146   //
2147
2148   if ((rangeFitPRF >    0) &&
2149       (rangeFitPRF <= 1.5)) {
2150     fRangeFitPRF = rangeFitPRF;
2151   }
2152   else {
2153     AliInfo("rangeFitPRF must be between 0 and 1.0");
2154   }
2155
2156 }
2157
2158 //_____________________________________________________________________________
2159 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2160 {
2161   //
2162   // Minimum entries for fitting
2163   //
2164
2165   if (minEntries >    0) {
2166     fMinEntries = minEntries;
2167   }
2168   else {
2169     AliInfo("fMinEntries must be >= 0.");
2170   }
2171
2172 }
2173
2174 //_____________________________________________________________________________
2175 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2176 {
2177   //
2178   // Rebin with rebin time less bins the Ch histo
2179   // You can set here rebin that should divide the number of bins of CH histo
2180   //
2181
2182   if (rebin > 0) {
2183     fRebin = rebin;
2184     AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2185   }
2186   else {
2187     AliInfo("You have to choose a positiv value!");
2188   }
2189
2190 }
2191 //_____________________________________________________________________________
2192 Bool_t AliTRDCalibraFit::FillVectorFit()
2193 {
2194   //
2195   // For the Fit functions fill the vector Fit
2196   //
2197
2198   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2199
2200   Int_t ntotal = 1;
2201   if (GetStack(fCountDet) == 2) {
2202     ntotal = 1728;
2203   }
2204   else {
2205     ntotal = 2304;
2206   }
2207
2208   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2209   Float_t *coef = new Float_t[ntotal];
2210   for (Int_t i = 0; i < ntotal; i++) {
2211     coef[i] = fCurrentCoefDetector[i];
2212   }
2213
2214   Int_t detector = fCountDet;
2215   // Set
2216   fitInfo->SetCoef(coef);
2217   fitInfo->SetDetector(detector);
2218   fVectorFit.Add((TObject *) fitInfo);
2219
2220   return kTRUE;
2221
2222 }
2223 //_____________________________________________________________________________
2224 Bool_t AliTRDCalibraFit::FillVectorFit2()
2225 {
2226   //
2227   // For the Fit functions fill the vector Fit
2228   //
2229
2230   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2231
2232   Int_t ntotal = 1;
2233   if (GetStack(fCountDet) == 2) {
2234     ntotal = 1728;
2235   }
2236   else {
2237     ntotal = 2304;
2238   }
2239
2240   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2241   Float_t *coef = new Float_t[ntotal];
2242   for (Int_t i = 0; i < ntotal; i++) {
2243     coef[i] = fCurrentCoefDetector2[i];
2244   }
2245
2246   Int_t detector = fCountDet;
2247   // Set
2248   fitInfo->SetCoef(coef);
2249   fitInfo->SetDetector(detector);
2250   fVectorFit2.Add((TObject *) fitInfo);
2251
2252   return kTRUE;
2253
2254 }
2255 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2256 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2257 {
2258   //
2259   // Init the number of expected bins and fDect1[i] fDect2[i]
2260   //
2261
2262   gStyle->SetPalette(1);
2263   gStyle->SetOptStat(1111);
2264   gStyle->SetPadBorderMode(0);
2265   gStyle->SetCanvasColor(10);
2266   gStyle->SetPadLeftMargin(0.13);
2267   gStyle->SetPadRightMargin(0.01);
2268
2269   // Mode groups of pads: the total number of bins!
2270   CalculNumberOfBinsExpected(i);
2271
2272   // Quick verification that we have the good pad calibration mode!
2273   if (fNumberOfBinsExpected != nbins) {
2274     AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2275     return kFALSE;
2276   }
2277
2278   // Security for fDebug 3 and 4
2279   if ((fDebugLevel >= 3) &&
2280       ((fDet[0] >  5) ||
2281        (fDet[1] >  4) ||
2282        (fDet[2] > 17))) {
2283     AliInfo("This detector doesn't exit!");
2284     return kFALSE;
2285   }
2286
2287   // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2288   CalculDect1Dect2(i);
2289
2290
2291   return kTRUE;
2292 }
2293 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2294 Bool_t AliTRDCalibraFit::InitFitCH()
2295 {
2296   //
2297   // Init the fVectorFitCH for normalisation
2298   // Init the histo for debugging
2299   //
2300
2301   gDirectory = gROOT;
2302
2303   fScaleFitFactor = 0.0;
2304   fCurrentCoefDetector   = new Float_t[2304];
2305   for (Int_t k = 0; k < 2304; k++) {
2306     fCurrentCoefDetector[k] = 0.0;
2307   }
2308   fVectorFit.SetName("gainfactorscoefficients");
2309
2310   // fDebug == 0 nothing
2311   // fDebug == 1 and fFitVoir no histo
2312   if (fDebugLevel == 1) {
2313     if(!CheckFitVoir()) return kFALSE;
2314   }
2315   //Get the CalDet object
2316   if(fAccCDB){
2317     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2318     if (!cal) {
2319       AliInfo("Could not get calibDB");
2320       return kFALSE;
2321     }
2322     if(fCalDet) delete fCalDet;
2323     fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2324   }
2325   else{
2326     Float_t devalue = 1.0;
2327     if(fCalDet) delete fCalDet;
2328     fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2329     for(Int_t k = 0; k < 540; k++){
2330       fCalDet->SetValue(k,devalue);
2331     }
2332   }
2333   return kTRUE;
2334
2335 }
2336 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2337 Bool_t AliTRDCalibraFit::InitFitPH()
2338 {
2339   //
2340   // Init the arrays of results
2341   // Init the histos for debugging
2342   //
2343
2344   gDirectory = gROOT;
2345   fVectorFit.SetName("driftvelocitycoefficients");
2346   fVectorFit2.SetName("t0coefficients");
2347
2348   fCurrentCoefDetector   = new Float_t[2304];
2349   for (Int_t k = 0; k < 2304; k++) {
2350     fCurrentCoefDetector[k] = 0.0;
2351   }
2352
2353   fCurrentCoefDetector2   = new Float_t[2304];
2354   for (Int_t k = 0; k < 2304; k++) {
2355     fCurrentCoefDetector2[k] = 0.0;
2356   }
2357
2358   //fDebug == 0 nothing
2359   // fDebug == 1 and fFitVoir no histo
2360   if (fDebugLevel == 1) {
2361     if(!CheckFitVoir()) return kFALSE;
2362   }
2363   //Get the CalDet object
2364   if(fAccCDB){
2365     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2366     if (!cal) {
2367       AliInfo("Could not get calibDB");
2368       return kFALSE;
2369     }
2370     if(fCalDet) delete fCalDet;
2371     if(fCalDet2) delete fCalDet2;
2372     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2373     fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2374   }
2375   else{
2376     Float_t devalue  = 1.5;
2377     Float_t devalue2 = 0.0;
2378     if(fCalDet) delete fCalDet;
2379     if(fCalDet2) delete fCalDet2;
2380     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2381     fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2382     for(Int_t k = 0; k < 540; k++){
2383       fCalDet->SetValue(k,devalue);
2384       fCalDet2->SetValue(k,devalue2);
2385     }
2386   }
2387   return kTRUE;
2388 }
2389 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2390 Bool_t AliTRDCalibraFit::InitFitPRF()
2391 {
2392   //
2393   // Init the calibration mode (Nz, Nrphi), the histograms for
2394   // debugging the fit methods if fDebug > 0,
2395   //
2396
2397   gDirectory = gROOT;
2398   fVectorFit.SetName("prfwidthcoefficients");
2399
2400   fCurrentCoefDetector   = new Float_t[2304];
2401   for (Int_t k = 0; k < 2304; k++) {
2402     fCurrentCoefDetector[k] = 0.0;
2403   }
2404
2405   // fDebug == 0 nothing
2406   // fDebug == 1 and fFitVoir no histo
2407   if (fDebugLevel == 1) {
2408     if(!CheckFitVoir()) return kFALSE;
2409   }
2410   return kTRUE;
2411 }
2412 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2413 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2414 {
2415   //
2416   // Init the fCalDet, fVectorFit fCurrentCoefDetector
2417   //
2418
2419   gDirectory = gROOT;
2420
2421   fCurrentCoefDetector   = new Float_t[2304];
2422   fCurrentCoefDetector2  = new Float_t[2304];
2423   for (Int_t k = 0; k < 2304; k++) {
2424     fCurrentCoefDetector[k]  = 0.0;
2425     fCurrentCoefDetector2[k] = 0.0;
2426   }
2427
2428   //printf("test0\n");
2429
2430   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2431   if (!cal) {
2432     AliInfo("Could not get calibDB");
2433     return kFALSE;
2434   }
2435
2436   //Get the CalDet object
2437   if(fAccCDB){
2438     if(fCalDet) delete fCalDet;
2439     if(fCalDet2) delete fCalDet2;
2440     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2441     //printf("test1\n");
2442     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2443     //printf("test2\n");
2444     for(Int_t k = 0; k < 540; k++){
2445       fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2446     }
2447     //printf("test3\n");
2448   }
2449   else{
2450     Float_t devalue  = 1.5;
2451     Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2452     if(fCalDet) delete fCalDet;
2453     if(fCalDet2) delete fCalDet2;
2454     //printf("test1\n");
2455     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2456     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2457     //printf("test2\n");
2458     for(Int_t k = 0; k < 540; k++){
2459       fCalDet->SetValue(k,devalue);
2460       fCalDet2->SetValue(k,devalue2);
2461     }
2462     //printf("test3\n");
2463   }
2464   return kTRUE;
2465 }
2466
2467 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2468 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2469 {
2470   //
2471   // Init the current detector where we are fCountDet and the
2472   // next fCount for the functions Fit...
2473   //
2474
2475   // Loop on the Xbins of ch!!
2476   fCountDet = -1; // Current detector
2477   fCount    =  0; // To find the next detector
2478
2479   // If fDebug >= 3
2480   if (fDebugLevel >= 3) {
2481     // Set countdet to the detector
2482     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2483     // Set counter to write at the end of the detector
2484     fCount = fDect2;
2485     // Get the right calib objects
2486     SetCalROC(i);
2487   }
2488   if(fDebugLevel == 1) {
2489     fCountDet = 0;
2490     fCalibraMode->CalculXBins(fCountDet,i);
2491     while(fCalibraMode->GetXbins(i) <=fFitVoir){
2492       fCountDet++;
2493       fCalibraMode->CalculXBins(fCountDet,i);
2494     }
2495     fCount    = fCalibraMode->GetXbins(i);
2496     fCountDet--;
2497     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2498     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2499     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2500                                       ,(Int_t) GetStack(fCountDet)
2501                                       ,(Int_t) GetSector(fCountDet),i);
2502   }
2503 }
2504 //_______________________________________________________________________________
2505 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2506 {
2507   //
2508   // Calculate the number of bins expected (calibration groups)
2509   //
2510
2511   fNumberOfBinsExpected = 0;
2512   // All
2513   if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2514     fNumberOfBinsExpected = 1;
2515     return;
2516   }
2517   // Per supermodule
2518   if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2519     fNumberOfBinsExpected = 18;
2520     return;
2521   }
2522   // More
2523   fCalibraMode->ModePadCalibration(2,i);
2524   fCalibraMode->ModePadFragmentation(0,2,0,i);
2525   fCalibraMode->SetDetChamb2(i);
2526   if (fDebugLevel > 1) {
2527     AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2528   }
2529   fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2530   fCalibraMode->ModePadCalibration(0,i);
2531   fCalibraMode->ModePadFragmentation(0,0,0,i);
2532   fCalibraMode->SetDetChamb0(i);
2533   if (fDebugLevel > 1) {
2534     AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2535   }
2536   fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2537
2538 }
2539 //_______________________________________________________________________________
2540 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2541 {
2542   //
2543   // Calculate the range of fits
2544   //
2545
2546   fDect1 = -1;
2547   fDect2 = -1;
2548   if (fDebugLevel == 1) {
2549     fDect1 = fFitVoir;
2550     fDect2 = fDect1 +1;
2551   }
2552   if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2553     fDect1 = 0;
2554     fDect2 = fNumberOfBinsExpected;
2555   }
2556   if (fDebugLevel >= 3) {
2557     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2558     fCalibraMode->CalculXBins(fCountDet,i);
2559     fDect1 = fCalibraMode->GetXbins(i);
2560     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2561     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2562     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2563                                       ,(Int_t) GetStack(fCountDet)
2564                                       ,(Int_t) GetSector(fCountDet),i);
2565     // Set for the next detector
2566     fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2567   }
2568 }
2569 //_______________________________________________________________________________
2570 Bool_t AliTRDCalibraFit::CheckFitVoir()
2571 {
2572   //
2573   // Check if fFitVoir is in the range
2574   //
2575
2576   if (fFitVoir < fNumberOfBinsExpected) {
2577     AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2578   }
2579   else {
2580     AliInfo("fFitVoir is out of range of the histo!");
2581     return kFALSE;
2582   }
2583   return kTRUE;
2584 }
2585 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2586 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2587 {
2588   //
2589   // See if we are in a new detector and update the
2590   // variables fNfragZ and fNfragRphi if yes
2591   // Will never happen for only one detector (3 and 4)
2592   // Doesn't matter for 2
2593   //
2594   if (fCount == idect) {
2595     // On en est au detector (or first detector in the group)
2596     fCountDet += 1;
2597     AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2598     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2599     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2600     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2601                                        ,(Int_t) GetStack(fCountDet)
2602                                        ,(Int_t) GetSector(fCountDet),i);
2603     // Set for the next detector
2604     fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2605     // calib objects
2606     SetCalROC(i);
2607   }
2608 }
2609 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2610 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2611 {
2612   //
2613   // Reconstruct the min pad row, max pad row, min pad col and
2614   // max pad col of the calibration group for the Fit functions
2615   // idect is the calibration group inside the detector
2616   //
2617   if (fDebugLevel !=  1) {
2618     fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2619   }
2620   AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2621   AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2622 }
2623 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2624 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2625 {
2626   //
2627   // For the case where there are not enough entries in the histograms
2628   // of the calibration group, the value present in the choosen database
2629   // will be put. A negativ sign enables to know that a fit was not possible.
2630   //
2631
2632   if (fDebugLevel == 1) {
2633     AliInfo("The element has not enough statistic to be fitted");
2634   }
2635   else if (fNbDet > 0){
2636     Int_t firstdetector = fCountDet;
2637     Int_t lastdetector  = fCountDet+fNbDet;
2638     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2639                  ,idect,firstdetector,lastdetector));
2640     // loop over detectors
2641     for(Int_t det = firstdetector; det < lastdetector; det++){
2642
2643       //Set the calibration object again
2644       fCountDet = det;
2645       SetCalROC(0);
2646
2647       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2648       // Put them at 1
2649       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2650       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2651                                          ,(Int_t) GetStack(fCountDet)
2652                                          ,(Int_t) GetSector(fCountDet),0);
2653       // Reconstruct row min row max
2654       ReconstructFitRowMinRowMax(idect,0);
2655
2656       // Calcul the coef from the database choosen for the detector
2657       CalculChargeCoefMean(kFALSE);
2658
2659       //stack 2, not stack 2
2660       Int_t factor = 0;
2661       if(GetStack(fCountDet) == 2) factor = 12;
2662       else factor = 16;
2663
2664       // Fill the fCurrentCoefDetector with negative value to say: not fitted
2665       for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2666         for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2667           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2668         }
2669       }
2670
2671       //Put default value negative
2672       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2673       fCurrentCoefE   = 0.0;
2674
2675       // Fill the stuff
2676       FillVectorFit();
2677       // Debug
2678       if(fDebugLevel > 1){
2679
2680         if ( !fDebugStreamer ) {
2681           //debug stream
2682           TDirectory *backup = gDirectory;
2683           fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2684           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2685         }
2686
2687         Int_t   detector   = fCountDet;
2688         Int_t   caligroup  = idect;
2689         Short_t rowmin     = fCalibraMode->GetRowMin(0);
2690         Short_t rowmax     = fCalibraMode->GetRowMax(0);
2691         Short_t colmin     = fCalibraMode->GetColMin(0);
2692         Short_t colmax     = fCalibraMode->GetColMax(0);
2693         Float_t gf         = fCurrentCoef[0];
2694         Float_t gfs        = fCurrentCoef[1];
2695         Float_t gfE        = fCurrentCoefE;
2696
2697         (*fDebugStreamer) << "FillFillCH" <<
2698           "detector=" << detector <<
2699           "caligroup=" << caligroup <<
2700           "rowmin=" << rowmin <<
2701           "rowmax=" << rowmax <<
2702           "colmin=" << colmin <<
2703           "colmax=" << colmax <<
2704           "gf=" << gf <<
2705           "gfs=" << gfs <<
2706           "gfE=" << gfE <<
2707           "\n";
2708
2709       }
2710       // Reset
2711       for (Int_t k = 0; k < 2304; k++) {
2712         fCurrentCoefDetector[k] = 0.0;
2713       }
2714
2715     }// loop detector
2716     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2717   }
2718   else {
2719
2720     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2721                  ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2722
2723     // Calcul the coef from the database choosen
2724     CalculChargeCoefMean(kFALSE);
2725
2726     //stack 2, not stack 2
2727     Int_t factor = 0;
2728     if(GetStack(fCountDet) == 2) factor = 12;
2729     else factor = 16;
2730
2731     // Fill the fCurrentCoefDetector with negative value to say: not fitted
2732     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2733       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2734         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2735       }
2736     }
2737
2738     //Put default value negative
2739     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2740     fCurrentCoefE   = 0.0;
2741
2742     FillFillCH(idect);
2743   }
2744
2745   return kTRUE;
2746 }
2747
2748
2749 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2750 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2751 {
2752   //
2753   // For the case where there are not enough entries in the histograms
2754   // of the calibration group, the value present in the choosen database
2755   // will be put. A negativ sign enables to know that a fit was not possible.
2756   //
2757   if (fDebugLevel == 1) {
2758     AliInfo("The element has not enough statistic to be fitted");
2759   }
2760   else if (fNbDet > 0) {
2761
2762     Int_t firstdetector = fCountDet;
2763     Int_t lastdetector  = fCountDet+fNbDet;
2764     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2765                  ,idect,firstdetector,lastdetector));
2766     // loop over detectors
2767     for(Int_t det = firstdetector; det < lastdetector; det++){
2768
2769       //Set the calibration object again
2770       fCountDet = det;
2771       SetCalROC(1);
2772
2773       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2774       // Put them at 1
2775       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
2776       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2777                                          ,(Int_t) GetStack(fCountDet)
2778                                          ,(Int_t) GetSector(fCountDet),1);
2779       // Reconstruct row min row max
2780       ReconstructFitRowMinRowMax(idect,1);
2781
2782       // Calcul the coef from the database choosen for the detector
2783       CalculVdriftCoefMean();
2784       CalculT0CoefMean();
2785
2786       //stack 2, not stack 2
2787       Int_t factor = 0;
2788       if(GetStack(fCountDet) == 2) factor = 12;
2789       else factor = 16;
2790
2791       // Fill the fCurrentCoefDetector with negative value to say: not fitted
2792       for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2793         for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2794           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2795           fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2796         }
2797       }
2798
2799       //Put default value negative
2800       fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2801       fCurrentCoefE    = 0.0;
2802       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2803       fCurrentCoefE2   = 0.0;
2804
2805       // Fill the stuff
2806       FillVectorFit();
2807       FillVectorFit2();
2808       // Debug
2809       if(fDebugLevel > 1){
2810
2811         if ( !fDebugStreamer ) {
2812           //debug stream
2813           TDirectory *backup = gDirectory;
2814           fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
2815           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2816         }
2817
2818
2819         Int_t   detector     = fCountDet;
2820         Int_t   caligroup    = idect;
2821         Short_t rowmin       = fCalibraMode->GetRowMin(1);
2822         Short_t rowmax       = fCalibraMode->GetRowMax(1);
2823         Short_t colmin       = fCalibraMode->GetColMin(1);
2824         Short_t colmax       = fCalibraMode->GetColMax(1);
2825         Float_t vf           = fCurrentCoef[0];
2826         Float_t vs           = fCurrentCoef[1];
2827         Float_t vfE          = fCurrentCoefE;
2828         Float_t t0f          = fCurrentCoef2[0];
2829         Float_t t0s          = fCurrentCoef2[1];
2830         Float_t t0E          = fCurrentCoefE2;
2831
2832
2833
2834         (* fDebugStreamer) << "FillFillPH"<<
2835         "detector="<<detector<<
2836           "nentries="<<nentries<<
2837           "caligroup="<<caligroup<<
2838           "rowmin="<<rowmin<<
2839           "rowmax="<<rowmax<<
2840           "colmin="<<colmin<<
2841           "colmax="<<colmax<<
2842           "vf="<<vf<<
2843           "vs="<<vs<<
2844           "vfE="<<vfE<<
2845           "t0f="<<t0f<<
2846           "t0s="<<t0s<<
2847           "t0E="<<t0E<<
2848           "\n";
2849       }
2850       // Reset
2851       for (Int_t k = 0; k < 2304; k++) {
2852         fCurrentCoefDetector[k] = 0.0;
2853         fCurrentCoefDetector2[k] = 0.0;
2854       }
2855
2856     }// loop detector
2857     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2858   }
2859   else {
2860
2861     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2862                  ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
2863
2864     CalculVdriftCoefMean();
2865     CalculT0CoefMean();
2866
2867     //stack 2 and not stack 2
2868     Int_t factor = 0;
2869     if(GetStack(fCountDet) == 2) factor = 12;
2870     else factor = 16;
2871
2872
2873     // Fill the fCurrentCoefDetector 2
2874     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
2875       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
2876         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2877         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
2878       }
2879     }
2880
2881     // Put the default value
2882     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
2883     fCurrentCoefE    = 0.0;
2884     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
2885     fCurrentCoefE2   = 0.0;
2886
2887     FillFillPH(idect,nentries);
2888
2889   }
2890
2891   return kTRUE;
2892
2893 }
2894
2895
2896 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2897 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
2898 {
2899   //
2900   // For the case where there are not enough entries in the histograms
2901   // of the calibration group, the value present in the choosen database
2902   // will be put. A negativ sign enables to know that a fit was not possible.
2903   //
2904
2905   if (fDebugLevel == 1) {
2906     AliInfo("The element has not enough statistic to be fitted");
2907   }
2908   else if (fNbDet > 0){
2909
2910     Int_t firstdetector = fCountDet;
2911     Int_t lastdetector  = fCountDet+fNbDet;
2912     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2913                  ,idect,firstdetector,lastdetector));
2914
2915     // loop over detectors
2916     for(Int_t det = firstdetector; det < lastdetector; det++){
2917
2918       //Set the calibration object again
2919       fCountDet = det;
2920       SetCalROC(2);
2921
2922       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2923       // Put them at 1
2924       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
2925       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2926                                          ,(Int_t) GetStack(fCountDet)
2927                                          ,(Int_t) GetSector(fCountDet),2);
2928       // Reconstruct row min row max
2929       ReconstructFitRowMinRowMax(idect,2);
2930
2931       // Calcul the coef from the database choosen for the detector
2932       CalculPRFCoefMean();
2933
2934       //stack 2, not stack 2
2935       Int_t factor = 0;
2936       if(GetStack(fCountDet) == 2) factor = 12;
2937       else factor = 16;
2938
2939       // Fill the fCurrentCoefDetector with negative value to say: not fitted
2940       for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
2941         for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
2942           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2943         }
2944       }
2945
2946       //Put default value negative
2947       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2948       fCurrentCoefE   = 0.0;
2949
2950       // Fill the stuff
2951       FillVectorFit();
2952       // Debug
2953       if(fDebugLevel > 1){
2954
2955         if ( !fDebugStreamer ) {
2956           //debug stream
2957           TDirectory *backup = gDirectory;
2958           fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
2959           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2960         }
2961
2962         Int_t   detector     = fCountDet;
2963         Int_t   layer        = GetLayer(fCountDet);
2964         Int_t   caligroup    = idect;
2965         Short_t rowmin       = fCalibraMode->GetRowMin(2);
2966         Short_t rowmax       = fCalibraMode->GetRowMax(2);
2967         Short_t colmin       = fCalibraMode->GetColMin(2);
2968         Short_t colmax       = fCalibraMode->GetColMax(2);
2969         Float_t widf         = fCurrentCoef[0];
2970         Float_t wids         = fCurrentCoef[1];
2971         Float_t widfE        = fCurrentCoefE;
2972
2973         (* fDebugStreamer) << "FillFillPRF"<<
2974           "detector="<<detector<<
2975           "layer="<<layer<<
2976           "caligroup="<<caligroup<<
2977           "rowmin="<<rowmin<<
2978           "rowmax="<<rowmax<<
2979           "colmin="<<colmin<<
2980           "colmax="<<colmax<<
2981           "widf="<<widf<<
2982           "wids="<<wids<<
2983           "widfE="<<widfE<<
2984           "\n";
2985       }
2986       // Reset
2987       for (Int_t k = 0; k < 2304; k++) {
2988         fCurrentCoefDetector[k] = 0.0;
2989       }
2990
2991     }// loop detector
2992     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2993   }
2994   else {
2995
2996     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2997                  ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
2998
2999     CalculPRFCoefMean();
3000
3001     // stack 2 and not stack 2
3002     Int_t factor = 0;
3003     if(GetStack(fCountDet) == 2) factor = 12;
3004     else factor = 16;
3005
3006
3007     // Fill the fCurrentCoefDetector
3008     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3009       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3010         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3011       }
3012     }
3013
3014     // Put the default value
3015     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3016     fCurrentCoefE   = 0.0;
3017
3018     FillFillPRF(idect);
3019   }
3020
3021   return kTRUE;
3022
3023 }
3024 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3025 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3026 {
3027   //
3028   // For the case where there are not enough entries in the histograms
3029   // of the calibration group, the value present in the choosen database
3030   // will be put. A negativ sign enables to know that a fit was not possible.
3031   //
3032
3033   // Calcul the coef from the database choosen
3034   CalculVdriftLorentzCoef();
3035
3036   Int_t factor = 0;
3037   if(GetStack(fCountDet) == 2) factor = 1728;
3038   else factor = 2304;
3039
3040
3041   // Fill the fCurrentCoefDetector
3042   for (Int_t k = 0; k < factor; k++) {
3043     fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3044     // should be negative
3045     fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3046   }
3047
3048
3049   //Put default opposite sign
3050   fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3051   fCurrentCoefE    = 0.0;
3052   fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3053   fCurrentCoefE2 = 0.0;
3054
3055   FillFillLinearFitter();
3056
3057   return kTRUE;
3058 }
3059
3060 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3061 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3062 {
3063   //
3064   // Fill the coefficients found with the fits or other
3065   // methods from the Fit functions
3066   //
3067
3068   if (fDebugLevel != 1) {
3069     if (fNbDet > 0){
3070       Int_t firstdetector = fCountDet;
3071       Int_t lastdetector  = fCountDet+fNbDet;
3072       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3073                    ,idect,firstdetector,lastdetector));
3074       // loop over detectors
3075       for(Int_t det = firstdetector; det < lastdetector; det++){
3076
3077         //Set the calibration object again
3078         fCountDet = det;
3079         SetCalROC(0);
3080
3081         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3082         // Put them at 1
3083         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3084         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3085                                            ,(Int_t) GetStack(fCountDet)
3086                                            ,(Int_t) GetSector(fCountDet),0);
3087         // Reconstruct row min row max
3088         ReconstructFitRowMinRowMax(idect,0);
3089
3090         // Calcul the coef from the database choosen for the detector
3091         if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3092         else CalculChargeCoefMean(kTRUE);
3093
3094         //stack 2, not stack 2
3095         Int_t factor = 0;
3096         if(GetStack(fCountDet) == 2) factor = 12;
3097         else factor = 16;
3098
3099         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3100         Double_t coeftoput = 1.0;
3101         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3102         else coeftoput = fCurrentCoef[0];
3103         for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3104           for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3105             fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3106           }
3107         }
3108
3109         // Fill the stuff
3110         FillVectorFit();
3111         // Debug
3112         if(fDebugLevel > 1){
3113
3114           if ( !fDebugStreamer ) {
3115             //debug stream
3116             TDirectory *backup = gDirectory;
3117             fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3118             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3119           }
3120
3121           Int_t   detector   = fCountDet;
3122           Int_t   caligroup  = idect;
3123           Short_t rowmin     = fCalibraMode->GetRowMin(0);
3124           Short_t rowmax     = fCalibraMode->GetRowMax(0);
3125           Short_t colmin     = fCalibraMode->GetColMin(0);
3126           Short_t colmax     = fCalibraMode->GetColMax(0);
3127           Float_t gf         = fCurrentCoef[0];
3128           Float_t gfs        = fCurrentCoef[1];
3129           Float_t gfE        = fCurrentCoefE;
3130
3131           (*fDebugStreamer) << "FillFillCH" <<
3132             "detector=" << detector <<
3133             "caligroup=" << caligroup <<
3134             "rowmin=" << rowmin <<
3135             "rowmax=" << rowmax <<
3136             "colmin=" << colmin <<
3137             "colmax=" << colmax <<
3138             "gf=" << gf <<
3139             "gfs=" << gfs <<
3140             "gfE=" << gfE <<
3141             "\n";
3142
3143         }
3144         // Reset
3145         for (Int_t k = 0; k < 2304; k++) {
3146           fCurrentCoefDetector[k] = 0.0;
3147         }
3148
3149       }// loop detector
3150       //printf("Check the count now: fCountDet %d\n",fCountDet);
3151     }
3152     else{
3153
3154       Int_t factor = 0;
3155       if(GetStack(fCountDet) == 2) factor = 12;
3156       else factor = 16;
3157
3158       for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3159         for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3160           fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3161         }
3162       }
3163
3164       FillFillCH(idect);
3165     }
3166   }
3167
3168   return kTRUE;
3169
3170 }
3171 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3172 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3173 {
3174   //
3175   // Fill the coefficients found with the fits or other
3176   // methods from the Fit functions
3177   //
3178
3179   if (fDebugLevel != 1) {
3180     if (fNbDet > 0){
3181
3182       Int_t firstdetector = fCountDet;
3183       Int_t lastdetector  = fCountDet+fNbDet;
3184       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3185                    ,idect,firstdetector,lastdetector));
3186
3187       // loop over detectors
3188       for(Int_t det = firstdetector; det < lastdetector; det++){
3189
3190         //Set the calibration object again
3191         fCountDet = det;
3192         SetCalROC(1);
3193
3194         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3195         // Put them at 1
3196         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3197         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3198                                            ,(Int_t) GetStack(fCountDet)
3199                                            ,(Int_t) GetSector(fCountDet),1);
3200         // Reconstruct row min row max
3201         ReconstructFitRowMinRowMax(idect,1);
3202
3203         // Calcul the coef from the database choosen for the detector
3204         CalculVdriftCoefMean();
3205         CalculT0CoefMean();
3206
3207         //stack 2, not stack 2
3208         Int_t factor = 0;
3209         if(GetStack(fCountDet) == 2) factor = 12;
3210         else factor = 16;
3211
3212         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3213         Double_t coeftoput  = 1.5;
3214         Double_t coeftoput2 = 0.0;
3215
3216         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3217         else coeftoput = fCurrentCoef[0];
3218
3219         if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3220         else coeftoput2 = fCurrentCoef2[0];
3221
3222         for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3223           for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3224             fCurrentCoefDetector[(Int_t)(j*factor+k)]  = coeftoput;
3225             fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3226           }
3227         }
3228
3229         // Fill the stuff
3230         FillVectorFit();
3231         FillVectorFit2();
3232         // Debug
3233         if(fDebugLevel > 1){
3234
3235           if ( !fDebugStreamer ) {
3236             //debug stream
3237             TDirectory *backup = gDirectory;
3238             fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3239             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3240           }
3241
3242
3243           Int_t   detector     = fCountDet;
3244           Int_t   caligroup    = idect;
3245           Short_t rowmin       = fCalibraMode->GetRowMin(1);
3246           Short_t rowmax       = fCalibraMode->GetRowMax(1);
3247           Short_t colmin       = fCalibraMode->GetColMin(1);
3248           Short_t colmax       = fCalibraMode->GetColMax(1);
3249           Float_t vf           = fCurrentCoef[0];
3250           Float_t vs           = fCurrentCoef[1];
3251           Float_t vfE          = fCurrentCoefE;
3252           Float_t t0f          = fCurrentCoef2[0];
3253           Float_t t0s          = fCurrentCoef2[1];
3254           Float_t t0E          = fCurrentCoefE2;
3255
3256
3257
3258           (* fDebugStreamer) << "FillFillPH"<<
3259             "detector="<<detector<<
3260             "nentries="<<nentries<<
3261             "caligroup="<<caligroup<<
3262             "rowmin="<<rowmin<<
3263             "rowmax="<<rowmax<<
3264             "colmin="<<colmin<<
3265             "colmax="<<colmax<<
3266             "vf="<<vf<<
3267             "vs="<<vs<<
3268             "vfE="<<vfE<<
3269             "t0f="<<t0f<<
3270             "t0s="<<t0s<<
3271             "t0E="<<t0E<<
3272             "\n";
3273         }
3274         // Reset
3275         for (Int_t k = 0; k < 2304; k++) {
3276           fCurrentCoefDetector[k] = 0.0;
3277           fCurrentCoefDetector2[k] = 0.0;
3278         }
3279
3280       }// loop detector
3281       //printf("Check the count now: fCountDet %d\n",fCountDet);
3282     }
3283     else {
3284
3285       Int_t factor = 0;
3286       if(GetStack(fCountDet) == 2) factor = 12;
3287       else factor = 16;
3288
3289       for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3290         for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3291           fCurrentCoefDetector[(Int_t)(j*factor+k)]  = fCurrentCoef[0];
3292           fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3293         }
3294       }
3295
3296       FillFillPH(idect,nentries);
3297     }
3298   }
3299   return kTRUE;
3300 }
3301 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3302 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3303 {
3304   //
3305   // Fill the coefficients found with the fits or other
3306   // methods from the Fit functions
3307   //
3308
3309   if (fDebugLevel != 1) {
3310     if (fNbDet > 0){
3311
3312       Int_t firstdetector = fCountDet;
3313       Int_t lastdetector  = fCountDet+fNbDet;
3314       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3315                    ,idect,firstdetector,lastdetector));
3316
3317       // loop over detectors
3318       for(Int_t det = firstdetector; det < lastdetector; det++){
3319
3320         //Set the calibration object again
3321         fCountDet = det;
3322         SetCalROC(2);
3323
3324         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3325         // Put them at 1
3326         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3327         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3328                                            ,(Int_t) GetStack(fCountDet)
3329                                            ,(Int_t) GetSector(fCountDet),2);
3330         // Reconstruct row min row max
3331         ReconstructFitRowMinRowMax(idect,2);
3332
3333         // Calcul the coef from the database choosen for the detector
3334         CalculPRFCoefMean();
3335
3336         //stack 2, not stack 2
3337         Int_t factor = 0;
3338         if(GetStack(fCountDet) == 2) factor = 12;
3339         else factor = 16;
3340
3341         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3342         Double_t coeftoput = 1.0;
3343         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3344         else coeftoput = fCurrentCoef[0];
3345         for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3346           for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3347             fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3348           }
3349         }
3350
3351         // Fill the stuff
3352         FillVectorFit();
3353         // Debug
3354         if(fDebugLevel > 1){
3355
3356           if ( !fDebugStreamer ) {
3357             //debug stream
3358             TDirectory *backup = gDirectory;
3359             fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3360             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3361           }
3362
3363           Int_t   detector     = fCountDet;
3364           Int_t   layer        = GetLayer(fCountDet);
3365           Int_t   caligroup    = idect;
3366           Short_t rowmin       = fCalibraMode->GetRowMin(2);
3367           Short_t rowmax       = fCalibraMode->GetRowMax(2);
3368           Short_t colmin       = fCalibraMode->GetColMin(2);
3369           Short_t colmax       = fCalibraMode->GetColMax(2);
3370           Float_t widf         = fCurrentCoef[0];
3371           Float_t wids         = fCurrentCoef[1];
3372           Float_t widfE        = fCurrentCoefE;
3373
3374           (* fDebugStreamer) << "FillFillPRF"<<
3375             "detector="<<detector<<
3376             "layer="<<layer<<
3377             "caligroup="<<caligroup<<
3378             "rowmin="<<rowmin<<
3379             "rowmax="<<rowmax<<
3380             "colmin="<<colmin<<
3381             "colmax="<<colmax<<
3382             "widf="<<widf<<
3383             "wids="<<wids<<
3384             "widfE="<<widfE<<
3385             "\n";
3386         }
3387         // Reset
3388         for (Int_t k = 0; k < 2304; k++) {
3389           fCurrentCoefDetector[k] = 0.0;
3390         }
3391
3392       }// loop detector
3393       //printf("Check the count now: fCountDet %d\n",fCountDet);
3394     }
3395     else {
3396
3397       Int_t factor = 0;
3398       if(GetStack(fCountDet) == 2) factor = 12;
3399       else factor = 16;
3400
3401       // Pointer to the branch
3402       for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3403         for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3404           fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3405         }
3406       }
3407       FillFillPRF(idect);
3408     }
3409   }
3410
3411   return kTRUE;
3412
3413 }
3414 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3415 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3416 {
3417   //
3418   // Fill the coefficients found with the fits or other
3419   // methods from the Fit functions
3420   //
3421
3422   Int_t factor = 0;
3423   if(GetStack(fCountDet) == 2) factor = 1728;
3424   else factor = 2304;
3425
3426   // Pointer to the branch
3427   for (Int_t k = 0; k < factor; k++) {
3428     fCurrentCoefDetector[k]  = fCurrentCoef[0];
3429     fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3430   }
3431
3432   FillFillLinearFitter();
3433
3434   return kTRUE;
3435
3436 }
3437 //________________________________________________________________________________
3438 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3439 {
3440   //
3441   // DebugStream and fVectorFit
3442   //
3443
3444   // End of one detector
3445   if ((idect == (fCount-1))) {
3446     FillVectorFit();
3447     // Reset
3448     for (Int_t k = 0; k < 2304; k++) {
3449       fCurrentCoefDetector[k] = 0.0;
3450     }
3451   }
3452
3453   if(fDebugLevel > 1){
3454
3455     if ( !fDebugStreamer ) {
3456       //debug stream
3457       TDirectory *backup = gDirectory;
3458       fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3459       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3460     }
3461
3462     Int_t   detector   = fCountDet;
3463     Int_t   caligroup  = idect;
3464     Short_t rowmin     = fCalibraMode->GetRowMin(0);
3465     Short_t rowmax     = fCalibraMode->GetRowMax(0);
3466     Short_t colmin     = fCalibraMode->GetColMin(0);
3467     Short_t colmax     = fCalibraMode->GetColMax(0);
3468     Float_t gf         = fCurrentCoef[0];
3469     Float_t gfs        = fCurrentCoef[1];
3470     Float_t gfE        = fCurrentCoefE;
3471
3472     (*fDebugStreamer) << "FillFillCH" <<
3473       "detector=" << detector <<
3474       "caligroup=" << caligroup <<
3475       "rowmin=" << rowmin <<
3476       "rowmax=" << rowmax <<
3477       "colmin=" << colmin <<
3478       "colmax=" << colmax <<
3479       "gf=" << gf <<
3480       "gfs=" << gfs <<
3481       "gfE=" << gfE <<
3482       "\n";
3483
3484   }
3485 }
3486 //________________________________________________________________________________
3487 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3488 {
3489   //
3490   // DebugStream and fVectorFit and fVectorFit2
3491   //
3492
3493   // End of one detector
3494     if ((idect == (fCount-1))) {
3495       FillVectorFit();
3496       FillVectorFit2();
3497       // Reset
3498       for (Int_t k = 0; k < 2304; k++) {
3499         fCurrentCoefDetector[k] = 0.0;
3500         fCurrentCoefDetector2[k] = 0.0;
3501       }
3502     }
3503
3504     if(fDebugLevel > 1){
3505
3506       if ( !fDebugStreamer ) {
3507         //debug stream
3508         TDirectory *backup = gDirectory;
3509         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3510         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3511       }
3512
3513
3514       Int_t   detector     = fCountDet;
3515       Int_t   caligroup    = idect;
3516       Short_t rowmin       = fCalibraMode->GetRowMin(1);
3517       Short_t rowmax       = fCalibraMode->GetRowMax(1);
3518       Short_t colmin       = fCalibraMode->GetColMin(1);
3519       Short_t colmax       = fCalibraMode->GetColMax(1);
3520       Float_t vf           = fCurrentCoef[0];
3521       Float_t vs           = fCurrentCoef[1];
3522       Float_t vfE          = fCurrentCoefE;
3523       Float_t t0f          = fCurrentCoef2[0];
3524       Float_t t0s          = fCurrentCoef2[1];
3525       Float_t t0E          = fCurrentCoefE2;
3526
3527
3528
3529       (* fDebugStreamer) << "FillFillPH"<<
3530         "detector="<<detector<<
3531         "nentries="<<nentries<<
3532         "caligroup="<<caligroup<<
3533         "rowmin="<<rowmin<<
3534         "rowmax="<<rowmax<<
3535         "colmin="<<colmin<<
3536         "colmax="<<colmax<<
3537         "vf="<<vf<<
3538         "vs="<<vs<<
3539         "vfE="<<vfE<<
3540         "t0f="<<t0f<<
3541         "t0s="<<t0s<<
3542         "t0E="<<t0E<<
3543         "\n";
3544     }
3545
3546 }
3547 //________________________________________________________________________________
3548 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3549 {
3550   //
3551   // DebugStream and fVectorFit
3552   //
3553
3554     // End of one detector
3555     if ((idect == (fCount-1))) {
3556       FillVectorFit();
3557       // Reset
3558       for (Int_t k = 0; k < 2304; k++) {
3559         fCurrentCoefDetector[k] = 0.0;
3560       }
3561     }
3562
3563
3564     if(fDebugLevel > 1){
3565
3566       if ( !fDebugStreamer ) {
3567         //debug stream
3568         TDirectory *backup = gDirectory;
3569         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3570         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3571       }
3572
3573       Int_t   detector     = fCountDet;
3574       Int_t   layer        = GetLayer(fCountDet);
3575       Int_t   caligroup    = idect;
3576       Short_t rowmin       = fCalibraMode->GetRowMin(2);
3577       Short_t rowmax       = fCalibraMode->GetRowMax(2);
3578       Short_t colmin       = fCalibraMode->GetColMin(2);
3579       Short_t colmax       = fCalibraMode->GetColMax(2);
3580       Float_t widf         = fCurrentCoef[0];
3581       Float_t wids         = fCurrentCoef[1];
3582       Float_t widfE        = fCurrentCoefE;
3583
3584       (* fDebugStreamer) << "FillFillPRF"<<
3585         "detector="<<detector<<
3586         "layer="<<layer<<
3587         "caligroup="<<caligroup<<
3588         "rowmin="<<rowmin<<
3589         "rowmax="<<rowmax<<
3590         "colmin="<<colmin<<
3591         "colmax="<<colmax<<
3592         "widf="<<widf<<
3593         "wids="<<wids<<
3594         "widfE="<<widfE<<
3595         "\n";
3596     }
3597
3598 }
3599 //________________________________________________________________________________
3600 void AliTRDCalibraFit::FillFillLinearFitter()
3601 {
3602   //
3603   // DebugStream and fVectorFit
3604   //
3605
3606   // End of one detector
3607   FillVectorFit();
3608   FillVectorFit2();
3609
3610
3611   // Reset
3612   for (Int_t k = 0; k < 2304; k++) {
3613   fCurrentCoefDetector[k]  = 0.0;
3614   fCurrentCoefDetector2[k] = 0.0;
3615   }
3616
3617
3618   if(fDebugLevel > 1){
3619
3620     if ( !fDebugStreamer ) {
3621       //debug stream
3622       TDirectory *backup = gDirectory;
3623       fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3624       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3625     }
3626
3627     //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3628     AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3629     Float_t rowmd            = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3630     Float_t r                = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3631     Float_t tiltangle        = padplane->GetTiltingAngle();
3632     Int_t   detector         = fCountDet;
3633     Int_t   stack            = GetStack(fCountDet);
3634     Int_t   layer            = GetLayer(fCountDet);
3635     Float_t vf               = fCurrentCoef[0];
3636     Float_t vs               = fCurrentCoef[1];
3637     Float_t vfE              = fCurrentCoefE;
3638     Float_t lorentzangler    = fCurrentCoef2[0];
3639     Float_t elorentzangler   = fCurrentCoefE2;
3640     Float_t lorentzangles    = fCurrentCoef2[1];
3641
3642     (* fDebugStreamer) << "FillFillLinearFitter"<<
3643       "detector="<<detector<<
3644       "stack="<<stack<<
3645       "layer="<<layer<<
3646       "rowmd="<<rowmd<<
3647       "r="<<r<<
3648       "tiltangle="<<tiltangle<<
3649       "vf="<<vf<<
3650       "vs="<<vs<<
3651       "vfE="<<vfE<<
3652       "lorentzangler="<<lorentzangler<<
3653       "Elorentzangler="<<elorentzangler<<
3654       "lorentzangles="<<lorentzangles<<
3655       "\n";
3656   }
3657
3658 }
3659 //
3660 //____________Calcul Coef Mean_________________________________________________
3661 //
3662 //_____________________________________________________________________________
3663 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3664 {
3665   //
3666   // For the detector Dect calcul the mean time 0
3667   // for the calibration group idect from the choosen database
3668   //
3669
3670   fCurrentCoef2[1] = 0.0;
3671   if(fDebugLevel != 1){
3672     if(((fCalibraMode->GetNz(1) > 0) ||
3673        (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1)    != 10) && (fCalibraMode->GetNz(1)    != 100))) {
3674
3675       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3676         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3677           fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3678         }
3679       }
3680
3681       fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3682
3683     }
3684     else {
3685
3686       if(!fAccCDB){
3687         fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3688       }
3689       else{
3690
3691         for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3692           for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3693             fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3694           }
3695         }
3696         fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3697
3698       }
3699     }
3700   }
3701   return kTRUE;
3702 }
3703
3704 //_____________________________________________________________________________
3705 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3706 {
3707   //
3708   // For the detector Dect calcul the mean gain factor
3709   // for the calibration group idect from the choosen database
3710   //
3711
3712   fCurrentCoef[1] = 0.0;
3713   if(fDebugLevel != 1){
3714     if (((fCalibraMode->GetNz(0)    > 0) ||
3715         (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0)    != 10) && (fCalibraMode->GetNz(0)    != 100))) {
3716       for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3717         for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3718           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3719           if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3720         }
3721       }
3722       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3723     }
3724     else {
3725       //Per detectors
3726       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3727       if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3728     }
3729   }
3730   return kTRUE;
3731 }
3732 //_____________________________________________________________________________
3733 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3734 {
3735   //
3736   // For the detector Dect calcul the mean sigma of pad response
3737   // function for the calibration group idect from the choosen database
3738   //
3739
3740   fCurrentCoef[1] = 0.0;
3741   if(fDebugLevel != 1){
3742     for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3743       for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3744         fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3745       }
3746     }
3747     fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3748   }
3749   return kTRUE;
3750 }
3751 //_____________________________________________________________________________
3752 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3753 {
3754   //
3755   // For the detector dect calcul the mean drift velocity for the
3756   // calibration group idect from the choosen database
3757   //
3758
3759   fCurrentCoef[1] = 0.0;
3760   if(fDebugLevel != 1){
3761     if (((fCalibraMode->GetNz(1)    > 0) ||
3762         (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1)    != 10) && (fCalibraMode->GetNz(1)    != 100))) {
3763
3764       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3765         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3766           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3767         }
3768       }
3769
3770       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3771
3772     }
3773     else {
3774       //per detectors
3775       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3776     }
3777   }
3778   return kTRUE;
3779 }
3780 //_____________________________________________________________________________
3781 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
3782 {
3783   //
3784   // For the detector fCountDet, mean drift velocity and tan lorentzangle
3785   //
3786
3787   fCurrentCoef[1]  = fCalDet->GetValue(fCountDet);
3788   fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3789
3790   return kTRUE;
3791 }
3792 //_____________________________________________________________________________
3793 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
3794 {
3795   //
3796   // Default width of the PRF if there is no database as reference
3797   //
3798   switch(layer)
3799     {
3800       // default database
3801       //case 0:  return 0.515;
3802       //case 1:  return 0.502;
3803       //case 2:  return 0.491;
3804       //case 3:  return 0.481;
3805       //case 4:  return 0.471;
3806       //case 5:  return 0.463;
3807       //default: return 0.0;
3808
3809       // fit database
3810     case 0:  return 0.538429;
3811     case 1:  return 0.524302;
3812     case 2:  return 0.511591;
3813     case 3:  return 0.500140;
3814     case 4:  return 0.489821;
3815     case 5:  return 0.480524;
3816     default: return 0.0;
3817   }
3818 }
3819 //________________________________________________________________________________
3820 void AliTRDCalibraFit::SetCalROC(Int_t i)
3821 {
3822   //
3823   // Set the calib object for fCountDet
3824   //
3825
3826   Float_t value = 0.0;
3827
3828   //Get the CalDet object
3829   if(fAccCDB){
3830     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
3831     if (!cal) {
3832       AliInfo("Could not get calibDB");
3833       return;
3834     }
3835     switch (i)
3836       {
3837       case 0:
3838         if( fCalROC ){
3839           fCalROC->~AliTRDCalROC();
3840           new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3841         }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
3842         break;
3843       case 1:
3844         if( fCalROC ){
3845           fCalROC->~AliTRDCalROC();
3846           new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3847         }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
3848         if( fCalROC2 ){
3849           fCalROC2->~AliTRDCalROC();
3850           new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3851         }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
3852         break;
3853       case 2:
3854         if( fCalROC ){
3855           fCalROC->~AliTRDCalROC();
3856           new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3857         }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
3858         break;
3859       default: return;
3860       }
3861   }
3862   else{
3863     switch (i)
3864       {
3865       case 0:
3866         if(fCalROC) delete fCalROC;
3867         fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3868         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3869           fCalROC->SetValue(k,1.0);
3870         }
3871         break;
3872       case 1:
3873         if(fCalROC)  delete fCalROC;
3874         if(fCalROC2) delete fCalROC2;
3875         fCalROC  = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3876         fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3877         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3878           fCalROC->SetValue(k,1.0);
3879           fCalROC2->SetValue(k,0.0);
3880         }
3881         break;
3882       case 2:
3883         if(fCalROC) delete fCalROC;
3884         value = GetPRFDefault(GetLayer(fCountDet));
3885         fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
3886         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
3887           fCalROC->SetValue(k,value);
3888         }
3889         break;
3890       default: return;
3891       }
3892   }
3893
3894 }
3895 //____________Fit Methods______________________________________________________
3896
3897 //_____________________________________________________________________________
3898 void AliTRDCalibraFit::FitPente(TH1* projPH)
3899 {
3900   //
3901   // Slope methode for the drift velocity
3902   //
3903
3904   // Constants
3905   const Float_t kDrWidth = AliTRDgeometry::DrThick();
3906   Int_t binmax           = 0;
3907   Int_t binmin           = 0;
3908   fPhd[0]                = 0.0;
3909   fPhd[1]                = 0.0;
3910   fPhd[2]                = 0.0;
3911   Int_t ju               = 0;
3912   fCurrentCoefE          = 0.0;
3913   fCurrentCoefE2         = 0.0;
3914   fCurrentCoef[0]        = 0.0;
3915   fCurrentCoef2[0]       = 0.0;
3916   TLine *line            = new TLine();
3917
3918   // Some variables
3919   TAxis   *xpph    = projPH->GetXaxis();
3920   Int_t    nbins   = xpph->GetNbins();
3921   Double_t lowedge = xpph->GetBinLowEdge(1);
3922   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
3923   Double_t widbins = (upedge - lowedge) / nbins;
3924   Double_t limit   = upedge + 0.5 * widbins;
3925   Bool_t put = kTRUE;
3926
3927   // Beginning of the signal
3928   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
3929   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
3930     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3931   }
3932   binmax = (Int_t) pentea->GetMaximumBin();
3933   if (binmax <= 1) {
3934     binmax = 2;
3935     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3936   }
3937   if (binmax >= nbins) {
3938     binmax = nbins-1;
3939     put = kFALSE;
3940     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3941   }
3942   pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
3943   Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
3944   Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
3945   Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
3946   Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
3947   if (l3P2am != 0) {
3948     fPhd[0] = -(l3P1am / (2 * l3P2am));
3949   }
3950   if(!fTakeTheMaxPH){
3951     if((l3P1am != 0.0) && (l3P2am != 0.0)){
3952       fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
3953     }
3954   }
3955   // Amplification region
3956   binmax = 0;
3957   ju     = 0;
3958   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
3959     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
3960       binmax = kbin;
3961       ju     = 1;
3962     }
3963   }
3964   if (binmax <= 1) {
3965     binmax = 2;
3966     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3967   }
3968   if (binmax >= nbins) {
3969     binmax = nbins-1;
3970     put = kFALSE;
3971     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
3972   }
3973   projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
3974   Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
3975   Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
3976   Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
3977   Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
3978   if (l3P2amf != 0) {
3979     fPhd[1] = -(l3P1amf / (2 * l3P2amf));
3980   }
3981   if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
3982     fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
3983   }
3984   if(fTakeTheMaxPH){
3985     fCurrentCoefE2 = fCurrentCoefE;
3986   }
3987   // Drift region
3988   TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
3989   for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
3990     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
3991   }
3992   binmin = 0;
3993   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
3994   if (binmin <= 1) {
3995     binmin = 2;
3996     AliInfo("Put the binmax from 1 to 2 to enable the fit");
3997   }
3998   if (binmin >= nbins) {
3999     binmin = nbins-1;
4000     put = kFALSE;
4001     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4002   }
4003   pente->Fit("pol2"
4004             ,"0MR"
4005             ,""
4006             ,TMath::Max(pente->GetBinCenter(binmin-1),             0.0)
4007             ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4008   Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4009   Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4010   Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4011   Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4012   if (l3P2dr != 0) {
4013     fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4014   }
4015   if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
4016     fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4017   }
4018   Float_t fPhdt0  = 0.0;
4019   Float_t t0Shift = 0.0;
4020   if(fTakeTheMaxPH) {
4021     fPhdt0 = fPhd[1];
4022     t0Shift = fT0Shift1;
4023   }
4024   else {
4025     fPhdt0 = fPhd[0];
4026     t0Shift = fT0Shift0;
4027   }
4028
4029   if ((fPhd[2] > fPhd[0]) &&
4030       (fPhd[2] > fPhd[1]) &&
4031       (fPhd[1] > fPhd[0]) &&
4032       (put)) {
4033     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4034     fNumberFitSuccess++;
4035
4036     if (fPhdt0 >= 0.0) {
4037       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4038       if (fCurrentCoef2[0] < -1.0) {
4039         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4040       }
4041     }
4042     else {
4043       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4044     }
4045
4046   }
4047   else {
4048     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
4049     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4050   }
4051
4052   if (fDebugLevel == 1) {
4053     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4054     cpentei->cd();
4055     projPH->Draw();
4056     line->SetLineColor(2);
4057     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4058     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4059     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4060     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
4061     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
4062     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
4063     AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4064     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4065     cpentei2->cd();
4066     pentea->Draw();
4067     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4068     cpentei3->cd();
4069     pente->Draw();
4070   }
4071   else {
4072     delete pentea;
4073     delete pente;
4074   }
4075 }
4076 //_____________________________________________________________________________
4077 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4078 {
4079   //
4080   // Slope methode but with polynomes de Lagrange
4081   //
4082
4083   // Constants
4084   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4085   Int_t binmax      = 0;
4086   Int_t binmin      = 0;
4087   Double_t    *x    = new Double_t[5];
4088   Double_t    *y    = new Double_t[5];
4089   x[0]              = 0.0;
4090   x[1]              = 0.0;
4091   x[2]              = 0.0;
4092   x[3]              = 0.0;
4093   x[4]              = 0.0;
4094   y[0]              = 0.0;
4095   y[1]              = 0.0;
4096   y[2]              = 0.0;
4097   y[3]              = 0.0;
4098   y[4]              = 0.0;
4099   fPhd[0]           = 0.0;
4100   fPhd[1]           = 0.0;
4101   fPhd[2]           = 0.0;
4102   Int_t ju          = 0;
4103   fCurrentCoefE     = 0.0;
4104   fCurrentCoefE2    = 1.0;
4105   fCurrentCoef[0]   = 0.0;
4106   fCurrentCoef2[0]  = 0.0;
4107   TLine *line = new TLine();
4108   TF1 * polynome = 0x0;
4109   TF1 * polynomea = 0x0;
4110   TF1 * polynomeb = 0x0;
4111   Double_t *c = 0x0;
4112
4113   // Some variables
4114   TAxis   *xpph    = projPH->GetXaxis();
4115   Int_t    nbins   = xpph->GetNbins();
4116   Double_t lowedge = xpph->GetBinLowEdge(1);
4117   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
4118   Double_t widbins = (upedge - lowedge) / nbins;
4119   Double_t limit   = upedge + 0.5 * widbins;
4120
4121
4122   Bool_t put = kTRUE;
4123
4124   // Beginning of the signal
4125   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4126   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
4127     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4128   }
4129
4130   binmax = (Int_t) pentea->GetMaximumBin();
4131
4132   Double_t minnn = 0.0;
4133   Double_t maxxx = 0.0;
4134
4135   Int_t kase = nbins-binmax;
4136
4137   switch(kase)
4138     {
4139     case 0:
4140       put = kFALSE;
4141       break;
4142     case 1:
4143       minnn = pentea->GetBinCenter(binmax-2);
4144       maxxx = pentea->GetBinCenter(binmax);
4145       x[0] = pentea->GetBinCenter(binmax-2);
4146       x[1] = pentea->GetBinCenter(binmax-1);
4147       x[2] = pentea->GetBinCenter(binmax);
4148       y[0] = pentea->GetBinContent(binmax-2);
4149       y[1] = pentea->GetBinContent(binmax-1);
4150       y[2] = pentea->GetBinContent(binmax);
4151       c = CalculPolynomeLagrange2(x,y);
4152       AliInfo("At the limit for beginning!");
4153       break;
4154     case 2:
4155       minnn = pentea->GetBinCenter(binmax-2);
4156       maxxx = pentea->GetBinCenter(binmax+1);
4157       x[0] = pentea->GetBinCenter(binmax-2);
4158       x[1] = pentea->GetBinCenter(binmax-1);
4159       x[2] = pentea->GetBinCenter(binmax);
4160       x[3] = pentea->GetBinCenter(binmax+1);
4161       y[0] = pentea->GetBinContent(binmax-2);
4162       y[1] = pentea->GetBinContent(binmax-1);
4163       y[2] = pentea->GetBinContent(binmax);
4164       y[3] = pentea->GetBinContent(binmax+1);
4165       c = CalculPolynomeLagrange3(x,y);
4166       break;
4167     default:
4168       switch(binmax){
4169       case 0:
4170         put = kFALSE;
4171         break;
4172       case 1:
4173         minnn = pentea->GetBinCenter(binmax);
4174         maxxx = pentea->GetBinCenter(binmax+2);
4175         x[0] = pentea->GetBinCenter(binmax);
4176         x[1] = pentea->GetBinCenter(binmax+1);
4177         x[2] = pentea->GetBinCenter(binmax+2);
4178         y[0] = pentea->GetBinContent(binmax);
4179         y[1] = pentea->GetBinContent(binmax+1);
4180         y[2] = pentea->GetBinContent(binmax+2);
4181         c = CalculPolynomeLagrange2(x,y);
4182         break;
4183       case 2:
4184         minnn = pentea->GetBinCenter(binmax-1);
4185         maxxx = pentea->GetBinCenter(binmax+2);
4186         x[0] = pentea->GetBinCenter(binmax-1);
4187         x[1] = pentea->GetBinCenter(binmax);
4188         x[2] = pentea->GetBinCenter(binmax+1);
4189         x[3] = pentea->GetBinCenter(binmax+2);
4190         y[0] = pentea->GetBinContent(binmax-1);
4191         y[1] = pentea->GetBinContent(binmax);
4192         y[2] = pentea->GetBinContent(binmax+1);
4193         y[3] = pentea->GetBinContent(binmax+2);
4194         c = CalculPolynomeLagrange3(x,y);
4195         break;
4196       default:
4197         minnn = pentea->GetBinCenter(binmax-2);
4198         maxxx = pentea->GetBinCenter(binmax+2);
4199         x[0] = pentea->GetBinCenter(binmax-2);
4200         x[1] = pentea->GetBinCenter(binmax-1);
4201         x[2] = pentea->GetBinCenter(binmax);
4202         x[3] = pentea->GetBinCenter(binmax+1);
4203         x[4] = pentea->GetBinCenter(binmax+2);
4204         y[0] = pentea->GetBinContent(binmax-2);
4205         y[1] = pentea->GetBinContent(binmax-1);
4206         y[2] = pentea->GetBinContent(binmax);
4207         y[3] = pentea->GetBinContent(binmax+1);
4208         y[4] = pentea->GetBinContent(binmax+2);
4209         c = CalculPolynomeLagrange4(x,y);
4210         break;
4211       }
4212       break;
4213     }
4214
4215
4216   if(put) {
4217     polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4218     polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4219
4220     Double_t step = (maxxx-minnn)/10000;
4221     Double_t l = minnn;
4222     Double_t maxvalue = 0.0;
4223     Double_t placemaximum = minnn;
4224     for(Int_t o = 0; o < 10000; o++){
4225       if(o == 0) maxvalue = polynomeb->Eval(l);
4226       if(maxvalue < (polynomeb->Eval(l))){
4227         maxvalue = polynomeb->Eval(l);
4228         placemaximum = l;
4229       }
4230       l += step;
4231     }
4232     fPhd[0] = placemaximum;
4233   }
4234
4235   // Amplification region
4236   binmax = 0;
4237   ju     = 0;
4238   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4239     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4240       binmax = kbin;
4241       ju     = 1;
4242     }
4243   }
4244
4245   Double_t minn = 0.0;
4246   Double_t maxx = 0.0;
4247   x[0] = 0.0;
4248   x[1] = 0.0;
4249   x[2] = 0.0;
4250   x[3] = 0.0;
4251   x[4] = 0.0;
4252   y[0] = 0.0;
4253   y[1] = 0.0;
4254   y[2] = 0.0;
4255   y[3] = 0.0;
4256   y[4] = 0.0;
4257
4258   Int_t    kase1 = nbins - binmax;
4259
4260   //Determination of minn and maxx
4261   //case binmax = nbins
4262   //pol2
4263   switch(kase1)
4264     {
4265     case 0:
4266       minn = projPH->GetBinCenter(binmax-2);
4267       maxx = projPH->GetBinCenter(binmax);
4268       x[0] = projPH->GetBinCenter(binmax-2);
4269       x[1] = projPH->GetBinCenter(binmax-1);
4270       x[2] = projPH->GetBinCenter(binmax);
4271       y[0] = projPH->GetBinContent(binmax-2);
4272       y[1] = projPH->GetBinContent(binmax-1);
4273       y[2] = projPH->GetBinContent(binmax);
4274       c = CalculPolynomeLagrange2(x,y);
4275       //AliInfo("At the limit for the drift!");
4276       break;
4277     case 1:
4278       minn = projPH->GetBinCenter(binmax-2);
4279       maxx = projPH->GetBinCenter(binmax+1);
4280       x[0] = projPH->GetBinCenter(binmax-2);
4281       x[1] = projPH->GetBinCenter(binmax-1);
4282       x[2] = projPH->GetBinCenter(binmax);
4283       x[3] = projPH->GetBinCenter(binmax+1);
4284       y[0] = projPH->GetBinContent(binmax-2);
4285       y[1] = projPH->GetBinContent(binmax-1);
4286       y[2] = projPH->GetBinContent(binmax);
4287       y[3] = projPH->GetBinContent(binmax+1);
4288       c = CalculPolynomeLagrange3(x,y);
4289       break;
4290     default:
4291       switch(binmax)
4292         {
4293         case 0:
4294           put = kFALSE;
4295           break;
4296         case 1:
4297           minn = projPH->GetBinCenter(binmax);
4298           maxx = projPH->GetBinCenter(binmax+2);
4299           x[0] = projPH->GetBinCenter(binmax);
4300           x[1] = projPH->GetBinCenter(binmax+1);
4301           x[2] = projPH->GetBinCenter(binmax+2);
4302           y[0] = projPH->GetBinContent(binmax);
4303           y[1] = projPH->GetBinContent(binmax+1);
4304           y[2] = projPH->GetBinContent(binmax+2);
4305           c = CalculPolynomeLagrange2(x,y);
4306           break;
4307         case 2:
4308           minn = projPH->GetBinCenter(binmax-1);
4309           maxx = projPH->GetBinCenter(binmax+2);
4310           x[0] = projPH->GetBinCenter(binmax-1);
4311           x[1] = projPH->GetBinCenter(binmax);
4312           x[2] = projPH->GetBinCenter(binmax+1);
4313           x[3] = projPH->GetBinCenter(binmax+2);
4314           y[0] = projPH->GetBinContent(binmax-1);
4315           y[1] = projPH->GetBinContent(binmax);
4316           y[2] = projPH->GetBinContent(binmax+1);
4317           y[3] = projPH->GetBinContent(binmax+2);
4318           c = CalculPolynomeLagrange3(x,y);
4319           break;
4320         default:
4321           minn = projPH->GetBinCenter(binmax-2);
4322           maxx = projPH->GetBinCenter(binmax+2);
4323           x[0] = projPH->GetBinCenter(binmax-2);
4324           x[1] = projPH->GetBinCenter(binmax-1);
4325           x[2] = projPH->GetBinCenter(binmax);
4326           x[3] = projPH->GetBinCenter(binmax+1);
4327           x[4] = projPH->GetBinCenter(binmax+2);
4328           y[0] = projPH->GetBinContent(binmax-2);
4329           y[1] = projPH->GetBinContent(binmax-1);
4330           y[2] = projPH->GetBinContent(binmax);
4331           y[3] = projPH->GetBinContent(binmax+1);
4332           y[4] = projPH->GetBinContent(binmax+2);
4333           c = CalculPolynomeLagrange4(x,y);
4334           break;
4335         }
4336       break;
4337     }
4338
4339   if(put) {
4340     polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4341     polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4342
4343     Double_t step = (maxx-minn)/1000;
4344     Double_t l = minn;
4345     Double_t maxvalue = 0.0;
4346     Double_t placemaximum = minn;
4347     for(Int_t o = 0; o < 1000; o++){
4348       if(o == 0) maxvalue = polynomea->Eval(l);
4349       if(maxvalue < (polynomea->Eval(l))){
4350         maxvalue = polynomea->Eval(l);
4351         placemaximum = l;
4352       }
4353       l += step;
4354     }
4355     fPhd[1] = placemaximum;
4356   }
4357
4358   // Drift region
4359   TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4360   for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
4361     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4362   }
4363   binmin = 0;
4364   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4365
4366   //should not happen
4367   if (binmin <= 1) {
4368     binmin = 2;
4369     put = 1;
4370     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4371   }
4372
4373   //check
4374   if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4375     AliInfo("Too many fluctuations at the end!");
4376     put = kFALSE;
4377   }
4378   if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4379     AliInfo("Too many fluctuations at the end!");
4380     put = kFALSE;
4381   }
4382   if(pente->GetBinContent(binmin+1)==0){
4383     AliInfo("No entries for the next bin!");
4384     pente->SetBinContent(binmin,0);
4385     if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4386   }
4387
4388
4389   x[0] = 0.0;
4390   x[1] = 0.0;
4391   x[2] = 0.0;
4392   x[3] = 0.0;
4393   x[4] = 0.0;
4394   y[0] = 0.0;
4395   y[1] = 0.0;
4396   y[2] = 0.0;
4397   y[3] = 0.0;
4398   y[4] = 0.0;
4399   Double_t min = 0.0;
4400   Double_t max = 0.0;
4401   Bool_t case1 = kFALSE;
4402   Bool_t case2 = kFALSE;
4403   Bool_t case4 = kFALSE;
4404
4405   //Determination of min and max
4406   //case binmin <= nbins-3
4407   //pol4 case 3
4408   if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4409     min = pente->GetBinCenter(binmin-2);
4410     max = pente->GetBinCenter(binmin+2);
4411     x[0] = pente->GetBinCenter(binmin-2);
4412     x[1] = pente->GetBinCenter(binmin-1);
4413     x[2] = pente->GetBinCenter(binmin);
4414     x[3] = pente->GetBinCenter(binmin+1);
4415     x[4] = pente->GetBinCenter(binmin+2);
4416     y[0] = pente->GetBinContent(binmin-2);
4417     y[1] = pente->GetBinContent(binmin-1);
4418     y[2] = pente->GetBinContent(binmin);
4419     y[3] = pente->GetBinContent(binmin+1);
4420     y[4] = pente->GetBinContent(binmin+2);
4421     //Calcul the polynome de Lagrange
4422     c = CalculPolynomeLagrange4(x,y);
4423     //richtung +/-
4424     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4425        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4426       //AliInfo("polynome 4 false 1");
4427       put = kFALSE;
4428     }
4429     if(((binmin+3) <= (nbins-1)) &&
4430        (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4431        ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4432        (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4433       AliInfo("polynome 4 false 2");
4434       put = kFALSE;
4435     }
4436     // poly 3
4437     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4438        (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4439       //AliInfo("polynome 4 case 1");
4440       case1 = kTRUE;
4441     }
4442     if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4443        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4444       //AliInfo("polynome 4 case 4");
4445       case4 = kTRUE;
4446     }
4447
4448   }
4449   //case binmin = nbins-2
4450   //pol3 case 1
4451   if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4452      (case1)){
4453     min = pente->GetBinCenter(binmin-2);
4454     max = pente->GetBinCenter(binmin+1);
4455     x[0] = pente->GetBinCenter(binmin-2);
4456     x[1] = pente->GetBinCenter(binmin-1);
4457     x[2] = pente->GetBinCenter(binmin);
4458     x[3] = pente->GetBinCenter(binmin+1);
4459     y[0] = pente->GetBinContent(binmin-2);
4460     y[1] = pente->GetBinContent(binmin-1);
4461     y[2] = pente->GetBinContent(binmin);
4462     y[3] = pente->GetBinContent(binmin+1);
4463     //Calcul the polynome de Lagrange
4464     c = CalculPolynomeLagrange3(x,y);
4465     //richtung +: nothing
4466     //richtung -
4467     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4468       //AliInfo("polynome 3- case 2");
4469       case2 = kTRUE;
4470     }
4471   }
4472   //pol3 case 4
4473   if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4474      (case4)){
4475     min = pente->GetBinCenter(binmin-1);
4476     max = pente->GetBinCenter(binmin+2);
4477     x[0] = pente->GetBinCenter(binmin-1);
4478     x[1] = pente->GetBinCenter(binmin);
4479     x[2] = pente->GetBinCenter(binmin+1);
4480     x[3] = pente->GetBinCenter(binmin+2);
4481     y[0] = pente->GetBinContent(binmin-1);
4482     y[1] = pente->GetBinContent(binmin);
4483     y[2] = pente->GetBinContent(binmin+1);
4484     y[3] = pente->GetBinContent(binmin+2);
4485     //Calcul the polynome de Lagrange
4486     c = CalculPolynomeLagrange3(x,y);
4487     //richtung +
4488     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4489       //AliInfo("polynome 3+ case 2");
4490       case2 = kTRUE;
4491     }
4492   }
4493   //pol2 case 5
4494   if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4495     min = pente->GetBinCenter(binmin);
4496     max = pente->GetBinCenter(binmin+2);
4497     x[0] = pente->GetBinCenter(binmin);
4498     x[1] = pente->GetBinCenter(binmin+1);
4499     x[2] = pente->GetBinCenter(binmin+2);
4500     y[0] = pente->GetBinContent(binmin);
4501     y[1] = pente->GetBinContent(binmin+1);
4502     y[2] = pente->GetBinContent(binmin+2);
4503     //Calcul the polynome de Lagrange
4504     c = CalculPolynomeLagrange2(x,y);
4505     //richtung +
4506     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4507       //AliInfo("polynome 2+ false");
4508       put = kFALSE;
4509     }
4510   }
4511   //pol2 case 2
4512   if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4513      (case2)){
4514     min = pente->GetBinCenter(binmin-1);
4515     max = pente->GetBinCenter(binmin+1);
4516     x[0] = pente->GetBinCenter(binmin-1);
4517     x[1] = pente->GetBinCenter(binmin);
4518     x[2] = pente->GetBinCenter(binmin+1);
4519     y[0] = pente->GetBinContent(binmin-1);
4520     y[1] = pente->GetBinContent(binmin);
4521     y[2] = pente->GetBinContent(binmin+1);
4522     //Calcul the polynome de Lagrange
4523     c = CalculPolynomeLagrange2(x,y);
4524     //richtung +: nothing
4525     //richtung -: nothing
4526   }
4527   //case binmin = nbins-1
4528   //pol2 case 0
4529   if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4530     min = pente->GetBinCenter(binmin-2);
4531     max = pente->GetBinCenter(binmin);
4532     x[0] = pente->GetBinCenter(binmin-2);
4533     x[1] = pente->GetBinCenter(binmin-1);
4534     x[2] = pente->GetBinCenter(binmin);
4535     y[0] = pente->GetBinContent(binmin-2);
4536     y[1] = pente->GetBinContent(binmin-1);
4537     y[2] = pente->GetBinContent(binmin);
4538     //Calcul the polynome de Lagrange
4539     c = CalculPolynomeLagrange2(x,y);
4540     //AliInfo("At the limit for the drift!");
4541     //fluctuation too big!
4542     //richtung +: nothing
4543     //richtung -
4544     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4545       //AliInfo("polynome 2- false ");
4546       put = kFALSE;
4547     }
4548   }
4549   if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4550     put = kFALSE;
4551     AliInfo("At the limit for the drift and not usable!");
4552   }
4553
4554   //pass
4555   if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4556     put = kFALSE;
4557     AliInfo("For the drift...problem!");
4558   }
4559   //pass but should not happen
4560   if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4561     put = kFALSE;
4562     AliInfo("For the drift...problem!");
4563   }
4564
4565   if(put) {
4566     polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4567     polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4568     //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4569     Double_t step = (max-min)/1000;
4570     Double_t l = min;
4571     Double_t minvalue = 0.0;
4572     Double_t placeminimum = min;
4573     for(Int_t o = 0; o < 1000; o++){
4574       if(o == 0) minvalue = polynome->Eval(l);
4575       if(minvalue > (polynome->Eval(l))){
4576         minvalue = polynome->Eval(l);
4577         placeminimum = l;
4578       }
4579       l += step;
4580     }
4581     fPhd[2] = placeminimum;
4582   }
4583   //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4584   if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4585   if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4586
4587   Float_t fPhdt0  = 0.0;
4588   Float_t t0Shift = 0.0;
4589   if(fTakeTheMaxPH) {
4590     fPhdt0 = fPhd[1];
4591     t0Shift = fT0Shift1;
4592   }
4593   else {
4594     fPhdt0 = fPhd[0];
4595     t0Shift = fT0Shift0;
4596   }
4597
4598   if ((fPhd[2] > fPhd[0]) &&
4599       (fPhd[2] > fPhd[1]) &&
4600       (fPhd[1] > fPhd[0]) &&
4601       (put)) {
4602     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4603     fNumberFitSuccess++;
4604     if (fPhdt0 >= 0.0) {
4605       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4606       if (fCurrentCoef2[0] < -1.0) {
4607         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4608       }
4609     }
4610     else {
4611       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4612     }
4613   }
4614   else {
4615     if((fPhd[1] > fPhd[0]) &&
4616        (put)) {
4617       if (fPhdt0 >= 0.0) {
4618         fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4619         if (fCurrentCoef2[0] < -1.0) {
4620           fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4621         }
4622       }
4623       else {
4624         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4625       }
4626     }
4627     else{
4628       fCurrentCoef[0]      = -TMath::Abs(fCurrentCoef[1]);
4629       fCurrentCoef2[0]     = fCurrentCoef2[1] + 100.0;
4630       //printf("Fit failed!\n");
4631     }
4632   }
4633
4634   if (fDebugLevel == 1) {
4635     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4636     cpentei->cd();
4637     projPH->Draw();
4638     line->SetLineColor(2);
4639     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4640     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4641     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4642     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
4643     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
4644     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
4645     AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4646     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4647     cpentei2->cd();
4648     pentea->Draw();
4649     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4650     cpentei3->cd();
4651     pente->Draw();
4652   }
4653   else {
4654     if(pentea) delete pentea;
4655     if(pente) delete pente;
4656     if(polynome) delete polynome;
4657     if(polynomea) delete polynomea;
4658     if(polynomeb) delete polynomeb;
4659     if(x) delete [] x;
4660     if(y) delete [] y;
4661     if(c) delete [] c;
4662     if(line) delete line;
4663
4664   }
4665
4666   //Provisoire
4667   //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4668   //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4669
4670   projPH->SetDirectory(0);
4671
4672 }
4673
4674 //_____________________________________________________________________________
4675 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4676 {
4677   //
4678   // Fit methode for the drift velocity
4679   //
4680
4681   // Constants
4682   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4683
4684   // Some variables
4685   TAxis   *xpph   = projPH->GetXaxis();
4686   Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4687
4688   TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4689   fPH->SetParameter(0,0.469);     // Scaling
4690   fPH->SetParameter(1,0.18);      // Start
4691   fPH->SetParameter(2,0.0857325); // AR
4692   fPH->SetParameter(3,1.89);      // DR
4693   fPH->SetParameter(4,0.08);      // QA/QD
4694   fPH->SetParameter(5,0.0);       // Baseline
4695
4696   TLine *line = new TLine();
4697
4698   fCurrentCoef[0]     = 0.0;
4699   fCurrentCoef2[0]    = 0.0;
4700   fCurrentCoefE       = 0.0;
4701   fCurrentCoefE2      = 0.0;
4702
4703   if (idect%fFitPHPeriode == 0) {
4704
4705     AliInfo(Form("The detector %d will be fitted",idect));
4706     fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4707     fPH->SetParameter(1,fPhd[0] - 0.1);                                                                 // Start
4708     fPH->SetParameter(2,fPhd[1] - fPhd[0]);                                                             // AR
4709     fPH->SetParameter(3,fPhd[2] - fPhd[1]);                                                             // DR
4710     fPH->SetParameter(4,0.225);                                                                         // QA/QD
4711     fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4712
4713     if (fDebugLevel != 1) {
4714       projPH->Fit(fPH,"0M","",0.0,upedge);
4715     }
4716     else {
4717       TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4718       cpente->cd();
4719       projPH->Fit(fPH,"M+","",0.0,upedge);
4720       projPH->Draw("E0");
4721       line->SetLineColor(4);
4722       line->DrawLine(fPH->GetParameter(1)
4723                     ,0
4724                     ,fPH->GetParameter(1)
4725                     ,projPH->GetMaximum());
4726       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4727                     ,0
4728                     ,fPH->GetParameter(1)+fPH->GetParameter(2)
4729                     ,projPH->GetMaximum());
4730       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4731                     ,0
4732                     ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4733                     ,projPH->GetMaximum());
4734     }
4735
4736     if (fPH->GetParameter(3) != 0) {
4737       fNumberFitSuccess++;
4738       fCurrentCoef[0]    = kDrWidth / (fPH->GetParameter(3));
4739       fCurrentCoefE      = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4740       fCurrentCoef2[0]   = fPH->GetParameter(1);
4741       fCurrentCoefE2     = fPH->GetParError(1);
4742     }
4743     else {
4744       fCurrentCoef[0]     = -TMath::Abs(fCurrentCoef[1]);
4745       fCurrentCoef2[0]    = fCurrentCoef2[1] + 100.0;
4746     }
4747
4748   }
4749   else {
4750
4751     // Put the default value
4752     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
4753     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4754   }
4755
4756   if (fDebugLevel != 1) {
4757     delete fPH;
4758   }
4759
4760 }
4761 //_____________________________________________________________________________
4762 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
4763 {
4764   //
4765   // Fit methode for the sigma of the pad response function
4766   //
4767
4768   TVectorD param(3);
4769
4770   fCurrentCoef[0]  = 0.0;
4771   fCurrentCoefE = 0.0;
4772
4773   Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
4774
4775   if(ret == -4){
4776     fCurrentCoef[0] = -fCurrentCoef[1];
4777     return kFALSE;
4778   }
4779   else {
4780     fNumberFitSuccess++;
4781     fCurrentCoef[0] = param[2];
4782     fCurrentCoefE   = ret;
4783     return kTRUE;
4784   }
4785 }
4786 //_____________________________________________________________________________
4787 Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, Bool_t bError)
4788 {
4789   //
4790   // Fit methode for the sigma of the pad response function
4791   //
4792
4793   //We should have at least 3 points
4794   if(nBins <=3) return -4.0;
4795
4796   TLinearFitter fitter(3,"pol2");
4797   fitter.StoreData(kFALSE);
4798   fitter.ClearPoints();
4799   TVectorD  par(3);
4800   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
4801   Float_t entries = 0;
4802   Int_t   nbbinwithentries = 0;
4803   for (Int_t i=0; i<nBins; i++){
4804     entries+=arraye[i];
4805     if(arraye[i] > 15) nbbinwithentries++;
4806     //printf("entries for i %d: %f\n",i,arraye[i]);
4807   }
4808   if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
4809   //printf("entries %f\n",entries);
4810   //printf("nbbinwithentries %d\n",nbbinwithentries);
4811
4812   Int_t npoints=0;
4813   Float_t errorm = 0.0;
4814   Float_t errorn = 0.0;
4815   Float_t error  = 0.0;
4816
4817   //
4818   for (Int_t ibin=0;ibin<nBins; ibin++){
4819       Float_t entriesI = arraye[ibin];
4820       Float_t valueI   = arraym[ibin];
4821       Double_t xcenter = 0.0;
4822       Float_t  val     = 0.0;
4823       if ((entriesI>15) && (valueI>0.0)){
4824         xcenter = xMin+(ibin+0.5)*binWidth;
4825         errorm   = 0.0;
4826         errorn   = 0.0;
4827         error    = 0.0;
4828         if(!bError){
4829           if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
4830           if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
4831           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4832         }
4833         else{
4834           if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
4835           if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
4836           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
4837         }
4838         if(error == 0.0) continue;
4839         val      = TMath::Log(Float_t(valueI));
4840         fitter.AddPoint(&xcenter,val,error);
4841         npoints++;
4842       }
4843
4844       if(fDebugLevel > 1){
4845
4846       if ( !fDebugStreamer ) {
4847         //debug stream
4848         TDirectory *backup = gDirectory;
4849         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4850         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
4851       }
4852
4853       Int_t    detector     = fCountDet;
4854       Int_t    layer        = GetLayer(fCountDet);
4855       Int_t    group        = ibin;
4856
4857       (* fDebugStreamer) << "FitGausMIFill"<<
4858         "detector="<<detector<<
4859         "layer="<<layer<<
4860         "nbins="<<nBins<<
4861         "group="<<group<<
4862         "entriesI="<<entriesI<<
4863         "valueI="<<valueI<<
4864         "val="<<val<<
4865         "xcenter="<<xcenter<<
4866         "errorm="<<errorm<<
4867         "errorn="<<errorn<<
4868         "error="<<error<<
4869         "bError="<<bError<<
4870         "\n";
4871     }
4872
4873   }
4874
4875   if(npoints <=3) return -4.0;
4876
4877   Double_t chi2 = 0;
4878   if (npoints>3){
4879     fitter.Eval();
4880     fitter.GetParameters(par);
4881     chi2 = fitter.GetChisquare()/Float_t(npoints);
4882
4883
4884     if (!param)  param  = new TVectorD(3);
4885     if(par[2] == 0.0) return -4.0;
4886     Double_t  x      = TMath::Sqrt(TMath::Abs(-2*par[2]));
4887     Double_t deltax = (fitter.GetParError(2))/x;
4888     Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
4889     chi2 = errorparam2;
4890
4891     (*param)[1] = par[1]/(-2.*par[2]);
4892     (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
4893     Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
4894     if ( lnparam0>307 ) return -4;
4895     (*param)[0] = TMath::Exp(lnparam0);
4896
4897     if(fDebugLevel > 1){
4898
4899       if ( !fDebugStreamer ) {
4900         //debug stream
4901         TDirectory *backup = gDirectory;
4902         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
4903         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
4904       }
4905
4906       Int_t    detector     = fCountDet;
4907       Int_t    layer        = GetLayer(fCountDet);
4908
4909
4910       (* fDebugStreamer) << "FitGausMIFit"<<
4911         "detector="<<detector<<
4912         "layer="<<layer<<
4913         "nbins="<<nBins<<
4914         "errorsigma="<<chi2<<
4915         "mean="<<(*param)[1]<<
4916         "sigma="<<(*param)[2]<<
4917         "constant="<<(*param)[0]<<
4918         "\n";
4919     }
4920   }
4921
4922   if((chi2/(*param)[2]) > 0.1){
4923     if(bError){
4924       chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
4925     }
4926     else return -4.0;
4927   }
4928
4929   if(fDebugLevel == 1){
4930     TString name("PRF");
4931     name += (Int_t)xMin;
4932     name += (Int_t)xMax;
4933     TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
4934     c1->cd();
4935     name += "histo";
4936     TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
4937     for(Int_t k = 0; k < nBins; k++){
4938       histo->SetBinContent(k+1,arraym[k]);
4939       histo->SetBinError(k+1,arrayme[k]);
4940     }
4941     histo->Draw();
4942     name += "functionf";
4943     TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
4944     f1->SetParameter(0, (*param)[0]);
4945     f1->SetParameter(1, (*param)[1]);
4946     f1->SetParameter(2, (*param)[2]);
4947     f1->Draw("same");
4948   }
4949
4950
4951   return chi2;
4952
4953 }
4954 //_____________________________________________________________________________
4955 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
4956 {
4957   //
4958   // Fit methode for the sigma of the pad response function
4959   //
4960
4961   fCurrentCoef[0]  = 0.0;
4962   fCurrentCoefE = 0.0;
4963
4964   if (fDebugLevel != 1) {
4965     projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
4966   }
4967   else {
4968     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4969     cfit->cd();
4970     projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
4971     projPRF->Draw();
4972   }
4973   fCurrentCoef[0]  = projPRF->GetFunction("gaus")->GetParameter(2);
4974   fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
4975   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4976   else {
4977     fNumberFitSuccess++;
4978   }
4979 }
4980 //_____________________________________________________________________________
4981 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
4982 {
4983   //
4984   // Fit methode for the sigma of the pad response function
4985   //
4986   fCurrentCoef[0]   = 0.0;
4987   fCurrentCoefE  = 0.0;
4988   if (fDebugLevel == 1) {
4989     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
4990     cfit->cd();
4991     projPRF->Draw();
4992   }
4993   fCurrentCoef[0] = projPRF->GetRMS();
4994   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
4995   else {
4996     fNumberFitSuccess++;
4997   }
4998 }
4999 //_____________________________________________________________________________
5000 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5001 {
5002   //
5003   // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5004   //
5005
5006   TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5007
5008
5009   Int_t   nbins    = (Int_t)(nybins/(2*nbg));
5010   Float_t lowedge  = -3.0*nbg;
5011   Float_t upedge   = lowedge + 3.0;
5012   Int_t   offset   = 0;
5013   Int_t   npoints  = 0;
5014   Double_t xvalues = -0.2*nbg+0.1;
5015   Double_t y       = 0.0;
5016   Int_t   total    = 2*nbg;
5017
5018
5019   for(Int_t k = 0; k < total; k++){
5020     if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5021       npoints++;
5022       y = fCurrentCoef[0]*fCurrentCoef[0];
5023       linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5024     }
5025
5026     if(fDebugLevel > 1){
5027
5028       if ( !fDebugStreamer ) {
5029         //debug stream
5030         TDirectory *backup = gDirectory;
5031         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5032         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5033       }
5034
5035       Int_t    detector     = fCountDet;
5036       Int_t    layer        = GetLayer(fCountDet);
5037       Int_t    nbtotal      = total;
5038       Int_t    group        = k;
5039       Float_t  low          = lowedge;
5040       Float_t  up           = upedge;
5041       Float_t  tnp          = xvalues;
5042       Float_t  wid          = fCurrentCoef[0];
5043       Float_t  widfE        = fCurrentCoefE;
5044
5045       (* fDebugStreamer) << "FitTnpRange0"<<
5046         "detector="<<detector<<
5047         "layer="<<layer<<
5048         "nbtotal="<<nbtotal<<
5049         "group="<<group<<
5050         "low="<<low<<
5051         "up="<<up<<
5052         "offset="<<offset<<
5053         "tnp="<<tnp<<
5054         "wid="<<wid<<
5055         "widfE="<<widfE<<
5056         "\n";
5057     }
5058
5059     offset  += nbins;
5060     lowedge += 3.0;
5061     upedge  += 3.0;
5062     xvalues += 0.2;
5063
5064   }
5065
5066   fCurrentCoefE = 0.0;
5067   fCurrentCoef[0] = 0.0;
5068
5069   //printf("npoints\n",npoints);
5070
5071   if(npoints < 3){
5072     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5073   }
5074   else{
5075
5076     TVectorD pars0;
5077     linearfitter.Eval();
5078     linearfitter.GetParameters(pars0);
5079     Double_t pointError0  =  TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5080     Double_t errorsx0     =  linearfitter.GetParError(2)*pointError0;
5081     Double_t min0         = 0.0;
5082     Double_t ermin0       = 0.0;
5083     //Double_t prfe0      = 0.0;
5084     Double_t prf0         = 0.0;
5085     if((pars0[2] > 0.0) && (pars0[1] != 0.0)) {
5086       min0 = -pars0[1]/(2*pars0[2]);
5087       ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5088       prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5089       if(prf0 > 0.0) {
5090         /*
5091           prfe0 = linearfitter->GetParError(0)*pointError0
5092           +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5093           +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5094           prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5095           fCurrentCoefE   = (Float_t) prfe0;
5096         */
5097         fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5098       }
5099       else{
5100         fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5101       }
5102     }
5103     else {
5104       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5105     }
5106
5107     if(fDebugLevel > 1){
5108
5109       if ( !fDebugStreamer ) {
5110         //debug stream
5111         TDirectory *backup = gDirectory;
5112         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5113         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5114       }
5115
5116       Int_t    detector     = fCountDet;
5117       Int_t    layer        = GetLayer(fCountDet);
5118       Int_t    nbtotal      = total;
5119       Double_t colsize[6]   = {0.635,0.665,0.695,0.725,0.755,0.785};
5120       Double_t sigmax       = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5121
5122       (* fDebugStreamer) << "FitTnpRange1"<<
5123         "detector="<<detector<<
5124         "layer="<<layer<<
5125         "nbtotal="<<nbtotal<<
5126         "par0="<<pars0[0]<<
5127         "par1="<<pars0[1]<<
5128         "par2="<<pars0[2]<<
5129         "npoints="<<npoints<<
5130         "sigmax="<<sigmax<<
5131         "tan="<<min0<<
5132         "sigmaprf="<<fCurrentCoef[0]<<
5133         "sigprf="<<fCurrentCoef[1]<<
5134         "\n";
5135     }
5136
5137   }
5138
5139 }
5140 //_____________________________________________________________________________
5141 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5142 {
5143   //
5144   // Only mean methode for the gain factor
5145   //
5146
5147   fCurrentCoef[0] = mean;
5148   fCurrentCoefE   = 0.0;
5149   if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5150   if (fDebugLevel == 1) {
5151     TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5152     cpmean->cd();
5153     projch->Draw();
5154   }
5155   CalculChargeCoefMean(kTRUE);
5156   fNumberFitSuccess++;
5157 }
5158 //_____________________________________________________________________________
5159 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5160 {
5161   //
5162   // mean w methode for the gain factor
5163   //
5164
5165   //Number of bins
5166   Int_t nybins = projch->GetNbinsX();
5167
5168   //The weight function
5169   Double_t a = 0.00228515;
5170   Double_t b = -0.00231487;
5171   Double_t c = 0.00044298;
5172   Double_t d = -0.00379239;
5173   Double_t e = 0.00338349;
5174
5175 //   0 |0.00228515
5176 //    1 |-0.00231487
5177 //    2 |0.00044298
5178 //    3 |-0.00379239
5179 //    4 |0.00338349
5180
5181
5182
5183   //A arbitrary error for the moment
5184   fCurrentCoefE = 0.0;
5185   fCurrentCoef[0] = 0.0;
5186
5187   //Calcul
5188   Double_t sumw = 0.0;
5189   Double_t sum = 0.0;
5190   Float_t sumAll   = (Float_t) nentries;
5191   Int_t sumCurrent = 0;
5192   for(Int_t k = 0; k <nybins; k++){
5193     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5194     if (fraction>0.95) break;
5195     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5196       e*fraction*fraction*fraction*fraction;
5197     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5198     sum  += weight*projch->GetBinContent(k+1);
5199     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5200     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5201   }
5202   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5203
5204   if (fDebugLevel == 1) {
5205     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5206     cpmeanw->cd();
5207     projch->Draw();
5208   }
5209   fNumberFitSuccess++;
5210   CalculChargeCoefMean(kTRUE);
5211 }
5212 //_____________________________________________________________________________
5213 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5214 {
5215   //
5216   // mean w methode for the gain factor
5217   //
5218
5219   //Number of bins
5220   Int_t nybins = projch->GetNbinsX();
5221
5222   //The weight function
5223   Double_t a = 0.00228515;
5224   Double_t b = -0.00231487;
5225   Double_t c = 0.00044298;
5226   Double_t d = -0.00379239;
5227   Double_t e = 0.00338349;
5228
5229 //   0 |0.00228515
5230 //    1 |-0.00231487
5231 //    2 |0.00044298
5232 //    3 |-0.00379239
5233 //    4 |0.00338349
5234
5235
5236
5237   //A arbitrary error for the moment
5238   fCurrentCoefE = 0.0;
5239   fCurrentCoef[0] = 0.0;
5240
5241   //Calcul
5242   Double_t sumw = 0.0;
5243   Double_t sum = 0.0;
5244   Int_t sumCurrent = 0;
5245   for(Int_t k = 0; k <nybins; k++){
5246     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5247     if (fraction>0.95) break;
5248     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5249       e*fraction*fraction*fraction*fraction;
5250     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5251     sum  += weight*projch->GetBinContent(k+1);
5252     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5253     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5254   }
5255   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5256
5257   if (fDebugLevel == 1) {
5258     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5259     cpmeanw->cd();
5260     projch->Draw();
5261   }
5262   fNumberFitSuccess++;
5263 }
5264 //_____________________________________________________________________________
5265 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5266 {
5267   //
5268   // Fit methode for the gain factor
5269   //
5270
5271   fCurrentCoef[0]  = 0.0;
5272   fCurrentCoefE    = 0.0;
5273   Double_t chisqrl = 0.0;
5274   Double_t chisqrg = 0.0;
5275   Double_t chisqr  = 0.0;
5276   TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5277
5278   projch->Fit("landau","0",""
5279              ,(Double_t) mean/fBeginFitCharge
5280              ,projch->GetBinCenter(projch->GetNbinsX()));
5281   Double_t l3P0         = projch->GetFunction("landau")->GetParameter(0);
5282   Double_t l3P1         = projch->GetFunction("landau")->GetParameter(1);
5283   Double_t l3P2         = projch->GetFunction("landau")->GetParameter(2);
5284   chisqrl = projch->GetFunction("landau")->GetChisquare();
5285
5286   projch->Fit("gaus","0",""
5287               ,(Double_t) mean/fBeginFitCharge
5288               ,projch->GetBinCenter(projch->GetNbinsX()));
5289   Double_t g3P0         = projch->GetFunction("gaus")->GetParameter(0);
5290   Double_t g3P2         = projch->GetFunction("gaus")->GetParameter(2);
5291   chisqrg = projch->GetFunction("gaus")->GetChisquare();
5292
5293   fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5294   if (fDebugLevel != 1) {
5295     projch->Fit("fLandauGaus","0",""
5296                 ,(Double_t) mean/fBeginFitCharge
5297                 ,projch->GetBinCenter(projch->GetNbinsX()));
5298     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5299   }
5300   else  {
5301     TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5302     cp->cd();
5303     projch->Fit("fLandauGaus","+",""
5304                 ,(Double_t) mean/fBeginFitCharge
5305                 ,projch->GetBinCenter(projch->GetNbinsX()));
5306     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5307     projch->Draw();
5308     fLandauGaus->Draw("same");
5309   }
5310
5311   if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5312     //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5313     fNumberFitSuccess++;
5314     CalculChargeCoefMean(kTRUE);
5315     fCurrentCoef[0]  = projch->GetFunction("fLandauGaus")->GetParameter(1);
5316     fCurrentCoefE    = projch->GetFunction("fLandauGaus")->GetParError(1);
5317   }
5318   else {
5319     CalculChargeCoefMean(kFALSE);
5320     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5321   }
5322
5323   if (fDebugLevel != 1) {
5324     delete fLandauGaus;
5325   }
5326
5327 }
5328 //_____________________________________________________________________________
5329 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5330 {
5331   //
5332   // Fit methode for the gain factor more time consuming
5333   //
5334
5335
5336   //Some parameters to initialise
5337   Double_t widthLandau, widthGaus, mPV, integral;
5338   Double_t chisquarel = 0.0;
5339   Double_t chisquareg = 0.0;
5340   projch->Fit("landau","0M+",""
5341               ,(Double_t) mean/6
5342               ,projch->GetBinCenter(projch->GetNbinsX()));
5343   widthLandau  = projch->GetFunction("landau")->GetParameter(2);
5344   chisquarel = projch->GetFunction("landau")->GetChisquare();
5345   projch->Fit("gaus","0M+",""
5346               ,(Double_t) mean/6
5347               ,projch->GetBinCenter(projch->GetNbinsX()));
5348   widthGaus    = projch->GetFunction("gaus")->GetParameter(2);
5349   chisquareg = projch->GetFunction("gaus")->GetChisquare();
5350
5351   mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5352   integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5353
5354   // Setting fit range and start values
5355   Double_t fr[2];
5356   //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5357   //Double_t sv[4]   = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5358   Double_t sv[4]   = { widthLandau, mPV, integral, widthGaus};
5359   Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5360   Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5361   Double_t fp[4]   = { 1.0, 1.0, 1.0, 1.0 };
5362   Double_t fpe[4]  = { 1.0, 1.0, 1.0, 1.0 };
5363   fr[0]            = 0.3 * mean;
5364   fr[1]            = 3.0 * mean;
5365   fCurrentCoef[0]  = 0.0;
5366   fCurrentCoefE    = 0.0;
5367
5368   Double_t chisqr;
5369   Int_t    ndf;
5370   TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5371                                 ,&pllo[0],&plhi[0]
5372                                 ,&fp[0],&fpe[0]
5373                                 ,&chisqr,&ndf);
5374
5375   Double_t projchPeak;
5376   Double_t projchFWHM;
5377   LanGauPro(fp,projchPeak,projchFWHM);
5378
5379   if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5380     //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5381     fNumberFitSuccess++;
5382     CalculChargeCoefMean(kTRUE);
5383     fCurrentCoef[0]  = fp[1];
5384     fCurrentCoefE = fpe[1];
5385     //chargeCoefE2 = chisqr;
5386   }
5387   else {
5388     CalculChargeCoefMean(kFALSE);
5389     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5390   }
5391   if (fDebugLevel == 1) {
5392     AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5393     TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5394     cpy->cd();
5395     projch->Draw();
5396     fitsnr->Draw("same");
5397   }
5398   else {
5399     delete fitsnr;
5400   }
5401 }
5402 //_____________________________________________________________________________
5403 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
5404 {
5405   //
5406   // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5407   //
5408   Double_t *c = new Double_t[5];
5409   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5410   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5411   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5412
5413   c[4] = 0.0;
5414   c[3] = 0.0;
5415   c[2] = x0+x1+x2;
5416   c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5417   c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5418
5419   return c;
5420
5421
5422 }
5423
5424 //_____________________________________________________________________________
5425 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
5426 {
5427   //
5428   // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5429   //
5430   Double_t *c = new Double_t[5];
5431   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5432   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5433   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5434   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5435
5436   c[4] = 0.0;
5437   c[3] = x0+x1+x2+x3;
5438   c[2] = -(x0*(x[1]+x[2]+x[3])
5439            +x1*(x[0]+x[2]+x[3])
5440            +x2*(x[0]+x[1]+x[3])
5441            +x3*(x[0]+x[1]+x[2]));
5442   c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5443           +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5444           +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5445           +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5446
5447   c[0] = -(x0*x[1]*x[2]*x[3]
5448           +x1*x[0]*x[2]*x[3]
5449           +x2*x[0]*x[1]*x[3]
5450           +x3*x[0]*x[1]*x[2]);
5451
5452
5453   return c;
5454
5455
5456 }
5457
5458 //_____________________________________________________________________________
5459 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
5460 {
5461   //
5462   // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5463   //
5464   Double_t *c = new Double_t[5];
5465   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5466   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5467   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5468   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5469   Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5470
5471
5472   c[4] = x0+x1+x2+x3+x4;
5473   c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
5474            +x1*(x[0]+x[2]+x[3]+x[4])
5475            +x2*(x[0]+x[1]+x[3]+x[4])
5476            +x3*(x[0]+x[1]+x[2]+x[4])
5477            +x4*(x[0]+x[1]+x[2]+x[3]));
5478   c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5479           +x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5480           +x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
5481           +x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
5482           +x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
5483
5484   c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5485           +x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
5486           +x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
5487           +x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
5488           +x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
5489
5490   c[0] = (x0*x[1]*x[2]*x[3]*x[4]
5491           +x1*x[0]*x[2]*x[3]*x[4]
5492           +x2*x[0]*x[1]*x[3]*x[4]
5493           +x3*x[0]*x[1]*x[2]*x[4]
5494           +x4*x[0]*x[1]*x[2]*x[3]);
5495
5496   return c;
5497
5498
5499 }
5500 //_____________________________________________________________________________
5501 void AliTRDCalibraFit::NormierungCharge()
5502 {
5503   //
5504   // Normalisation of the gain factor resulting for the fits
5505   //
5506
5507   // Calcul of the mean of choosen method by fFitChargeNDB
5508   Double_t sum         = 0.0;
5509   //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5510   for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5511     Int_t    total    = 0;
5512     Int_t    detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5513     Float_t *coef     = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5514     //printf("detector %d coef[0] %f\n",detector,coef[0]);
5515     if (GetStack(detector) == 2) {
5516       total = 1728;
5517     }
5518     if (GetStack(detector) != 2) {
5519       total = 2304;
5520     }
5521     for (Int_t j = 0; j < total; j++) {
5522       if (coef[j] >= 0) {
5523         sum += coef[j];
5524       }
5525     }
5526   }
5527
5528   if (sum > 0) {
5529     fScaleFitFactor = fScaleFitFactor / sum;
5530   }
5531   else {
5532     fScaleFitFactor = 1.0;
5533   }
5534
5535   //methode de boeuf mais bon...
5536   Double_t scalefactor = fScaleFitFactor;
5537
5538   if(fDebugLevel > 1){
5539
5540     if ( !fDebugStreamer ) {
5541       //debug stream
5542       TDirectory *backup = gDirectory;
5543       fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5544       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5545     }
5546     (* fDebugStreamer) << "NormierungCharge"<<
5547       "scalefactor="<<scalefactor<<
5548       "\n";
5549     }
5550 }
5551 //_____________________________________________________________________________
5552 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5553 {
5554   //
5555   // Rebin of the 1D histo for the gain calibration if needed.
5556   // you have to choose fRebin, divider of fNumberBinCharge
5557   //
5558
5559  TAxis *xhist  = hist->GetXaxis();
5560  TH1I  *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5561                                         ,xhist->GetBinLowEdge(1)
5562                                         ,xhist->GetBinUpEdge(xhist->GetNbins()));
5563
5564  AliInfo(Form("fRebin: %d",fRebin));
5565  Int_t i = 1;
5566  for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5567    Double_t sum = 0.0;
5568    for (Int_t ji = i; ji < i+fRebin; ji++) {
5569      sum += hist->GetBinContent(ji);
5570    }
5571    sum = sum / fRebin;
5572    rehist->SetBinContent(k,sum);
5573    i += fRebin;
5574  }
5575
5576  return rehist;
5577
5578 }
5579
5580 //_____________________________________________________________________________
5581 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5582 {
5583   //
5584   // Rebin of the 1D histo for the gain calibration if needed
5585   // you have to choose fRebin divider of fNumberBinCharge
5586   //
5587
5588   TAxis *xhist  = hist->GetXaxis();
5589   TH1F  *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5590                                          ,xhist->GetBinLowEdge(1)
5591                                          ,xhist->GetBinUpEdge(xhist->GetNbins()));
5592
5593   AliInfo(Form("fRebin: %d",fRebin));
5594   Int_t i = 1;
5595   for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5596     Double_t sum = 0.0;
5597     for (Int_t ji = i; ji < i+fRebin; ji++) {
5598       sum += hist->GetBinContent(ji);
5599     }
5600     sum = sum/fRebin;
5601     rehist->SetBinContent(k,sum);
5602     i += fRebin;
5603   }
5604
5605   return rehist;
5606
5607 }
5608 //
5609 //____________Some basic geometry function_____________________________________
5610 //
5611
5612 //_____________________________________________________________________________
5613 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5614 {
5615   //
5616   // Reconstruct the plane number from the detector number
5617   //
5618
5619   return ((Int_t) (d % 6));
5620
5621 }
5622
5623 //_____________________________________________________________________________
5624 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5625 {
5626   //
5627   // Reconstruct the stack number from the detector number
5628   //
5629   const Int_t kNlayer = 6;
5630
5631   return ((Int_t) (d % 30) / kNlayer);
5632
5633 }
5634
5635 //_____________________________________________________________________________
5636 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5637 {
5638   //
5639   // Reconstruct the sector number from the detector number
5640   //
5641   Int_t fg = 30;
5642
5643   return ((Int_t) (d / fg));
5644
5645 }
5646
5647 //
5648 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5649 //
5650 //_______________________________________________________________________________
5651 void AliTRDCalibraFit::ResetVectorFit()
5652 {
5653   //
5654   // Reset the VectorFits
5655   //
5656
5657   fVectorFit.SetOwner();
5658   fVectorFit.Clear();
5659   fVectorFit2.SetOwner();
5660   fVectorFit2.Clear();
5661
5662 }
5663 //
5664 //____________Private Functions________________________________________________
5665 //
5666
5667 //_____________________________________________________________________________
5668 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5669 {
5670   //
5671   // Function for the fit
5672   //
5673
5674   //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5675
5676   //PARAMETERS FOR FIT PH
5677   // PASAv.4
5678   //fAsymmGauss->SetParameter(0,0.113755);
5679   //fAsymmGauss->SetParameter(1,0.350706);
5680   //fAsymmGauss->SetParameter(2,0.0604244);
5681   //fAsymmGauss->SetParameter(3,7.65596);
5682   //fAsymmGauss->SetParameter(4,1.00124);
5683   //fAsymmGauss->SetParameter(5,0.870597);  // No tail cancelation
5684
5685   Double_t xx = x[0];
5686
5687   if (xx < par[1]) {
5688     return par[5];
5689   }
5690
5691   Double_t dx       = 0.005;
5692   Double_t xs       = par[1];
5693   Double_t ss       = 0.0;
5694   Double_t paras[2] = { 0.0, 0.0 };
5695
5696   while (xs < xx) {
5697     if ((xs >= par[1]) &&
5698         (xs < (par[1]+par[2]))) {
5699       //fAsymmGauss->SetParameter(0,par[0]);
5700       //fAsymmGauss->SetParameter(1,xs);
5701       //ss += fAsymmGauss->Eval(xx);
5702       paras[0] = par[0];
5703       paras[1] = xs;
5704       ss += AsymmGauss(&xx,paras);
5705     }
5706     if ((xs >= (par[1]+par[2])) &&
5707         (xs <  (par[1]+par[2]+par[3]))) {
5708       //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5709       //fAsymmGauss->SetParameter(1,xs);
5710       //ss += fAsymmGauss->Eval(xx);
5711       paras[0] = par[0]*par[4];
5712       paras[1] = xs;
5713       ss += AsymmGauss(&xx,paras);
5714     }
5715     xs += dx;
5716   }
5717
5718   return ss + par[5];
5719
5720 }
5721
5722 //_____________________________________________________________________________
5723 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5724 {
5725   //
5726   // Function for the fit
5727   //
5728
5729   //par[0] = normalization
5730   //par[1] = mean
5731   //par[2] = sigma
5732   //norm0  = 1
5733   //par[3] = lambda0
5734   //par[4] = norm1
5735   //par[5] = lambda1
5736
5737   Double_t par1save = par[1];
5738   //Double_t par2save = par[2];
5739   Double_t par2save = 0.0604244;
5740   //Double_t par3save = par[3];
5741   Double_t par3save = 7.65596;
5742   //Double_t par5save = par[5];
5743   Double_t par5save = 0.870597;
5744   Double_t dx       = x[0] - par1save;
5745
5746   Double_t  sigma2  = par2save*par2save;
5747   Double_t  sqrt2   = TMath::Sqrt(2.0);
5748   Double_t  exp1    = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5749                                * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5750   Double_t  exp2    = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5751                                * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5752
5753   //return par[0]*(exp1+par[4]*exp2);
5754   return par[0] * (exp1 + 1.00124 * exp2);
5755
5756 }
5757
5758 //_____________________________________________________________________________
5759 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5760 {
5761   //
5762   // Sum Landau + Gaus with identical mean
5763   //
5764
5765   Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5766   //Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5767   Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5768   Double_t val       = valLandau + valGaus;
5769
5770   return val;
5771
5772 }
5773
5774 //_____________________________________________________________________________
5775 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
5776 {
5777   //
5778   // Function for the fit
5779   //
5780   // Fit parameters:
5781   // par[0]=Width (scale) parameter of Landau density
5782   // par[1]=Most Probable (MP, location) parameter of Landau density
5783   // par[2]=Total area (integral -inf to inf, normalization constant)
5784   // par[3]=Width (sigma) of convoluted Gaussian function
5785   //
5786   // In the Landau distribution (represented by the CERNLIB approximation),
5787   // the maximum is located at x=-0.22278298 with the location parameter=0.
5788   // This shift is corrected within this function, so that the actual
5789   // maximum is identical to the MP parameter.
5790   //
5791
5792   // Numeric constants
5793   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
5794   Double_t mpshift  = -0.22278298;       // Landau maximum location
5795
5796   // Control constants
5797   Double_t np       = 100.0;             // Number of convolution steps
5798   Double_t sc       =   5.0;             // Convolution extends to +-sc Gaussian sigmas
5799
5800   // Variables
5801   Double_t xx;
5802   Double_t mpc;
5803   Double_t fland;
5804   Double_t sum = 0.0;
5805   Double_t xlow;
5806   Double_t xupp;
5807   Double_t step;
5808   Double_t i;
5809
5810   // MP shift correction
5811   mpc = par[1] - mpshift * par[0];
5812
5813   // Range of convolution integral
5814   xlow = x[0] - sc * par[3];
5815   xupp = x[0] + sc * par[3];
5816
5817   step = (xupp - xlow) / np;
5818
5819   // Convolution integral of Landau and Gaussian by sum
5820   for (i = 1.0; i <= np/2; i++) {
5821
5822     xx    = xlow + (i-.5) * step;
5823     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5824     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
5825
5826     xx    = xupp - (i-.5) * step;
5827     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
5828     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
5829
5830   }
5831
5832   return (par[2] * step * sum * invsq2pi / par[3]);
5833
5834 }
5835 //_____________________________________________________________________________
5836 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
5837                                       , const Double_t *parlimitslo, const Double_t *parlimitshi
5838                                       , Double_t *fitparams, Double_t *fiterrors
5839                                       , Double_t *chiSqr, Int_t *ndf) const
5840 {
5841   //
5842   // Function for the fit
5843   //
5844
5845   Int_t i;
5846   Char_t funname[100];
5847
5848   TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
5849   if (ffitold) {
5850     delete ffitold;
5851   }
5852
5853   TF1 *ffit    = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
5854   ffit->SetParameters(startvalues);
5855   ffit->SetParNames("Width","MP","Area","GSigma");
5856
5857   for (i = 0; i < 4; i++) {
5858     ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
5859   }
5860
5861   his->Fit(funname,"RB0");                   // Fit within specified range, use ParLimits, do not plot
5862
5863   ffit->GetParameters(fitparams);            // Obtain fit parameters
5864   for (i = 0; i < 4; i++) {
5865     fiterrors[i] = ffit->GetParError(i);     // Obtain fit parameter errors
5866   }
5867   chiSqr[0] = ffit->GetChisquare();          // Obtain chi^2
5868   ndf[0]    = ffit->GetNDF();                // Obtain ndf
5869
5870   return (ffit);                             // Return fit function
5871
5872 }
5873
5874 //_____________________________________________________________________________
5875 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
5876 {
5877   //
5878   // Function for the fit
5879   //
5880
5881   Double_t p;
5882   Double_t x;
5883   Double_t fy;
5884   Double_t fxr;
5885   Double_t fxl;
5886   Double_t step;
5887   Double_t l;
5888   Double_t lold;
5889
5890   Int_t    i        = 0;
5891   Int_t    maxcalls = 10000;
5892
5893   // Search for maximum
5894   p    = params[1] - 0.1 * params[0];
5895   step = 0.05 * params[0];
5896   lold = -2.0;
5897   l    = -1.0;
5898
5899   while ((l != lold) && (i < maxcalls)) {
5900     i++;
5901     lold = l;
5902     x    = p + step;
5903     l    = LanGauFun(&x,params);
5904     if (l < lold) {
5905       step = -step / 10.0;
5906     }
5907     p += step;
5908   }
5909
5910   if (i == maxcalls) {
5911     return (-1);
5912   }
5913   maxx = x;
5914   fy = l / 2.0;
5915
5916   // Search for right x location of fy
5917   p    = maxx + params[0];
5918   step = params[0];
5919   lold = -2.0;
5920   l    = -1e300;
5921   i    = 0;
5922
5923   while ( (l != lold) && (i < maxcalls) ) {
5924     i++;
5925
5926     lold = l;
5927     x = p + step;
5928     l = TMath::Abs(LanGauFun(&x,params) - fy);
5929
5930     if (l > lold)
5931       step = -step/10;
5932
5933     p += step;
5934   }
5935
5936   if (i == maxcalls)
5937     return (-2);
5938
5939   fxr = x;
5940
5941
5942   // Search for left x location of fy
5943
5944   p = maxx - 0.5 * params[0];
5945   step = -params[0];
5946   lold = -2.0;
5947   l    = -1.0e300;
5948   i    = 0;
5949
5950   while ((l != lold) && (i < maxcalls)) {
5951     i++;
5952     lold = l;
5953     x    = p + step;
5954     l    = TMath::Abs(LanGauFun(&x,params) - fy);
5955     if (l > lold) {
5956       step = -step / 10.0;
5957     }
5958     p += step;
5959   }
5960
5961   if (i == maxcalls) {
5962     return (-3);
5963   }
5964
5965   fxl  = x;
5966   fwhm = fxr - fxl;
5967
5968   return (0);
5969 }
5970 //_____________________________________________________________________________
5971 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
5972 {
5973   //
5974   // Gaus with identical mean
5975   //
5976
5977   Double_t gauss   = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
5978
5979   return gauss;
5980
5981 }
5982