TRD/AliTRDCalibraFit.cxx

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