TRD/AliTRDCalibraFit.cxx

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