TRD/AliTRDCalibraFit.cxx

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