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