TRD/AliTRDCalibraFit.cxx

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