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