TRD/AliTRDCalibraFit.cxx

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