TRD/AliTRDCalibraFit.cxx

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