TRD/AliTRDCalibraFit.cxx

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