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     object->SetValue(detector,mean);
1948   }
1949
1950   return object;
1951 }
1952 //_____________________________________________________________________________
1953 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector)
1954 {
1955   //
1956   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
1957   // It takes the min value of the coefficients per detector
1958   // This object has to be written in the database
1959   //
1960
1961   // Create the DetObject
1962   AliTRDCalDet *object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
1963
1964   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
1965   if(loop != 540) AliInfo("The Vector Fit is not complete!");
1966   Int_t detector = -1;
1967   Float_t value  = 0.0;
1968
1969   for (Int_t k = 0; k < loop; k++) {
1970     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
1971     Float_t min  = 100.0;
1972     if(perdetector){
1973       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
1974       // check successful
1975       if(value > 70.0) value = value-100.0;
1976       //
1977       min = value;
1978     }
1979     else{
1980       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
1981       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
1982       for (Int_t row = 0; row < rowMax; row++) {
1983         for (Int_t col = 0; col < colMax; col++) {
1984           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
1985           // check successful
1986           if(value > 70.0) value = value-100.0;
1987           //
1988           if(min > value) min = value;
1989         } // Col
1990       } // Row
1991     }
1992     object->SetValue(detector,min);
1993   }
1994
1995   return object;
1996
1997 }
1998 //_____________________________________________________________________________
1999 AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit)
2000 {
2001   //
2002   // It creates the AliTRDCalDet object from the AliTRDFitInfo2
2003   // It takes the min value of the coefficients per detector
2004   // This object has to be written in the database
2005   //
2006
2007   // Create the DetObject
2008   AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
2009
2010
2011   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2012   if(loop != 540) AliInfo("The Vector Fit is not complete!");
2013   Int_t detector = -1;
2014   Float_t value  = 0.0;
2015
2016   for (Int_t k = 0; k < loop; k++) {
2017     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2018     /*
2019       Int_t rowMax    = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
2020       Int_t colMax    = fGeo->GetColMax(GetLayer(detector));
2021       Float_t min  = 100.0;
2022       for (Int_t row = 0; row < rowMax; row++) {
2023       for (Int_t col = 0; col < colMax; col++) {
2024       value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2025       mean += -TMath::Abs(value);
2026       count++;
2027       } // Col
2028       } // Row
2029       if(count > 0) mean = mean/count;
2030     */
2031     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2032     object->SetValue(detector,-TMath::Abs(value));
2033   }
2034
2035   return object;
2036
2037 }
2038 //_____________________________________________________________________________
2039 TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
2040 {
2041   //
2042   // It Creates the AliTRDCalPad object from AliTRDFitInfo
2043   // You need first to create the object for the detectors,
2044   // where the mean value is put.
2045   // This object has to be written in the database
2046   //
2047
2048   // Create the DetObject
2049   AliTRDCalPad *object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
2050
2051   if(!vectorFit){
2052     for(Int_t k = 0; k < 540; k++){
2053       AliTRDCalROC *calROC = object->GetCalROC(k);
2054       Int_t nchannels = calROC->GetNchannels();
2055       for(Int_t ch = 0; ch < nchannels; ch++){
2056         calROC->SetValue(ch,1.0);
2057       }
2058     }
2059   }
2060   else{
2061
2062     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2063     if(loop != 540) AliInfo("The Vector Fit is not complete!");
2064     Int_t detector = -1;
2065     Float_t value  = 0.0;
2066
2067     for (Int_t k = 0; k < loop; k++) {
2068       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2069       AliTRDCalROC *calROC = object->GetCalROC(detector);
2070       Float_t mean         = detobject->GetValue(detector);
2071       if(TMath::Abs(mean) <= 0.0000000001) continue;
2072       Int_t rowMax    = calROC->GetNrows();
2073       Int_t colMax    = calROC->GetNcols();
2074       for (Int_t row = 0; row < rowMax; row++) {
2075         for (Int_t col = 0; col < colMax; col++) {
2076           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2077           if(value > 0) value = value*scaleFitFactor;
2078           calROC->SetValue(col,row,TMath::Abs(value)/mean);
2079         } // Col
2080       } // Row
2081     }
2082   }
2083
2084   return object;
2085 }
2086 //_____________________________________________________________________________
2087 TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2088 {
2089   //
2090   // It Creates the AliTRDCalPad object from AliTRDFitInfo
2091   // You need first to create the object for the detectors,
2092   // where the mean value is put.
2093   // This object has to be written in the database
2094   //
2095
2096   // Create the DetObject
2097   AliTRDCalPad *object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
2098
2099   if(!vectorFit){
2100     for(Int_t k = 0; k < 540; k++){
2101       AliTRDCalROC *calROC = object->GetCalROC(k);
2102       Int_t nchannels = calROC->GetNchannels();
2103       for(Int_t ch = 0; ch < nchannels; ch++){
2104         calROC->SetValue(ch,1.0);
2105       }
2106     }
2107   }
2108   else {
2109
2110     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2111     if(loop != 540) AliInfo("The Vector Fit is not complete!");
2112     Int_t detector = -1;
2113     Float_t value  = 0.0;
2114
2115     for (Int_t k = 0; k < loop; k++) {
2116       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2117       AliTRDCalROC *calROC = object->GetCalROC(detector);
2118       Float_t mean         = detobject->GetValue(detector);
2119       if(mean == 0) continue;
2120       Int_t rowMax    = calROC->GetNrows();
2121       Int_t colMax    = calROC->GetNcols();
2122       for (Int_t row = 0; row < rowMax; row++) {
2123         for (Int_t col = 0; col < colMax; col++) {
2124           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2125           calROC->SetValue(col,row,TMath::Abs(value)/mean);
2126         } // Col
2127       } // Row
2128     }
2129   }
2130   return object;
2131
2132 }
2133 //_____________________________________________________________________________
2134 TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject)
2135 {
2136   //
2137   // It Creates the AliTRDCalPad object from AliTRDFitInfo2
2138   // You need first to create the object for the detectors,
2139   // where the mean value is put.
2140   // This object has to be written in the database
2141   //
2142
2143   // Create the DetObject
2144   AliTRDCalPad *object = new AliTRDCalPad("LocalT0","T0 (local variations)");
2145
2146   if(!vectorFit){
2147     for(Int_t k = 0; k < 540; k++){
2148       AliTRDCalROC *calROC = object->GetCalROC(k);
2149       Int_t nchannels = calROC->GetNchannels();
2150       for(Int_t ch = 0; ch < nchannels; ch++){
2151         calROC->SetValue(ch,0.0);
2152       }
2153     }
2154   }
2155   else {
2156
2157     Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2158     if(loop != 540) AliInfo("The Vector Fit is not complete!");
2159     Int_t detector = -1;
2160     Float_t value  = 0.0;
2161
2162     for (Int_t k = 0; k < loop; k++) {
2163       detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2164       AliTRDCalROC *calROC = object->GetCalROC(detector);
2165       Float_t min          = detobject->GetValue(detector);
2166       Int_t rowMax    = calROC->GetNrows();
2167       Int_t colMax    = calROC->GetNcols();
2168       for (Int_t row = 0; row < rowMax; row++) {
2169         for (Int_t col = 0; col < colMax; col++) {
2170           value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2171           // check successful
2172           if(value > 70.0) value = value - 100.0;
2173           //
2174           calROC->SetValue(col,row,value-min);
2175         } // Col
2176       } // Row
2177     }
2178   }
2179   return object;
2180
2181 }
2182 //_____________________________________________________________________________
2183 TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit)
2184 {
2185   //
2186   // It Creates the AliTRDCalPad object from AliTRDFitInfo
2187   // This object has to be written in the database
2188   //
2189
2190   // Create the DetObject
2191   AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
2192
2193   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2194   if(loop != 540) AliInfo("The Vector Fit is not complete!");
2195   Int_t detector = -1;
2196   Float_t value  = 0.0;
2197
2198   for (Int_t k = 0; k < loop; k++) {
2199     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2200     AliTRDCalROC *calROC = object->GetCalROC(detector);
2201     Int_t rowMax    = calROC->GetNrows();
2202     Int_t colMax    = calROC->GetNcols();
2203     for (Int_t row = 0; row < rowMax; row++) {
2204       for (Int_t col = 0; col < colMax; col++) {
2205         value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
2206         calROC->SetValue(col,row,TMath::Abs(value));
2207       } // Col
2208     } // Row
2209   }
2210
2211   return object;
2212
2213 }
2214 //_____________________________________________________________________________
2215 AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean)
2216 {
2217   //
2218   // It Creates the AliTRDCalDet object from AliTRDFitInfo
2219   // 0 successful fit 1 not successful fit
2220   // mean is the mean value over the successful fit
2221   // do not use it for t0: no meaning
2222   //
2223
2224   // Create the CalObject
2225   AliTRDCalDet *object = new AliTRDCalDet(name,name);
2226   mean = 0.0;
2227   Int_t count = 0;
2228
2229   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2230   if(loop != 540) {
2231     AliInfo("The Vector Fit is not complete! We initialise all outliers");
2232     for(Int_t k = 0; k < 540; k++){
2233       object->SetValue(k,1.0);
2234     }
2235   }
2236   Int_t detector = -1;
2237   Float_t value  = 0.0;
2238
2239   for (Int_t k = 0; k < loop; k++) {
2240     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2241     value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
2242     if(value <= 0) object->SetValue(detector,1.0);
2243     else {
2244       object->SetValue(detector,0.0);
2245       mean += value;
2246       count++;
2247     }
2248   }
2249   if(count > 0) mean /= count;
2250   return object;
2251 }
2252 //_____________________________________________________________________________
2253 TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean)
2254 {
2255   //
2256   // It Creates the AliTRDCalPad object from AliTRDFitInfo
2257   // 0 not successful fit 1 successful fit
2258   // mean mean value over the successful fit
2259   //
2260
2261   // Create the CalObject
2262   AliTRDCalPad *object = new AliTRDCalPad(name,name);
2263   mean = 0.0;
2264   Int_t count = 0;
2265
2266   Int_t loop = (Int_t) vectorFit->GetEntriesFast();
2267   if(loop != 540) {
2268     AliInfo("The Vector Fit is not complete! We initialise all outliers");
2269     for(Int_t k = 0; k < 540; k++){
2270       AliTRDCalROC *calROC = object->GetCalROC(k);
2271       Int_t nchannels = calROC->GetNchannels();
2272       for(Int_t ch = 0; ch < nchannels; ch++){
2273         calROC->SetValue(ch,1.0);
2274       }
2275     }
2276   }
2277   Int_t detector = -1;
2278   Float_t value  = 0.0;
2279
2280   for (Int_t k = 0; k < loop; k++) {
2281     detector  = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
2282     AliTRDCalROC *calROC = object->GetCalROC(detector);
2283     Int_t nchannels    = calROC->GetNchannels();
2284     for (Int_t ch = 0; ch < nchannels; ch++) {
2285       value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[ch];
2286       if(value <= 0) calROC->SetValue(ch,1.0);
2287       else {
2288         calROC->SetValue(ch,0.0);
2289         mean += value;
2290         count++;
2291       }
2292     } // channels
2293   }
2294   if(count > 0) mean /= count;
2295   return object;
2296 }
2297 //_____________________________________________________________________________
2298 void AliTRDCalibraFit::SetPeriodeFitPH(Int_t periodeFitPH)
2299 {
2300   //
2301   // Set FitPH if 1 then each detector will be fitted
2302   //
2303
2304   if (periodeFitPH > 0) {
2305     fFitPHPeriode   = periodeFitPH;
2306   }
2307   else {
2308     AliInfo("periodeFitPH must be higher than 0!");
2309   }
2310
2311 }
2312 //_____________________________________________________________________________
2313 void AliTRDCalibraFit::SetBeginFitCharge(Float_t beginFitCharge)
2314 {
2315   //
2316   // The fit of the deposited charge distribution begins at
2317   // histo->Mean()/beginFitCharge
2318   // You can here set beginFitCharge
2319   //
2320
2321   if (beginFitCharge > 0) {
2322     fBeginFitCharge = beginFitCharge;
2323   }
2324   else {
2325     AliInfo("beginFitCharge must be strict positif!");
2326   }
2327
2328 }
2329
2330 //_____________________________________________________________________________
2331 void AliTRDCalibraFit::SetT0Shift0(Float_t t0Shift)
2332 {
2333   //
2334   // The t0 calculated with the maximum positif slope is shift from t0Shift0
2335   // You can here set t0Shift0
2336   //
2337
2338   if (t0Shift > 0) {
2339     fT0Shift0 = t0Shift;
2340   }
2341   else {
2342     AliInfo("t0Shift0 must be strict positif!");
2343   }
2344
2345 }
2346
2347 //_____________________________________________________________________________
2348 void AliTRDCalibraFit::SetT0Shift1(Float_t t0Shift)
2349 {
2350   //
2351   // The t0 calculated with the maximum of the amplification region is shift from t0Shift1
2352   // You can here set t0Shift1
2353   //
2354
2355   if (t0Shift > 0) {
2356     fT0Shift1 = t0Shift;
2357   }
2358   else {
2359     AliInfo("t0Shift must be strict positif!");
2360   }
2361
2362 }
2363
2364 //_____________________________________________________________________________
2365 void AliTRDCalibraFit::SetRangeFitPRF(Float_t rangeFitPRF)
2366 {
2367   //
2368   // The fit of the PRF is from -rangeFitPRF to rangeFitPRF
2369   // You can here set rangeFitPRF
2370   //
2371
2372   if ((rangeFitPRF >    0) &&
2373       (rangeFitPRF <= 1.5)) {
2374     fRangeFitPRF = rangeFitPRF;
2375   }
2376   else {
2377     AliInfo("rangeFitPRF must be between 0 and 1.0");
2378   }
2379
2380 }
2381
2382 //_____________________________________________________________________________
2383 void AliTRDCalibraFit::SetMinEntries(Int_t minEntries)
2384 {
2385   //
2386   // Minimum entries for fitting
2387   //
2388
2389   if (minEntries >    0) {
2390     fMinEntries = minEntries;
2391   }
2392   else {
2393     AliInfo("fMinEntries must be >= 0.");
2394   }
2395
2396 }
2397
2398 //_____________________________________________________________________________
2399 void AliTRDCalibraFit::SetRebin(Short_t rebin)
2400 {
2401   //
2402   // Rebin with rebin time less bins the Ch histo
2403   // You can set here rebin that should divide the number of bins of CH histo
2404   //
2405
2406   if (rebin > 0) {
2407     fRebin = rebin;
2408     AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
2409   }
2410   else {
2411     AliInfo("You have to choose a positiv value!");
2412   }
2413
2414 }
2415 //_____________________________________________________________________________
2416 Bool_t AliTRDCalibraFit::FillVectorFit()
2417 {
2418   //
2419   // For the Fit functions fill the vector Fit
2420   //
2421
2422   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2423
2424   Int_t ntotal = 1;
2425   if (GetStack(fCountDet) == 2) {
2426     ntotal = 1728;
2427   }
2428   else {
2429     ntotal = 2304;
2430   }
2431
2432   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2433   Float_t *coef = new Float_t[ntotal];
2434   for (Int_t i = 0; i < ntotal; i++) {
2435     coef[i] = fCurrentCoefDetector[i];
2436   }
2437
2438   Int_t detector = fCountDet;
2439   // Set
2440   fitInfo->SetCoef(coef);
2441   fitInfo->SetDetector(detector);
2442   fVectorFit.Add((TObject *) fitInfo);
2443
2444   return kTRUE;
2445
2446 }
2447 //_____________________________________________________________________________
2448 Bool_t AliTRDCalibraFit::FillVectorFit2()
2449 {
2450   //
2451   // For the Fit functions fill the vector Fit
2452   //
2453
2454   AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
2455
2456   Int_t ntotal = 1;
2457   if (GetStack(fCountDet) == 2) {
2458     ntotal = 1728;
2459   }
2460   else {
2461     ntotal = 2304;
2462   }
2463
2464   //printf("For the detector %d , ntotal %d and fCoefCH[0] %f\n",countdet,ntotal,fCoefCH[0]);
2465   Float_t *coef = new Float_t[ntotal];
2466   for (Int_t i = 0; i < ntotal; i++) {
2467     coef[i] = fCurrentCoefDetector2[i];
2468   }
2469
2470   Int_t detector = fCountDet;
2471   // Set
2472   fitInfo->SetCoef(coef);
2473   fitInfo->SetDetector(detector);
2474   fVectorFit2.Add((TObject *) fitInfo);
2475
2476   return kTRUE;
2477
2478 }
2479 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2480 Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i)
2481 {
2482   //
2483   // Init the number of expected bins and fDect1[i] fDect2[i]
2484   //
2485
2486   gStyle->SetPalette(1);
2487   gStyle->SetOptStat(1111);
2488   gStyle->SetPadBorderMode(0);
2489   gStyle->SetCanvasColor(10);
2490   gStyle->SetPadLeftMargin(0.13);
2491   gStyle->SetPadRightMargin(0.01);
2492
2493   // Mode groups of pads: the total number of bins!
2494   CalculNumberOfBinsExpected(i);
2495
2496   // Quick verification that we have the good pad calibration mode!
2497   if (fNumberOfBinsExpected != nbins) {
2498     AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins));
2499     return kFALSE;
2500   }
2501
2502   // Security for fDebug 3 and 4
2503   if ((fDebugLevel >= 3) &&
2504       ((fDet[0] >  5) ||
2505        (fDet[1] >  4) ||
2506        (fDet[2] > 17))) {
2507     AliInfo("This detector doesn't exit!");
2508     return kFALSE;
2509   }
2510
2511   // Determine fDet1 and fDet2 and set the fNfragZ and fNfragRphi for debug 3 and 4
2512   CalculDect1Dect2(i);
2513
2514
2515   return kTRUE;
2516 }
2517 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2518 Bool_t AliTRDCalibraFit::InitFitCH()
2519 {
2520   //
2521   // Init the fVectorFitCH for normalisation
2522   // Init the histo for debugging
2523   //
2524
2525   gDirectory = gROOT;
2526
2527   fScaleFitFactor = 0.0;
2528   fCurrentCoefDetector   = new Float_t[2304];
2529   for (Int_t k = 0; k < 2304; k++) {
2530     fCurrentCoefDetector[k] = 0.0;
2531   }
2532   fVectorFit.SetName("gainfactorscoefficients");
2533
2534   // fDebug == 0 nothing
2535   // fDebug == 1 and fFitVoir no histo
2536   if (fDebugLevel == 1) {
2537     if(!CheckFitVoir()) return kFALSE;
2538   }
2539   //Get the CalDet object
2540   if(fAccCDB){
2541     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2542     if (!cal) {
2543       AliInfo("Could not get calibDB");
2544       return kFALSE;
2545     }
2546     if(fCalDet) delete fCalDet;
2547     fCalDet = new AliTRDCalDet(*(cal->GetGainFactorDet()));
2548   }
2549   else{
2550     Float_t devalue = 1.0;
2551     if(fCalDet) delete fCalDet;
2552     fCalDet = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
2553     for(Int_t k = 0; k < 540; k++){
2554       fCalDet->SetValue(k,devalue);
2555     }
2556   }
2557   return kTRUE;
2558
2559 }
2560 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2561 Bool_t AliTRDCalibraFit::InitFitPH()
2562 {
2563   //
2564   // Init the arrays of results
2565   // Init the histos for debugging
2566   //
2567
2568   gDirectory = gROOT;
2569   fVectorFit.SetName("driftvelocitycoefficients");
2570   fVectorFit2.SetName("t0coefficients");
2571
2572   fCurrentCoefDetector   = new Float_t[2304];
2573   for (Int_t k = 0; k < 2304; k++) {
2574     fCurrentCoefDetector[k] = 0.0;
2575   }
2576
2577   fCurrentCoefDetector2   = new Float_t[2304];
2578   for (Int_t k = 0; k < 2304; k++) {
2579     fCurrentCoefDetector2[k] = 0.0;
2580   }
2581
2582   //fDebug == 0 nothing
2583   // fDebug == 1 and fFitVoir no histo
2584   if (fDebugLevel == 1) {
2585     if(!CheckFitVoir()) return kFALSE;
2586   }
2587   //Get the CalDet object
2588   if(fAccCDB){
2589     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2590     if (!cal) {
2591       AliInfo("Could not get calibDB");
2592       return kFALSE;
2593     }
2594     if(fCalDet) delete fCalDet;
2595     if(fCalDet2) delete fCalDet2;
2596     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2597     fCalDet2 = new AliTRDCalDet(*(cal->GetT0Det()));
2598   }
2599   else{
2600     Float_t devalue  = 1.5;
2601     Float_t devalue2 = 0.0;
2602     if(fCalDet) delete fCalDet;
2603     if(fCalDet2) delete fCalDet2;
2604     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2605     fCalDet2 = new AliTRDCalDet("ChamberT0","T0 (detector value)");
2606     for(Int_t k = 0; k < 540; k++){
2607       fCalDet->SetValue(k,devalue);
2608       fCalDet2->SetValue(k,devalue2);
2609     }
2610   }
2611   return kTRUE;
2612 }
2613 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2614 Bool_t AliTRDCalibraFit::InitFitPRF()
2615 {
2616   //
2617   // Init the calibration mode (Nz, Nrphi), the histograms for
2618   // debugging the fit methods if fDebug > 0,
2619   //
2620
2621   gDirectory = gROOT;
2622   fVectorFit.SetName("prfwidthcoefficients");
2623
2624   fCurrentCoefDetector   = new Float_t[2304];
2625   for (Int_t k = 0; k < 2304; k++) {
2626     fCurrentCoefDetector[k] = 0.0;
2627   }
2628
2629   // fDebug == 0 nothing
2630   // fDebug == 1 and fFitVoir no histo
2631   if (fDebugLevel == 1) {
2632     if(!CheckFitVoir()) return kFALSE;
2633   }
2634   return kTRUE;
2635 }
2636 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2637 Bool_t AliTRDCalibraFit::InitFitLinearFitter()
2638 {
2639   //
2640   // Init the fCalDet, fVectorFit fCurrentCoefDetector
2641   //
2642
2643   gDirectory = gROOT;
2644
2645   fCurrentCoefDetector   = new Float_t[2304];
2646   fCurrentCoefDetector2  = new Float_t[2304];
2647   for (Int_t k = 0; k < 2304; k++) {
2648     fCurrentCoefDetector[k]  = 0.0;
2649     fCurrentCoefDetector2[k] = 0.0;
2650   }
2651
2652   //printf("test0\n");
2653
2654   AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
2655   if (!cal) {
2656     AliInfo("Could not get calibDB");
2657     return kFALSE;
2658   }
2659
2660   //Get the CalDet object
2661   if(fAccCDB){
2662     if(fCalDet) delete fCalDet;
2663     if(fCalDet2) delete fCalDet2;
2664     fCalDet  = new AliTRDCalDet(*(cal->GetVdriftDet()));
2665     //printf("test1\n");
2666     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2667     //printf("test2\n");
2668     for(Int_t k = 0; k < 540; k++){
2669       fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
2670     }
2671     //printf("test3\n");
2672   }
2673   else{
2674     Float_t devalue  = 1.5;
2675     Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
2676     if(fCalDet) delete fCalDet;
2677     if(fCalDet2) delete fCalDet2;
2678     //printf("test1\n");
2679     fCalDet  = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
2680     fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
2681     //printf("test2\n");
2682     for(Int_t k = 0; k < 540; k++){
2683       fCalDet->SetValue(k,devalue);
2684       fCalDet2->SetValue(k,devalue2);
2685     }
2686     //printf("test3\n");
2687   }
2688   return kTRUE;
2689 }
2690
2691 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2692 void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
2693 {
2694   //
2695   // Init the current detector where we are fCountDet and the
2696   // next fCount for the functions Fit...
2697   //
2698
2699   // Loop on the Xbins of ch!!
2700   fCountDet = -1; // Current detector
2701   fCount    =  0; // To find the next detector
2702
2703   // If fDebug >= 3
2704   if (fDebugLevel >= 3) {
2705     // Set countdet to the detector
2706     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2707     // Set counter to write at the end of the detector
2708     fCount = fDect2;
2709     // Get the right calib objects
2710     SetCalROC(i);
2711   }
2712   if(fDebugLevel == 1) {
2713     fCountDet = 0;
2714     fCalibraMode->CalculXBins(fCountDet,i);
2715     if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
2716       while(fCalibraMode->GetXbins(i) <=fFitVoir){
2717         fCountDet++;
2718         fCalibraMode->CalculXBins(fCountDet,i);
2719         //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
2720       }
2721     }
2722     else {
2723       fCountDet++;
2724     }
2725     fCount    = fCalibraMode->GetXbins(i);
2726     fCountDet--;
2727     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2728     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2729     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2730                                       ,(Int_t) GetStack(fCountDet)
2731                                       ,(Int_t) GetSector(fCountDet),i);
2732   }
2733 }
2734 //_______________________________________________________________________________
2735 void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i)
2736 {
2737   //
2738   // Calculate the number of bins expected (calibration groups)
2739   //
2740
2741   fNumberOfBinsExpected = 0;
2742   // All
2743   if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){
2744     fNumberOfBinsExpected = 1;
2745     return;
2746   }
2747   // Per supermodule
2748   if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){
2749     fNumberOfBinsExpected = 18;
2750     return;
2751   }
2752   // More
2753   fCalibraMode->ModePadCalibration(2,i);
2754   fCalibraMode->ModePadFragmentation(0,2,0,i);
2755   fCalibraMode->SetDetChamb2(i);
2756   if (fDebugLevel > 1) {
2757     AliInfo(Form("For the chamber 2: %d",fCalibraMode->GetDetChamb2(i)));
2758   }
2759   fNumberOfBinsExpected += 6 * 18 * fCalibraMode->GetDetChamb2(i);
2760   fCalibraMode->ModePadCalibration(0,i);
2761   fCalibraMode->ModePadFragmentation(0,0,0,i);
2762   fCalibraMode->SetDetChamb0(i);
2763   if (fDebugLevel > 1) {
2764     AliInfo(Form("For the other chamber 0: %d",fCalibraMode->GetDetChamb0(i)));
2765   }
2766   fNumberOfBinsExpected += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
2767
2768 }
2769 //_______________________________________________________________________________
2770 void AliTRDCalibraFit::CalculDect1Dect2(Int_t i)
2771 {
2772   //
2773   // Calculate the range of fits
2774   //
2775
2776   fDect1 = -1;
2777   fDect2 = -1;
2778   if (fDebugLevel == 1) {
2779     fDect1 = fFitVoir;
2780     fDect2 = fDect1 +1;
2781   }
2782   if ((fDebugLevel == 2) || (fDebugLevel == 0)) {
2783     fDect1 = 0;
2784     fDect2 = fNumberOfBinsExpected;
2785   }
2786   if (fDebugLevel >= 3) {
2787     fCountDet = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
2788     fCalibraMode->CalculXBins(fCountDet,i);
2789     fDect1 = fCalibraMode->GetXbins(i);
2790     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2791     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2792     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2793                                       ,(Int_t) GetStack(fCountDet)
2794                                       ,(Int_t) GetSector(fCountDet),i);
2795     // Set for the next detector
2796     fDect2 = fDect1 + fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2797   }
2798 }
2799 //_______________________________________________________________________________
2800 Bool_t AliTRDCalibraFit::CheckFitVoir()
2801 {
2802   //
2803   // Check if fFitVoir is in the range
2804   //
2805
2806   if (fFitVoir < fNumberOfBinsExpected) {
2807     AliInfo(Form("We will see the fit of the object %d",fFitVoir));
2808   }
2809   else {
2810     AliInfo("fFitVoir is out of range of the histo!");
2811     return kFALSE;
2812   }
2813   return kTRUE;
2814 }
2815 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2816 void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
2817 {
2818   //
2819   // See if we are in a new detector and update the
2820   // variables fNfragZ and fNfragRphi if yes
2821   // Will never happen for only one detector (3 and 4)
2822   // Doesn't matter for 2
2823   //
2824   if (fCount == idect) {
2825     // On en est au detector (or first detector in the group)
2826     fCountDet += 1;
2827     AliDebug(2,Form("We are at the detector %d\n",fCountDet));
2828     // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2829     fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i);
2830     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2831                                        ,(Int_t) GetStack(fCountDet)
2832                                        ,(Int_t) GetSector(fCountDet),i);
2833     // Set for the next detector
2834     fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i);
2835     // calib objects
2836     SetCalROC(i);
2837   }
2838 }
2839 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2840 void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
2841 {
2842   //
2843   // Reconstruct the min pad row, max pad row, min pad col and
2844   // max pad col of the calibration group for the Fit functions
2845   // idect is the calibration group inside the detector
2846   //
2847   if (fDebugLevel !=  1) {
2848     fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i);
2849   }
2850   AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))));
2851   AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)));
2852 }
2853 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2854 Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect)
2855 {
2856   //
2857   // For the case where there are not enough entries in the histograms
2858   // of the calibration group, the value present in the choosen database
2859   // will be put. A negativ sign enables to know that a fit was not possible.
2860   //
2861
2862   if (fDebugLevel == 1) {
2863     AliInfo("The element has not enough statistic to be fitted");
2864   }
2865   else if (fNbDet > 0){
2866     Int_t firstdetector = fCountDet;
2867     Int_t lastdetector  = fCountDet+fNbDet;
2868     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2869                  ,idect,firstdetector,lastdetector));
2870     // loop over detectors
2871     for(Int_t det = firstdetector; det < lastdetector; det++){
2872
2873       //Set the calibration object again
2874       fCountDet = det;
2875       SetCalROC(0);
2876
2877       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
2878       // Put them at 1
2879       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
2880       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
2881                                          ,(Int_t) GetStack(fCountDet)
2882                                          ,(Int_t) GetSector(fCountDet),0);
2883       // Reconstruct row min row max
2884       ReconstructFitRowMinRowMax(idect,0);
2885
2886       // Calcul the coef from the database choosen for the detector
2887       CalculChargeCoefMean(kFALSE);
2888
2889       //stack 2, not stack 2
2890       Int_t factor = 0;
2891       if(GetStack(fCountDet) == 2) factor = 12;
2892       else factor = 16;
2893
2894       // Fill the fCurrentCoefDetector with negative value to say: not fitted
2895       for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2896         for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2897           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2898         }
2899       }
2900
2901       //Put default value negative
2902       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2903       fCurrentCoefE   = 0.0;
2904
2905       // Fill the stuff
2906       FillVectorFit();
2907       // Debug
2908       if(fDebugLevel > 1){
2909
2910         if ( !fDebugStreamer ) {
2911           //debug stream
2912           TDirectory *backup = gDirectory;
2913           fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
2914           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2915         }
2916
2917         Int_t   detector   = fCountDet;
2918         Int_t   caligroup  = idect;
2919         Short_t rowmin     = fCalibraMode->GetRowMin(0);
2920         Short_t rowmax     = fCalibraMode->GetRowMax(0);
2921         Short_t colmin     = fCalibraMode->GetColMin(0);
2922         Short_t colmax     = fCalibraMode->GetColMax(0);
2923         Float_t gf         = fCurrentCoef[0];
2924         Float_t gfs        = fCurrentCoef[1];
2925         Float_t gfE        = fCurrentCoefE;
2926
2927         (*fDebugStreamer) << "FillFillCH" <<
2928           "detector=" << detector <<
2929           "caligroup=" << caligroup <<
2930           "rowmin=" << rowmin <<
2931           "rowmax=" << rowmax <<
2932           "colmin=" << colmin <<
2933           "colmax=" << colmax <<
2934           "gf=" << gf <<
2935           "gfs=" << gfs <<
2936           "gfE=" << gfE <<
2937           "\n";
2938
2939       }
2940       // Reset
2941       for (Int_t k = 0; k < 2304; k++) {
2942         fCurrentCoefDetector[k] = 0.0;
2943       }
2944
2945     }// loop detector
2946     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
2947   }
2948   else {
2949
2950     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
2951                  ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
2952
2953     // Calcul the coef from the database choosen
2954     CalculChargeCoefMean(kFALSE);
2955
2956     //stack 2, not stack 2
2957     Int_t factor = 0;
2958     if(GetStack(fCountDet) == 2) factor = 12;
2959     else factor = 16;
2960
2961     // Fill the fCurrentCoefDetector with negative value to say: not fitted
2962     for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
2963       for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
2964         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
2965       }
2966     }
2967
2968     //Put default value negative
2969     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
2970     fCurrentCoefE   = 0.0;
2971
2972     FillFillCH(idect);
2973   }
2974
2975   return kTRUE;
2976 }
2977
2978
2979 //____________Functions for initialising the AliTRDCalibraFit in the code_________
2980 Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries)
2981 {
2982   //
2983   // For the case where there are not enough entries in the histograms
2984   // of the calibration group, the value present in the choosen database
2985   // will be put. A negativ sign enables to know that a fit was not possible.
2986   //
2987   if (fDebugLevel == 1) {
2988     AliInfo("The element has not enough statistic to be fitted");
2989   }
2990   else if (fNbDet > 0) {
2991
2992     Int_t firstdetector = fCountDet;
2993     Int_t lastdetector  = fCountDet+fNbDet;
2994     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
2995                  ,idect,firstdetector,lastdetector));
2996     // loop over detectors
2997     for(Int_t det = firstdetector; det < lastdetector; det++){
2998
2999       //Set the calibration object again
3000       fCountDet = det;
3001       SetCalROC(1);
3002
3003       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3004       // Put them at 1
3005       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3006       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3007                                          ,(Int_t) GetStack(fCountDet)
3008                                          ,(Int_t) GetSector(fCountDet),1);
3009       // Reconstruct row min row max
3010       ReconstructFitRowMinRowMax(idect,1);
3011
3012       // Calcul the coef from the database choosen for the detector
3013       CalculVdriftCoefMean();
3014       CalculT0CoefMean();
3015
3016       //stack 2, not stack 2
3017       Int_t factor = 0;
3018       if(GetStack(fCountDet) == 2) factor = 12;
3019       else factor = 16;
3020
3021       // Fill the fCurrentCoefDetector with negative value to say: not fitted
3022       for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3023         for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3024           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3025           fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3026         }
3027       }
3028
3029       //Put default value negative
3030       fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3031       fCurrentCoefE    = 0.0;
3032       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3033       fCurrentCoefE2   = 0.0;
3034
3035       // Fill the stuff
3036       FillVectorFit();
3037       FillVectorFit2();
3038       // Debug
3039       if(fDebugLevel > 1){
3040
3041         if ( !fDebugStreamer ) {
3042           //debug stream
3043           TDirectory *backup = gDirectory;
3044           fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3045           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3046         }
3047
3048
3049         Int_t   detector     = fCountDet;
3050         Int_t   caligroup    = idect;
3051         Short_t rowmin       = fCalibraMode->GetRowMin(1);
3052         Short_t rowmax       = fCalibraMode->GetRowMax(1);
3053         Short_t colmin       = fCalibraMode->GetColMin(1);
3054         Short_t colmax       = fCalibraMode->GetColMax(1);
3055         Float_t vf           = fCurrentCoef[0];
3056         Float_t vs           = fCurrentCoef[1];
3057         Float_t vfE          = fCurrentCoefE;
3058         Float_t t0f          = fCurrentCoef2[0];
3059         Float_t t0s          = fCurrentCoef2[1];
3060         Float_t t0E          = fCurrentCoefE2;
3061
3062
3063
3064         (* fDebugStreamer) << "FillFillPH"<<
3065         "detector="<<detector<<
3066           "nentries="<<nentries<<
3067           "caligroup="<<caligroup<<
3068           "rowmin="<<rowmin<<
3069           "rowmax="<<rowmax<<
3070           "colmin="<<colmin<<
3071           "colmax="<<colmax<<
3072           "vf="<<vf<<
3073           "vs="<<vs<<
3074           "vfE="<<vfE<<
3075           "t0f="<<t0f<<
3076           "t0s="<<t0s<<
3077           "t0E="<<t0E<<
3078           "\n";
3079       }
3080       // Reset
3081       for (Int_t k = 0; k < 2304; k++) {
3082         fCurrentCoefDetector[k] = 0.0;
3083         fCurrentCoefDetector2[k] = 0.0;
3084       }
3085
3086     }// loop detector
3087     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3088   }
3089   else {
3090
3091     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3092                  ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
3093
3094     CalculVdriftCoefMean();
3095     CalculT0CoefMean();
3096
3097     //stack 2 and not stack 2
3098     Int_t factor = 0;
3099     if(GetStack(fCountDet) == 2) factor = 12;
3100     else factor = 16;
3101
3102
3103     // Fill the fCurrentCoefDetector 2
3104     for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3105       for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3106         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3107         fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0;
3108       }
3109     }
3110
3111     // Put the default value
3112     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3113     fCurrentCoefE    = 0.0;
3114     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
3115     fCurrentCoefE2   = 0.0;
3116
3117     FillFillPH(idect,nentries);
3118
3119   }
3120
3121   return kTRUE;
3122
3123 }
3124
3125
3126 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3127 Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect)
3128 {
3129   //
3130   // For the case where there are not enough entries in the histograms
3131   // of the calibration group, the value present in the choosen database
3132   // will be put. A negativ sign enables to know that a fit was not possible.
3133   //
3134
3135   if (fDebugLevel == 1) {
3136     AliInfo("The element has not enough statistic to be fitted");
3137   }
3138   else if (fNbDet > 0){
3139
3140     Int_t firstdetector = fCountDet;
3141     Int_t lastdetector  = fCountDet+fNbDet;
3142     AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
3143                  ,idect,firstdetector,lastdetector));
3144
3145     // loop over detectors
3146     for(Int_t det = firstdetector; det < lastdetector; det++){
3147
3148       //Set the calibration object again
3149       fCountDet = det;
3150       SetCalROC(2);
3151
3152       // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3153       // Put them at 1
3154       fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3155       fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3156                                          ,(Int_t) GetStack(fCountDet)
3157                                          ,(Int_t) GetSector(fCountDet),2);
3158       // Reconstruct row min row max
3159       ReconstructFitRowMinRowMax(idect,2);
3160
3161       // Calcul the coef from the database choosen for the detector
3162       CalculPRFCoefMean();
3163
3164       //stack 2, not stack 2
3165       Int_t factor = 0;
3166       if(GetStack(fCountDet) == 2) factor = 12;
3167       else factor = 16;
3168
3169       // Fill the fCurrentCoefDetector with negative value to say: not fitted
3170       for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3171         for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3172           fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3173         }
3174       }
3175
3176       //Put default value negative
3177       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3178       fCurrentCoefE   = 0.0;
3179
3180       // Fill the stuff
3181       FillVectorFit();
3182       // Debug
3183       if(fDebugLevel > 1){
3184
3185         if ( !fDebugStreamer ) {
3186           //debug stream
3187           TDirectory *backup = gDirectory;
3188           fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3189           if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3190         }
3191
3192         Int_t   detector     = fCountDet;
3193         Int_t   layer        = GetLayer(fCountDet);
3194         Int_t   caligroup    = idect;
3195         Short_t rowmin       = fCalibraMode->GetRowMin(2);
3196         Short_t rowmax       = fCalibraMode->GetRowMax(2);
3197         Short_t colmin       = fCalibraMode->GetColMin(2);
3198         Short_t colmax       = fCalibraMode->GetColMax(2);
3199         Float_t widf         = fCurrentCoef[0];
3200         Float_t wids         = fCurrentCoef[1];
3201         Float_t widfE        = fCurrentCoefE;
3202
3203         (* fDebugStreamer) << "FillFillPRF"<<
3204           "detector="<<detector<<
3205           "layer="<<layer<<
3206           "caligroup="<<caligroup<<
3207           "rowmin="<<rowmin<<
3208           "rowmax="<<rowmax<<
3209           "colmin="<<colmin<<
3210           "colmax="<<colmax<<
3211           "widf="<<widf<<
3212           "wids="<<wids<<
3213           "widfE="<<widfE<<
3214           "\n";
3215       }
3216       // Reset
3217       for (Int_t k = 0; k < 2304; k++) {
3218         fCurrentCoefDetector[k] = 0.0;
3219       }
3220
3221     }// loop detector
3222     AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet));
3223   }
3224   else {
3225
3226     AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
3227                  ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
3228
3229     CalculPRFCoefMean();
3230
3231     // stack 2 and not stack 2
3232     Int_t factor = 0;
3233     if(GetStack(fCountDet) == 2) factor = 12;
3234     else factor = 16;
3235
3236
3237     // Fill the fCurrentCoefDetector
3238     for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3239       for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3240         fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]);
3241       }
3242     }
3243
3244     // Put the default value
3245     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
3246     fCurrentCoefE   = 0.0;
3247
3248     FillFillPRF(idect);
3249   }
3250
3251   return kTRUE;
3252
3253 }
3254 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3255 Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter()
3256 {
3257   //
3258   // For the case where there are not enough entries in the histograms
3259   // of the calibration group, the value present in the choosen database
3260   // will be put. A negativ sign enables to know that a fit was not possible.
3261   //
3262
3263   // Calcul the coef from the database choosen
3264   CalculVdriftLorentzCoef();
3265
3266   Int_t factor = 0;
3267   if(GetStack(fCountDet) == 2) factor = 1728;
3268   else factor = 2304;
3269
3270
3271   // Fill the fCurrentCoefDetector
3272   for (Int_t k = 0; k < factor; k++) {
3273     fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
3274     // should be negative
3275     fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
3276   }
3277
3278
3279   //Put default opposite sign
3280   fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
3281   fCurrentCoefE    = 0.0;
3282   fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
3283   fCurrentCoefE2 = 0.0;
3284
3285   FillFillLinearFitter();
3286
3287   return kTRUE;
3288 }
3289
3290 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3291 Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
3292 {
3293   //
3294   // Fill the coefficients found with the fits or other
3295   // methods from the Fit functions
3296   //
3297
3298   if (fDebugLevel != 1) {
3299     if (fNbDet > 0){
3300       Int_t firstdetector = fCountDet;
3301       Int_t lastdetector  = fCountDet+fNbDet;
3302       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3303                    ,idect,firstdetector,lastdetector));
3304       // loop over detectors
3305       for(Int_t det = firstdetector; det < lastdetector; det++){
3306
3307         //Set the calibration object again
3308         fCountDet = det;
3309         SetCalROC(0);
3310
3311         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3312         // Put them at 1
3313         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0);
3314         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3315                                            ,(Int_t) GetStack(fCountDet)
3316                                            ,(Int_t) GetSector(fCountDet),0);
3317         // Reconstruct row min row max
3318         ReconstructFitRowMinRowMax(idect,0);
3319
3320         // Calcul the coef from the database choosen for the detector
3321         if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE);
3322         else CalculChargeCoefMean(kTRUE);
3323
3324         //stack 2, not stack 2
3325         Int_t factor = 0;
3326         if(GetStack(fCountDet) == 2) factor = 12;
3327         else factor = 16;
3328
3329         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3330         Double_t coeftoput = 1.0;
3331         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3332         else coeftoput = fCurrentCoef[0];
3333         for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3334           for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3335             fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3336           }
3337         }
3338
3339         // Fill the stuff
3340         FillVectorFit();
3341         // Debug
3342         if(fDebugLevel > 1){
3343
3344           if ( !fDebugStreamer ) {
3345             //debug stream
3346             TDirectory *backup = gDirectory;
3347             fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3348             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3349           }
3350
3351           Int_t   detector   = fCountDet;
3352           Int_t   caligroup  = idect;
3353           Short_t rowmin     = fCalibraMode->GetRowMin(0);
3354           Short_t rowmax     = fCalibraMode->GetRowMax(0);
3355           Short_t colmin     = fCalibraMode->GetColMin(0);
3356           Short_t colmax     = fCalibraMode->GetColMax(0);
3357           Float_t gf         = fCurrentCoef[0];
3358           Float_t gfs        = fCurrentCoef[1];
3359           Float_t gfE        = fCurrentCoefE;
3360
3361           (*fDebugStreamer) << "FillFillCH" <<
3362             "detector=" << detector <<
3363             "caligroup=" << caligroup <<
3364             "rowmin=" << rowmin <<
3365             "rowmax=" << rowmax <<
3366             "colmin=" << colmin <<
3367             "colmax=" << colmax <<
3368             "gf=" << gf <<
3369             "gfs=" << gfs <<
3370             "gfE=" << gfE <<
3371             "\n";
3372
3373         }
3374         // Reset
3375         for (Int_t k = 0; k < 2304; k++) {
3376           fCurrentCoefDetector[k] = 0.0;
3377         }
3378
3379       }// loop detector
3380       //printf("Check the count now: fCountDet %d\n",fCountDet);
3381     }
3382     else{
3383
3384       Int_t factor = 0;
3385       if(GetStack(fCountDet) == 2) factor = 12;
3386       else factor = 16;
3387
3388       for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) {
3389         for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) {
3390           fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3391         }
3392       }
3393
3394       FillFillCH(idect);
3395     }
3396   }
3397
3398   return kTRUE;
3399
3400 }
3401 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3402 Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries)
3403 {
3404   //
3405   // Fill the coefficients found with the fits or other
3406   // methods from the Fit functions
3407   //
3408
3409   if (fDebugLevel != 1) {
3410     if (fNbDet > 0){
3411
3412       Int_t firstdetector = fCountDet;
3413       Int_t lastdetector  = fCountDet+fNbDet;
3414       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3415                    ,idect,firstdetector,lastdetector));
3416
3417       // loop over detectors
3418       for(Int_t det = firstdetector; det < lastdetector; det++){
3419
3420         //Set the calibration object again
3421         fCountDet = det;
3422         SetCalROC(1);
3423
3424         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3425         // Put them at 1
3426         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1);
3427         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3428                                            ,(Int_t) GetStack(fCountDet)
3429                                            ,(Int_t) GetSector(fCountDet),1);
3430         // Reconstruct row min row max
3431         ReconstructFitRowMinRowMax(idect,1);
3432
3433         // Calcul the coef from the database choosen for the detector
3434         CalculVdriftCoefMean();
3435         CalculT0CoefMean();
3436
3437         //stack 2, not stack 2
3438         Int_t factor = 0;
3439         if(GetStack(fCountDet) == 2) factor = 12;
3440         else factor = 16;
3441
3442         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3443         Double_t coeftoput  = 1.5;
3444         Double_t coeftoput2 = 0.0;
3445
3446         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3447         else coeftoput = fCurrentCoef[0];
3448
3449         if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0;
3450         else coeftoput2 = fCurrentCoef2[0];
3451
3452         for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3453           for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3454             fCurrentCoefDetector[(Int_t)(j*factor+k)]  = coeftoput;
3455             fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2;
3456           }
3457         }
3458
3459         // Fill the stuff
3460         FillVectorFit();
3461         FillVectorFit2();
3462         // Debug
3463         if(fDebugLevel > 1){
3464
3465           if ( !fDebugStreamer ) {
3466             //debug stream
3467             TDirectory *backup = gDirectory;
3468             fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3469             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3470           }
3471
3472
3473           Int_t   detector     = fCountDet;
3474           Int_t   caligroup    = idect;
3475           Short_t rowmin       = fCalibraMode->GetRowMin(1);
3476           Short_t rowmax       = fCalibraMode->GetRowMax(1);
3477           Short_t colmin       = fCalibraMode->GetColMin(1);
3478           Short_t colmax       = fCalibraMode->GetColMax(1);
3479           Float_t vf           = fCurrentCoef[0];
3480           Float_t vs           = fCurrentCoef[1];
3481           Float_t vfE          = fCurrentCoefE;
3482           Float_t t0f          = fCurrentCoef2[0];
3483           Float_t t0s          = fCurrentCoef2[1];
3484           Float_t t0E          = fCurrentCoefE2;
3485
3486
3487
3488           (* fDebugStreamer) << "FillFillPH"<<
3489             "detector="<<detector<<
3490             "nentries="<<nentries<<
3491             "caligroup="<<caligroup<<
3492             "rowmin="<<rowmin<<
3493             "rowmax="<<rowmax<<
3494             "colmin="<<colmin<<
3495             "colmax="<<colmax<<
3496             "vf="<<vf<<
3497             "vs="<<vs<<
3498             "vfE="<<vfE<<
3499             "t0f="<<t0f<<
3500             "t0s="<<t0s<<
3501             "t0E="<<t0E<<
3502             "\n";
3503         }
3504         // Reset
3505         for (Int_t k = 0; k < 2304; k++) {
3506           fCurrentCoefDetector[k] = 0.0;
3507           fCurrentCoefDetector2[k] = 0.0;
3508         }
3509
3510       }// loop detector
3511       //printf("Check the count now: fCountDet %d\n",fCountDet);
3512     }
3513     else {
3514
3515       Int_t factor = 0;
3516       if(GetStack(fCountDet) == 2) factor = 12;
3517       else factor = 16;
3518
3519       for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) {
3520         for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) {
3521           fCurrentCoefDetector[(Int_t)(j*factor+k)]  = fCurrentCoef[0];
3522           fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0];
3523         }
3524       }
3525
3526       FillFillPH(idect,nentries);
3527     }
3528   }
3529   return kTRUE;
3530 }
3531 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3532 Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect)
3533 {
3534   //
3535   // Fill the coefficients found with the fits or other
3536   // methods from the Fit functions
3537   //
3538
3539   if (fDebugLevel != 1) {
3540     if (fNbDet > 0){
3541
3542       Int_t firstdetector = fCountDet;
3543       Int_t lastdetector  = fCountDet+fNbDet;
3544       AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
3545                    ,idect,firstdetector,lastdetector));
3546
3547       // loop over detectors
3548       for(Int_t det = firstdetector; det < lastdetector; det++){
3549
3550         //Set the calibration object again
3551         fCountDet = det;
3552         SetCalROC(2);
3553
3554         // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
3555         // Put them at 1
3556         fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2);
3557         fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet)
3558                                            ,(Int_t) GetStack(fCountDet)
3559                                            ,(Int_t) GetSector(fCountDet),2);
3560         // Reconstruct row min row max
3561         ReconstructFitRowMinRowMax(idect,2);
3562
3563         // Calcul the coef from the database choosen for the detector
3564         CalculPRFCoefMean();
3565
3566         //stack 2, not stack 2
3567         Int_t factor = 0;
3568         if(GetStack(fCountDet) == 2) factor = 12;
3569         else factor = 16;
3570
3571         // Fill the fCurrentCoefDetector with negative value to say: not fitted
3572         Double_t coeftoput = 1.0;
3573         if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]);
3574         else coeftoput = fCurrentCoef[0];
3575         for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3576           for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3577             fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput;
3578           }
3579         }
3580
3581         // Fill the stuff
3582         FillVectorFit();
3583         // Debug
3584         if(fDebugLevel > 1){
3585
3586           if ( !fDebugStreamer ) {
3587             //debug stream
3588             TDirectory *backup = gDirectory;
3589             fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3590             if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3591           }
3592
3593           Int_t   detector     = fCountDet;
3594           Int_t   layer        = GetLayer(fCountDet);
3595           Int_t   caligroup    = idect;
3596           Short_t rowmin       = fCalibraMode->GetRowMin(2);
3597           Short_t rowmax       = fCalibraMode->GetRowMax(2);
3598           Short_t colmin       = fCalibraMode->GetColMin(2);
3599           Short_t colmax       = fCalibraMode->GetColMax(2);
3600           Float_t widf         = fCurrentCoef[0];
3601           Float_t wids         = fCurrentCoef[1];
3602           Float_t widfE        = fCurrentCoefE;
3603
3604           (* fDebugStreamer) << "FillFillPRF"<<
3605             "detector="<<detector<<
3606             "layer="<<layer<<
3607             "caligroup="<<caligroup<<
3608             "rowmin="<<rowmin<<
3609             "rowmax="<<rowmax<<
3610             "colmin="<<colmin<<
3611             "colmax="<<colmax<<
3612             "widf="<<widf<<
3613             "wids="<<wids<<
3614             "widfE="<<widfE<<
3615             "\n";
3616         }
3617         // Reset
3618         for (Int_t k = 0; k < 2304; k++) {
3619           fCurrentCoefDetector[k] = 0.0;
3620         }
3621
3622       }// loop detector
3623       //printf("Check the count now: fCountDet %d\n",fCountDet);
3624     }
3625     else {
3626
3627       Int_t factor = 0;
3628       if(GetStack(fCountDet) == 2) factor = 12;
3629       else factor = 16;
3630
3631       // Pointer to the branch
3632       for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) {
3633         for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) {
3634           fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0];
3635         }
3636       }
3637       FillFillPRF(idect);
3638     }
3639   }
3640
3641   return kTRUE;
3642
3643 }
3644 //____________Functions for initialising the AliTRDCalibraFit in the code_________
3645 Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter()
3646 {
3647   //
3648   // Fill the coefficients found with the fits or other
3649   // methods from the Fit functions
3650   //
3651
3652   Int_t factor = 0;
3653   if(GetStack(fCountDet) == 2) factor = 1728;
3654   else factor = 2304;
3655
3656   // Pointer to the branch
3657   for (Int_t k = 0; k < factor; k++) {
3658     fCurrentCoefDetector[k]  = fCurrentCoef[0];
3659     fCurrentCoefDetector2[k] = fCurrentCoef2[0];
3660   }
3661
3662   FillFillLinearFitter();
3663
3664   return kTRUE;
3665
3666 }
3667 //________________________________________________________________________________
3668 void AliTRDCalibraFit::FillFillCH(Int_t idect)
3669 {
3670   //
3671   // DebugStream and fVectorFit
3672   //
3673
3674   // End of one detector
3675   if ((idect == (fCount-1))) {
3676     FillVectorFit();
3677     // Reset
3678     for (Int_t k = 0; k < 2304; k++) {
3679       fCurrentCoefDetector[k] = 0.0;
3680     }
3681   }
3682
3683   if(fDebugLevel > 1){
3684
3685     if ( !fDebugStreamer ) {
3686       //debug stream
3687       TDirectory *backup = gDirectory;
3688       fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
3689       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3690     }
3691
3692     Int_t   detector   = fCountDet;
3693     Int_t   caligroup  = idect;
3694     Short_t rowmin     = fCalibraMode->GetRowMin(0);
3695     Short_t rowmax     = fCalibraMode->GetRowMax(0);
3696     Short_t colmin     = fCalibraMode->GetColMin(0);
3697     Short_t colmax     = fCalibraMode->GetColMax(0);
3698     Float_t gf         = fCurrentCoef[0];
3699     Float_t gfs        = fCurrentCoef[1];
3700     Float_t gfE        = fCurrentCoefE;
3701
3702     (*fDebugStreamer) << "FillFillCH" <<
3703       "detector=" << detector <<
3704       "caligroup=" << caligroup <<
3705       "rowmin=" << rowmin <<
3706       "rowmax=" << rowmax <<
3707       "colmin=" << colmin <<
3708       "colmax=" << colmax <<
3709       "gf=" << gf <<
3710       "gfs=" << gfs <<
3711       "gfE=" << gfE <<
3712       "\n";
3713
3714   }
3715 }
3716 //________________________________________________________________________________
3717 void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries)
3718 {
3719   //
3720   // DebugStream and fVectorFit and fVectorFit2
3721   //
3722
3723   // End of one detector
3724     if ((idect == (fCount-1))) {
3725       FillVectorFit();
3726       FillVectorFit2();
3727       // Reset
3728       for (Int_t k = 0; k < 2304; k++) {
3729         fCurrentCoefDetector[k] = 0.0;
3730         fCurrentCoefDetector2[k] = 0.0;
3731       }
3732     }
3733
3734     if(fDebugLevel > 1){
3735
3736       if ( !fDebugStreamer ) {
3737         //debug stream
3738         TDirectory *backup = gDirectory;
3739         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root");
3740         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3741       }
3742
3743
3744       Int_t   detector     = fCountDet;
3745       Int_t   caligroup    = idect;
3746       Short_t rowmin       = fCalibraMode->GetRowMin(1);
3747       Short_t rowmax       = fCalibraMode->GetRowMax(1);
3748       Short_t colmin       = fCalibraMode->GetColMin(1);
3749       Short_t colmax       = fCalibraMode->GetColMax(1);
3750       Float_t vf           = fCurrentCoef[0];
3751       Float_t vs           = fCurrentCoef[1];
3752       Float_t vfE          = fCurrentCoefE;
3753       Float_t t0f          = fCurrentCoef2[0];
3754       Float_t t0s          = fCurrentCoef2[1];
3755       Float_t t0E          = fCurrentCoefE2;
3756
3757
3758
3759       (* fDebugStreamer) << "FillFillPH"<<
3760         "detector="<<detector<<
3761         "nentries="<<nentries<<
3762         "caligroup="<<caligroup<<
3763         "rowmin="<<rowmin<<
3764         "rowmax="<<rowmax<<
3765         "colmin="<<colmin<<
3766         "colmax="<<colmax<<
3767         "vf="<<vf<<
3768         "vs="<<vs<<
3769         "vfE="<<vfE<<
3770         "t0f="<<t0f<<
3771         "t0s="<<t0s<<
3772         "t0E="<<t0E<<
3773         "\n";
3774     }
3775
3776 }
3777 //________________________________________________________________________________
3778 void AliTRDCalibraFit::FillFillPRF(Int_t idect)
3779 {
3780   //
3781   // DebugStream and fVectorFit
3782   //
3783
3784     // End of one detector
3785     if ((idect == (fCount-1))) {
3786       FillVectorFit();
3787       // Reset
3788       for (Int_t k = 0; k < 2304; k++) {
3789         fCurrentCoefDetector[k] = 0.0;
3790       }
3791     }
3792
3793
3794     if(fDebugLevel > 1){
3795
3796       if ( !fDebugStreamer ) {
3797         //debug stream
3798         TDirectory *backup = gDirectory;
3799         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
3800         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3801       }
3802
3803       Int_t   detector     = fCountDet;
3804       Int_t   layer        = GetLayer(fCountDet);
3805       Int_t   caligroup    = idect;
3806       Short_t rowmin       = fCalibraMode->GetRowMin(2);
3807       Short_t rowmax       = fCalibraMode->GetRowMax(2);
3808       Short_t colmin       = fCalibraMode->GetColMin(2);
3809       Short_t colmax       = fCalibraMode->GetColMax(2);
3810       Float_t widf         = fCurrentCoef[0];
3811       Float_t wids         = fCurrentCoef[1];
3812       Float_t widfE        = fCurrentCoefE;
3813
3814       (* fDebugStreamer) << "FillFillPRF"<<
3815         "detector="<<detector<<
3816         "layer="<<layer<<
3817         "caligroup="<<caligroup<<
3818         "rowmin="<<rowmin<<
3819         "rowmax="<<rowmax<<
3820         "colmin="<<colmin<<
3821         "colmax="<<colmax<<
3822         "widf="<<widf<<
3823         "wids="<<wids<<
3824         "widfE="<<widfE<<
3825         "\n";
3826     }
3827
3828 }
3829 //________________________________________________________________________________
3830 void AliTRDCalibraFit::FillFillLinearFitter()
3831 {
3832   //
3833   // DebugStream and fVectorFit
3834   //
3835
3836   // End of one detector
3837   FillVectorFit();
3838   FillVectorFit2();
3839
3840
3841   // Reset
3842   for (Int_t k = 0; k < 2304; k++) {
3843   fCurrentCoefDetector[k]  = 0.0;
3844   fCurrentCoefDetector2[k] = 0.0;
3845   }
3846
3847
3848   if(fDebugLevel > 1){
3849
3850     if ( !fDebugStreamer ) {
3851       //debug stream
3852       TDirectory *backup = gDirectory;
3853       fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root");
3854       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
3855     }
3856
3857     //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
3858     AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
3859     Float_t rowmd            = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
3860     Float_t r                = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
3861     Float_t tiltangle        = padplane->GetTiltingAngle();
3862     Int_t   detector         = fCountDet;
3863     Int_t   stack            = GetStack(fCountDet);
3864     Int_t   layer            = GetLayer(fCountDet);
3865     Float_t vf               = fCurrentCoef[0];
3866     Float_t vs               = fCurrentCoef[1];
3867     Float_t vfE              = fCurrentCoefE;
3868     Float_t lorentzangler    = fCurrentCoef2[0];
3869     Float_t elorentzangler   = fCurrentCoefE2;
3870     Float_t lorentzangles    = fCurrentCoef2[1];
3871
3872     (* fDebugStreamer) << "FillFillLinearFitter"<<
3873       "detector="<<detector<<
3874       "stack="<<stack<<
3875       "layer="<<layer<<
3876       "rowmd="<<rowmd<<
3877       "r="<<r<<
3878       "tiltangle="<<tiltangle<<
3879       "vf="<<vf<<
3880       "vs="<<vs<<
3881       "vfE="<<vfE<<
3882       "lorentzangler="<<lorentzangler<<
3883       "Elorentzangler="<<elorentzangler<<
3884       "lorentzangles="<<lorentzangles<<
3885       "\n";
3886   }
3887
3888 }
3889 //
3890 //____________Calcul Coef Mean_________________________________________________
3891 //
3892 //_____________________________________________________________________________
3893 Bool_t AliTRDCalibraFit::CalculT0CoefMean()
3894 {
3895   //
3896   // For the detector Dect calcul the mean time 0
3897   // for the calibration group idect from the choosen database
3898   //
3899
3900   fCurrentCoef2[1] = 0.0;
3901   if(fDebugLevel != 1){
3902     if(((fCalibraMode->GetNz(1) > 0) ||
3903        (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1)    != 10) && (fCalibraMode->GetNz(1)    != 100))) {
3904
3905       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3906         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3907           fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3908         }
3909       }
3910
3911       fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
3912
3913     }
3914     else {
3915
3916       if(!fAccCDB){
3917         fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
3918       }
3919       else{
3920
3921         for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){
3922           for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){
3923             fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet));
3924           }
3925         }
3926         fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet))));
3927
3928       }
3929     }
3930   }
3931   return kTRUE;
3932 }
3933
3934 //_____________________________________________________________________________
3935 Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai)
3936 {
3937   //
3938   // For the detector Dect calcul the mean gain factor
3939   // for the calibration group idect from the choosen database
3940   //
3941
3942   fCurrentCoef[1] = 0.0;
3943   if(fDebugLevel != 1){
3944     if (((fCalibraMode->GetNz(0)    > 0) ||
3945         (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0)    != 10) && (fCalibraMode->GetNz(0)    != 100))) {
3946       for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) {
3947         for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) {
3948           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3949           if (vrai) fScaleFitFactor += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3950         }
3951       }
3952       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0)));
3953     }
3954     else {
3955       //Per detectors
3956       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
3957       if (vrai) fScaleFitFactor += ((Float_t) fCalDet->GetValue(fCountDet))*(fCalibraMode->GetColMax(0)-fCalibraMode->GetColMin(0))*(fCalibraMode->GetRowMax(0)-fCalibraMode->GetRowMin(0));
3958     }
3959   }
3960   return kTRUE;
3961 }
3962 //_____________________________________________________________________________
3963 Bool_t AliTRDCalibraFit::CalculPRFCoefMean()
3964 {
3965   //
3966   // For the detector Dect calcul the mean sigma of pad response
3967   // function for the calibration group idect from the choosen database
3968   //
3969
3970   fCurrentCoef[1] = 0.0;
3971   if(fDebugLevel != 1){
3972     for (Int_t row = fCalibraMode->GetRowMin(2); row < fCalibraMode->GetRowMax(2); row++) {
3973       for (Int_t col = fCalibraMode->GetColMin(2); col < fCalibraMode->GetColMax(2); col++) {
3974         fCurrentCoef[1] += (Float_t) fCalROC->GetValue(col,row);
3975       }
3976     }
3977     fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(2)-fCalibraMode->GetColMin(2))*(fCalibraMode->GetRowMax(2)-fCalibraMode->GetRowMin(2)));
3978   }
3979   return kTRUE;
3980 }
3981 //_____________________________________________________________________________
3982 Bool_t AliTRDCalibraFit::CalculVdriftCoefMean()
3983 {
3984   //
3985   // For the detector dect calcul the mean drift velocity for the
3986   // calibration group idect from the choosen database
3987   //
3988
3989   fCurrentCoef[1] = 0.0;
3990   if(fDebugLevel != 1){
3991     if (((fCalibraMode->GetNz(1)    > 0) ||
3992         (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1)    != 10) && (fCalibraMode->GetNz(1)    != 100))) {
3993
3994       for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) {
3995         for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) {
3996           fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet));
3997         }
3998       }
3999
4000       fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1)));
4001
4002     }
4003     else {
4004       //per detectors
4005       fCurrentCoef[1] = (Float_t) fCalDet->GetValue(fCountDet);
4006     }
4007   }
4008   return kTRUE;
4009 }
4010 //_____________________________________________________________________________
4011 Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef()
4012 {
4013   //
4014   // For the detector fCountDet, mean drift velocity and tan lorentzangle
4015   //
4016
4017   fCurrentCoef[1]  = fCalDet->GetValue(fCountDet);
4018   fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
4019
4020   return kTRUE;
4021 }
4022 //_____________________________________________________________________________
4023 Float_t AliTRDCalibraFit::GetPRFDefault(Int_t layer) const
4024 {
4025   //
4026   // Default width of the PRF if there is no database as reference
4027   //
4028   switch(layer)
4029     {
4030       // default database
4031       //case 0:  return 0.515;
4032       //case 1:  return 0.502;
4033       //case 2:  return 0.491;
4034       //case 3:  return 0.481;
4035       //case 4:  return 0.471;
4036       //case 5:  return 0.463;
4037       //default: return 0.0;
4038
4039       // fit database
4040     case 0:  return 0.538429;
4041     case 1:  return 0.524302;
4042     case 2:  return 0.511591;
4043     case 3:  return 0.500140;
4044     case 4:  return 0.489821;
4045     case 5:  return 0.480524;
4046     default: return 0.0;
4047   }
4048 }
4049 //________________________________________________________________________________
4050 void AliTRDCalibraFit::SetCalROC(Int_t i)
4051 {
4052   //
4053   // Set the calib object for fCountDet
4054   //
4055
4056   Float_t value = 0.0;
4057
4058   //Get the CalDet object
4059   if(fAccCDB){
4060     AliTRDcalibDB     *cal    = AliTRDcalibDB::Instance();
4061     if (!cal) {
4062       AliInfo("Could not get calibDB");
4063       return;
4064     }
4065     switch (i)
4066       {
4067       case 0:
4068         if( fCalROC ){
4069           fCalROC->~AliTRDCalROC();
4070           new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4071         }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet)));
4072         break;
4073       case 1:
4074         if( fCalROC ){
4075           fCalROC->~AliTRDCalROC();
4076           new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4077         }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet)));
4078         if( fCalROC2 ){
4079           fCalROC2->~AliTRDCalROC();
4080           new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4081         }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet)));
4082         break;
4083       case 2:
4084         if( fCalROC ){
4085           fCalROC->~AliTRDCalROC();
4086           new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4087         }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet)));
4088         break;
4089       default: return;
4090       }
4091   }
4092   else{
4093     switch (i)
4094       {
4095       case 0:
4096         if(fCalROC) delete fCalROC;
4097         fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4098         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4099           fCalROC->SetValue(k,1.0);
4100         }
4101         break;
4102       case 1:
4103         if(fCalROC)  delete fCalROC;
4104         if(fCalROC2) delete fCalROC2;
4105         fCalROC  = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4106         fCalROC2 = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4107         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4108           fCalROC->SetValue(k,1.0);
4109           fCalROC2->SetValue(k,0.0);
4110         }
4111         break;
4112       case 2:
4113         if(fCalROC) delete fCalROC;
4114         value = GetPRFDefault(GetLayer(fCountDet));
4115         fCalROC = new AliTRDCalROC(GetLayer(fCountDet),GetStack(fCountDet));
4116         for(Int_t k = 0; k < fCalROC->GetNchannels(); k++){
4117           fCalROC->SetValue(k,value);
4118         }
4119         break;
4120       default: return;
4121       }
4122   }
4123
4124 }
4125 //____________Fit Methods______________________________________________________
4126
4127 //_____________________________________________________________________________
4128 void AliTRDCalibraFit::FitPente(TH1* projPH)
4129 {
4130   //
4131   // Slope methode for the drift velocity
4132   //
4133
4134   // Constants
4135   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4136   Int_t binmax           = 0;
4137   Int_t binmin           = 0;
4138   fPhd[0]                = 0.0;
4139   fPhd[1]                = 0.0;
4140   fPhd[2]                = 0.0;
4141   Int_t ju               = 0;
4142   fCurrentCoefE          = 0.0;
4143   fCurrentCoefE2         = 0.0;
4144   fCurrentCoef[0]        = 0.0;
4145   fCurrentCoef2[0]       = 0.0;
4146   TLine *line            = new TLine();
4147
4148   // Some variables
4149   TAxis   *xpph    = projPH->GetXaxis();
4150   Int_t    nbins   = xpph->GetNbins();
4151   Double_t lowedge = xpph->GetBinLowEdge(1);
4152   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
4153   Double_t widbins = (upedge - lowedge) / nbins;
4154   Double_t limit   = upedge + 0.5 * widbins;
4155   Bool_t put = kTRUE;
4156
4157   // Beginning of the signal
4158   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4159   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
4160     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4161   }
4162   binmax = (Int_t) pentea->GetMaximumBin();
4163   if (binmax <= 1) {
4164     binmax = 2;
4165     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4166   }
4167   if (binmax >= nbins) {
4168     binmax = nbins-1;
4169     put = kFALSE;
4170     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4171   }
4172   pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
4173   Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
4174   Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
4175   Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
4176   Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
4177   if (TMath::Abs(l3P2am) > 0.00000001) {
4178     fPhd[0] = -(l3P1am / (2 * l3P2am));
4179   }
4180   if(!fTakeTheMaxPH){
4181     if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
4182       fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
4183     }
4184   }
4185   // Amplification region
4186   binmax = 0;
4187   ju     = 0;
4188   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4189     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4190       binmax = kbin;
4191       ju     = 1;
4192     }
4193   }
4194   if (binmax <= 1) {
4195     binmax = 2;
4196     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4197   }
4198   if (binmax >= nbins) {
4199     binmax = nbins-1;
4200     put = kFALSE;
4201     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4202   }
4203   projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
4204   Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
4205   Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
4206   Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
4207   Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
4208   if (TMath::Abs(l3P2amf) > 0.00000000001) {
4209     fPhd[1] = -(l3P1amf / (2 * l3P2amf));
4210   }
4211   if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
4212     fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
4213   }
4214   if(fTakeTheMaxPH){
4215     fCurrentCoefE2 = fCurrentCoefE;
4216   }
4217   // Drift region
4218   TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
4219   for (Int_t k = TMath::Min(binmax+4,projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
4220     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4221   }
4222   binmin = 0;
4223   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4224   if (binmin <= 1) {
4225     binmin = 2;
4226     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4227   }
4228   if (binmin >= nbins) {
4229     binmin = nbins-1;
4230     put = kFALSE;
4231     AliInfo("Put the binmax from nbins-1 to nbins-2 to enable the fit");
4232   }
4233   pente->Fit("pol2"
4234             ,"0MR"
4235             ,""
4236             ,TMath::Max(pente->GetBinCenter(binmin-1),             0.0)
4237             ,TMath::Min(pente->GetBinCenter(binmin+1),(Double_t) limit));
4238   Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
4239   Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
4240   Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
4241   Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
4242   if (TMath::Abs(l3P2dr) > 0.00000001) {
4243     fPhd[2] = -(l3P1dr / (2 * l3P2dr));
4244   }
4245   if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
4246     fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
4247   }
4248   Float_t fPhdt0  = 0.0;
4249   Float_t t0Shift = 0.0;
4250   if(fTakeTheMaxPH) {
4251     fPhdt0 = fPhd[1];
4252     t0Shift = fT0Shift1;
4253   }
4254   else {
4255     fPhdt0 = fPhd[0];
4256     t0Shift = fT0Shift0;
4257   }
4258
4259   if ((fPhd[2] > fPhd[0]) &&
4260       (fPhd[2] > fPhd[1]) &&
4261       (fPhd[1] > fPhd[0]) &&
4262       (put)) {
4263     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4264     fNumberFitSuccess++;
4265
4266     if (fPhdt0 >= 0.0) {
4267       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4268       if (fCurrentCoef2[0] < -1.0) {
4269         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4270       }
4271     }
4272     else {
4273       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4274     }
4275
4276   }
4277   else {
4278     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
4279     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4280   }
4281
4282   if (fDebugLevel == 1) {
4283     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4284     cpentei->cd();
4285     projPH->Draw();
4286     line->SetLineColor(2);
4287     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4288     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4289     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4290     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
4291     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
4292     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
4293     AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4294     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4295     cpentei2->cd();
4296     pentea->Draw();
4297     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4298     cpentei3->cd();
4299     pente->Draw();
4300   }
4301   else {
4302     delete pentea;
4303     delete pente;
4304   }
4305 }
4306 //_____________________________________________________________________________
4307 void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH)
4308 {
4309   //
4310   // Slope methode but with polynomes de Lagrange
4311   //
4312
4313   // Constants
4314   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4315   Int_t binmax      = 0;
4316   Int_t binmin      = 0;
4317   //Double_t    *x    = new Double_t[5];
4318   //Double_t    *y    = new Double_t[5];
4319   Double_t x[5];
4320   Double_t y[5];
4321   x[0]              = 0.0;
4322   x[1]              = 0.0;
4323   x[2]              = 0.0;
4324   x[3]              = 0.0;
4325   x[4]              = 0.0;
4326   y[0]              = 0.0;
4327   y[1]              = 0.0;
4328   y[2]              = 0.0;
4329   y[3]              = 0.0;
4330   y[4]              = 0.0;
4331   fPhd[0]           = 0.0;
4332   fPhd[1]           = 0.0;
4333   fPhd[2]           = 0.0;
4334   Int_t ju          = 0;
4335   fCurrentCoefE     = 0.0;
4336   fCurrentCoefE2    = 1.0;
4337   fCurrentCoef[0]   = 0.0;
4338   fCurrentCoef2[0]  = 0.0;
4339   TLine *line = new TLine();
4340   TF1 * polynome = 0x0;
4341   TF1 * polynomea = 0x0;
4342   TF1 * polynomeb = 0x0;
4343   Double_t c0 = 0.0;
4344   Double_t c1 = 0.0;
4345   Double_t c2 = 0.0;
4346   Double_t c3 = 0.0;
4347   Double_t c4 = 0.0;
4348
4349   // Some variables
4350   TAxis   *xpph    = projPH->GetXaxis();
4351   Int_t    nbins   = xpph->GetNbins();
4352   Double_t lowedge = xpph->GetBinLowEdge(1);
4353   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
4354   Double_t widbins = (upedge - lowedge) / nbins;
4355   Double_t limit   = upedge + 0.5 * widbins;
4356
4357
4358   Bool_t put = kTRUE;
4359
4360   // Beginning of the signal
4361   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4362   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
4363     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4364   }
4365
4366   binmax = (Int_t) pentea->GetMaximumBin();
4367
4368   Double_t minnn = 0.0;
4369   Double_t maxxx = 0.0;
4370
4371   Int_t kase = nbins-binmax;
4372
4373   switch(kase)
4374     {
4375     case 0:
4376       put = kFALSE;
4377       break;
4378     case 1:
4379       minnn = pentea->GetBinCenter(binmax-2);
4380       maxxx = pentea->GetBinCenter(binmax);
4381       x[0] = pentea->GetBinCenter(binmax-2);
4382       x[1] = pentea->GetBinCenter(binmax-1);
4383       x[2] = pentea->GetBinCenter(binmax);
4384       y[0] = pentea->GetBinContent(binmax-2);
4385       y[1] = pentea->GetBinContent(binmax-1);
4386       y[2] = pentea->GetBinContent(binmax);
4387       CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4388       AliInfo("At the limit for beginning!");
4389       break;
4390     case 2:
4391       minnn = pentea->GetBinCenter(binmax-2);
4392       maxxx = pentea->GetBinCenter(binmax+1);
4393       x[0] = pentea->GetBinCenter(binmax-2);
4394       x[1] = pentea->GetBinCenter(binmax-1);
4395       x[2] = pentea->GetBinCenter(binmax);
4396       x[3] = pentea->GetBinCenter(binmax+1);
4397       y[0] = pentea->GetBinContent(binmax-2);
4398       y[1] = pentea->GetBinContent(binmax-1);
4399       y[2] = pentea->GetBinContent(binmax);
4400       y[3] = pentea->GetBinContent(binmax+1);
4401       CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4402       break;
4403     default:
4404       switch(binmax){
4405       case 0:
4406         put = kFALSE;
4407         break;
4408       case 1:
4409         minnn = pentea->GetBinCenter(binmax);
4410         maxxx = pentea->GetBinCenter(binmax+2);
4411         x[0] = pentea->GetBinCenter(binmax);
4412         x[1] = pentea->GetBinCenter(binmax+1);
4413         x[2] = pentea->GetBinCenter(binmax+2);
4414         y[0] = pentea->GetBinContent(binmax);
4415         y[1] = pentea->GetBinContent(binmax+1);
4416         y[2] = pentea->GetBinContent(binmax+2);
4417         CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4418         break;
4419       case 2:
4420         minnn = pentea->GetBinCenter(binmax-1);
4421         maxxx = pentea->GetBinCenter(binmax+2);
4422         x[0] = pentea->GetBinCenter(binmax-1);
4423         x[1] = pentea->GetBinCenter(binmax);
4424         x[2] = pentea->GetBinCenter(binmax+1);
4425         x[3] = pentea->GetBinCenter(binmax+2);
4426         y[0] = pentea->GetBinContent(binmax-1);
4427         y[1] = pentea->GetBinContent(binmax);
4428         y[2] = pentea->GetBinContent(binmax+1);
4429         y[3] = pentea->GetBinContent(binmax+2);
4430         CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4431         break;
4432       default:
4433         minnn = pentea->GetBinCenter(binmax-2);
4434         maxxx = pentea->GetBinCenter(binmax+2);
4435         x[0] = pentea->GetBinCenter(binmax-2);
4436         x[1] = pentea->GetBinCenter(binmax-1);
4437         x[2] = pentea->GetBinCenter(binmax);
4438         x[3] = pentea->GetBinCenter(binmax+1);
4439         x[4] = pentea->GetBinCenter(binmax+2);
4440         y[0] = pentea->GetBinContent(binmax-2);
4441         y[1] = pentea->GetBinContent(binmax-1);
4442         y[2] = pentea->GetBinContent(binmax);
4443         y[3] = pentea->GetBinContent(binmax+1);
4444         y[4] = pentea->GetBinContent(binmax+2);
4445         CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4446         break;
4447       }
4448       break;
4449     }
4450
4451
4452   if(put) {
4453     polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4454     polynomeb->SetParameters(c0,c1,c2,c3,c4);
4455
4456     Double_t step = (maxxx-minnn)/10000;
4457     Double_t l = minnn;
4458     Double_t maxvalue = 0.0;
4459     Double_t placemaximum = minnn;
4460     for(Int_t o = 0; o < 10000; o++){
4461       if(o == 0) maxvalue = polynomeb->Eval(l);
4462       if(maxvalue < (polynomeb->Eval(l))){
4463         maxvalue = polynomeb->Eval(l);
4464         placemaximum = l;
4465       }
4466       l += step;
4467     }
4468     fPhd[0] = placemaximum;
4469   }
4470
4471   // Amplification region
4472   binmax = 0;
4473   ju     = 0;
4474   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4475     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4476       binmax = kbin;
4477       ju     = 1;
4478     }
4479   }
4480
4481   Double_t minn = 0.0;
4482   Double_t maxx = 0.0;
4483   x[0] = 0.0;
4484   x[1] = 0.0;
4485   x[2] = 0.0;
4486   x[3] = 0.0;
4487   x[4] = 0.0;
4488   y[0] = 0.0;
4489   y[1] = 0.0;
4490   y[2] = 0.0;
4491   y[3] = 0.0;
4492   y[4] = 0.0;
4493
4494   Int_t    kase1 = nbins - binmax;
4495
4496   //Determination of minn and maxx
4497   //case binmax = nbins
4498   //pol2
4499   switch(kase1)
4500     {
4501     case 0:
4502       minn = projPH->GetBinCenter(binmax-2);
4503       maxx = projPH->GetBinCenter(binmax);
4504       x[0] = projPH->GetBinCenter(binmax-2);
4505       x[1] = projPH->GetBinCenter(binmax-1);
4506       x[2] = projPH->GetBinCenter(binmax);
4507       y[0] = projPH->GetBinContent(binmax-2);
4508       y[1] = projPH->GetBinContent(binmax-1);
4509       y[2] = projPH->GetBinContent(binmax);
4510       CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4511       //AliInfo("At the limit for the drift!");
4512       break;
4513     case 1:
4514       minn = projPH->GetBinCenter(binmax-2);
4515       maxx = projPH->GetBinCenter(binmax+1);
4516       x[0] = projPH->GetBinCenter(binmax-2);
4517       x[1] = projPH->GetBinCenter(binmax-1);
4518       x[2] = projPH->GetBinCenter(binmax);
4519       x[3] = projPH->GetBinCenter(binmax+1);
4520       y[0] = projPH->GetBinContent(binmax-2);
4521       y[1] = projPH->GetBinContent(binmax-1);
4522       y[2] = projPH->GetBinContent(binmax);
4523       y[3] = projPH->GetBinContent(binmax+1);
4524       CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4525       break;
4526     default:
4527       switch(binmax)
4528         {
4529         case 0:
4530           put = kFALSE;
4531           break;
4532         case 1:
4533           minn = projPH->GetBinCenter(binmax);
4534           maxx = projPH->GetBinCenter(binmax+2);
4535           x[0] = projPH->GetBinCenter(binmax);
4536           x[1] = projPH->GetBinCenter(binmax+1);
4537           x[2] = projPH->GetBinCenter(binmax+2);
4538           y[0] = projPH->GetBinContent(binmax);
4539           y[1] = projPH->GetBinContent(binmax+1);
4540           y[2] = projPH->GetBinContent(binmax+2);
4541           CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4542           break;
4543         case 2:
4544           minn = projPH->GetBinCenter(binmax-1);
4545           maxx = projPH->GetBinCenter(binmax+2);
4546           x[0] = projPH->GetBinCenter(binmax-1);
4547           x[1] = projPH->GetBinCenter(binmax);
4548           x[2] = projPH->GetBinCenter(binmax+1);
4549           x[3] = projPH->GetBinCenter(binmax+2);
4550           y[0] = projPH->GetBinContent(binmax-1);
4551           y[1] = projPH->GetBinContent(binmax);
4552           y[2] = projPH->GetBinContent(binmax+1);
4553           y[3] = projPH->GetBinContent(binmax+2);
4554           CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4555           break;
4556         default:
4557           minn = projPH->GetBinCenter(binmax-2);
4558           maxx = projPH->GetBinCenter(binmax+2);
4559           x[0] = projPH->GetBinCenter(binmax-2);
4560           x[1] = projPH->GetBinCenter(binmax-1);
4561           x[2] = projPH->GetBinCenter(binmax);
4562           x[3] = projPH->GetBinCenter(binmax+1);
4563           x[4] = projPH->GetBinCenter(binmax+2);
4564           y[0] = projPH->GetBinContent(binmax-2);
4565           y[1] = projPH->GetBinContent(binmax-1);
4566           y[2] = projPH->GetBinContent(binmax);
4567           y[3] = projPH->GetBinContent(binmax+1);
4568           y[4] = projPH->GetBinContent(binmax+2);
4569           CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4570           break;
4571         }
4572       break;
4573     }
4574
4575   if(put) {
4576     polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4577     polynomea->SetParameters(c0,c1,c2,c3,c4);
4578
4579     Double_t step = (maxx-minn)/1000;
4580     Double_t l = minn;
4581     Double_t maxvalue = 0.0;
4582     Double_t placemaximum = minn;
4583     for(Int_t o = 0; o < 1000; o++){
4584       if(o == 0) maxvalue = polynomea->Eval(l);
4585       if(maxvalue < (polynomea->Eval(l))){
4586         maxvalue = polynomea->Eval(l);
4587         placemaximum = l;
4588       }
4589       l += step;
4590     }
4591     fPhd[1] = placemaximum;
4592   }
4593
4594   // Drift region
4595   TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4596   for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
4597     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4598   }
4599   binmin = 0;
4600   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4601
4602   //should not happen
4603   if (binmin <= 1) {
4604     binmin = 2;
4605     put = 1;
4606     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4607   }
4608
4609   //check
4610   if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4611     AliInfo("Too many fluctuations at the end!");
4612     put = kFALSE;
4613   }
4614   if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4615     AliInfo("Too many fluctuations at the end!");
4616     put = kFALSE;
4617   }
4618   if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4619     AliInfo("No entries for the next bin!");
4620     pente->SetBinContent(binmin,0);
4621     if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4622   }
4623
4624
4625   x[0] = 0.0;
4626   x[1] = 0.0;
4627   x[2] = 0.0;
4628   x[3] = 0.0;
4629   x[4] = 0.0;
4630   y[0] = 0.0;
4631   y[1] = 0.0;
4632   y[2] = 0.0;
4633   y[3] = 0.0;
4634   y[4] = 0.0;
4635   Double_t min = 0.0;
4636   Double_t max = 0.0;
4637   Bool_t case1 = kFALSE;
4638   Bool_t case2 = kFALSE;
4639   Bool_t case4 = kFALSE;
4640
4641   //Determination of min and max
4642   //case binmin <= nbins-3
4643   //pol4 case 3
4644   if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4645     min = pente->GetBinCenter(binmin-2);
4646     max = pente->GetBinCenter(binmin+2);
4647     x[0] = pente->GetBinCenter(binmin-2);
4648     x[1] = pente->GetBinCenter(binmin-1);
4649     x[2] = pente->GetBinCenter(binmin);
4650     x[3] = pente->GetBinCenter(binmin+1);
4651     x[4] = pente->GetBinCenter(binmin+2);
4652     y[0] = pente->GetBinContent(binmin-2);
4653     y[1] = pente->GetBinContent(binmin-1);
4654     y[2] = pente->GetBinContent(binmin);
4655     y[3] = pente->GetBinContent(binmin+1);
4656     y[4] = pente->GetBinContent(binmin+2);
4657     //Calcul the polynome de Lagrange
4658     CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
4659     //richtung +/-
4660     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4661        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4662       //AliInfo("polynome 4 false 1");
4663       put = kFALSE;
4664     }
4665     if(((binmin+3) <= (nbins-1)) &&
4666        (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4667        ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4668        (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4669       AliInfo("polynome 4 false 2");
4670       put = kFALSE;
4671     }
4672     // poly 3
4673     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4674        (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4675       //AliInfo("polynome 4 case 1");
4676       case1 = kTRUE;
4677     }
4678     if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4679        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4680       //AliInfo("polynome 4 case 4");
4681       case4 = kTRUE;
4682     }
4683
4684   }
4685   //case binmin = nbins-2
4686   //pol3 case 1
4687   if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4688      (case1)){
4689     min = pente->GetBinCenter(binmin-2);
4690     max = pente->GetBinCenter(binmin+1);
4691     x[0] = pente->GetBinCenter(binmin-2);
4692     x[1] = pente->GetBinCenter(binmin-1);
4693     x[2] = pente->GetBinCenter(binmin);
4694     x[3] = pente->GetBinCenter(binmin+1);
4695     y[0] = pente->GetBinContent(binmin-2);
4696     y[1] = pente->GetBinContent(binmin-1);
4697     y[2] = pente->GetBinContent(binmin);
4698     y[3] = pente->GetBinContent(binmin+1);
4699     //Calcul the polynome de Lagrange
4700     CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4701     //richtung +: nothing
4702     //richtung -
4703     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4704       //AliInfo("polynome 3- case 2");
4705       case2 = kTRUE;
4706     }
4707   }
4708   //pol3 case 4
4709   if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4710      (case4)){
4711     min = pente->GetBinCenter(binmin-1);
4712     max = pente->GetBinCenter(binmin+2);
4713     x[0] = pente->GetBinCenter(binmin-1);
4714     x[1] = pente->GetBinCenter(binmin);
4715     x[2] = pente->GetBinCenter(binmin+1);
4716     x[3] = pente->GetBinCenter(binmin+2);
4717     y[0] = pente->GetBinContent(binmin-1);
4718     y[1] = pente->GetBinContent(binmin);
4719     y[2] = pente->GetBinContent(binmin+1);
4720     y[3] = pente->GetBinContent(binmin+2);
4721     //Calcul the polynome de Lagrange
4722     CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
4723     //richtung +
4724     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4725       //AliInfo("polynome 3+ case 2");
4726       case2 = kTRUE;
4727     }
4728   }
4729   //pol2 case 5
4730   if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4731     min = pente->GetBinCenter(binmin);
4732     max = pente->GetBinCenter(binmin+2);
4733     x[0] = pente->GetBinCenter(binmin);
4734     x[1] = pente->GetBinCenter(binmin+1);
4735     x[2] = pente->GetBinCenter(binmin+2);
4736     y[0] = pente->GetBinContent(binmin);
4737     y[1] = pente->GetBinContent(binmin+1);
4738     y[2] = pente->GetBinContent(binmin+2);
4739     //Calcul the polynome de Lagrange
4740     CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4741     //richtung +
4742     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4743       //AliInfo("polynome 2+ false");
4744       put = kFALSE;
4745     }
4746   }
4747   //pol2 case 2
4748   if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4749      (case2)){
4750     min = pente->GetBinCenter(binmin-1);
4751     max = pente->GetBinCenter(binmin+1);
4752     x[0] = pente->GetBinCenter(binmin-1);
4753     x[1] = pente->GetBinCenter(binmin);
4754     x[2] = pente->GetBinCenter(binmin+1);
4755     y[0] = pente->GetBinContent(binmin-1);
4756     y[1] = pente->GetBinContent(binmin);
4757     y[2] = pente->GetBinContent(binmin+1);
4758     //Calcul the polynome de Lagrange
4759     CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4760     //richtung +: nothing
4761     //richtung -: nothing
4762   }
4763   //case binmin = nbins-1
4764   //pol2 case 0
4765   if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4766     min = pente->GetBinCenter(binmin-2);
4767     max = pente->GetBinCenter(binmin);
4768     x[0] = pente->GetBinCenter(binmin-2);
4769     x[1] = pente->GetBinCenter(binmin-1);
4770     x[2] = pente->GetBinCenter(binmin);
4771     y[0] = pente->GetBinContent(binmin-2);
4772     y[1] = pente->GetBinContent(binmin-1);
4773     y[2] = pente->GetBinContent(binmin);
4774     //Calcul the polynome de Lagrange
4775     CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
4776     //AliInfo("At the limit for the drift!");
4777     //fluctuation too big!
4778     //richtung +: nothing
4779     //richtung -
4780     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4781       //AliInfo("polynome 2- false ");
4782       put = kFALSE;
4783     }
4784   }
4785   if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4786     put = kFALSE;
4787     AliInfo("At the limit for the drift and not usable!");
4788   }
4789
4790   //pass
4791   if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4792     put = kFALSE;
4793     AliInfo("For the drift...problem!");
4794   }
4795   //pass but should not happen
4796   if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4797     put = kFALSE;
4798     AliInfo("For the drift...problem!");
4799   }
4800
4801   if(put) {
4802     polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4803     polynome->SetParameters(c0,c1,c2,c3,c4);
4804     //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4805     Double_t step = (max-min)/1000;
4806     Double_t l = min;
4807     Double_t minvalue = 0.0;
4808     Double_t placeminimum = min;
4809     for(Int_t o = 0; o < 1000; o++){
4810       if(o == 0) minvalue = polynome->Eval(l);
4811       if(minvalue > (polynome->Eval(l))){
4812         minvalue = polynome->Eval(l);
4813         placeminimum = l;
4814       }
4815       l += step;
4816     }
4817     fPhd[2] = placeminimum;
4818   }
4819   //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4820   if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4821   if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4822
4823   Float_t fPhdt0  = 0.0;
4824   Float_t t0Shift = 0.0;
4825   if(fTakeTheMaxPH) {
4826     fPhdt0 = fPhd[1];
4827     t0Shift = fT0Shift1;
4828   }
4829   else {
4830     fPhdt0 = fPhd[0];
4831     t0Shift = fT0Shift0;
4832   }
4833
4834   if ((fPhd[2] > fPhd[0]) &&
4835       (fPhd[2] > fPhd[1]) &&
4836       (fPhd[1] > fPhd[0]) &&
4837       (put)) {
4838     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4839     if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] =  -TMath::Abs(fCurrentCoef[1]);
4840     else fNumberFitSuccess++;
4841     if (fPhdt0 >= 0.0) {
4842       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4843       if (fCurrentCoef2[0] < -1.0) {
4844         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4845       }
4846     }
4847     else {
4848       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4849     }
4850   }
4851   else {
4852     //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4853     fCurrentCoef[0]      = -TMath::Abs(fCurrentCoef[1]);
4854
4855     if((fPhd[1] > fPhd[0]) &&
4856        (put)) {
4857       if (fPhdt0 >= 0.0) {
4858         fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4859         if (fCurrentCoef2[0] < -1.0) {
4860           fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4861         }
4862       }
4863       else {
4864         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4865       }
4866     }
4867     else{
4868       fCurrentCoef2[0]     = fCurrentCoef2[1] + 100.0;
4869       //printf("Fit failed!\n");
4870     }
4871   }
4872
4873   if (fDebugLevel == 1) {
4874     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4875     cpentei->cd();
4876     projPH->Draw();
4877     line->SetLineColor(2);
4878     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4879     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4880     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4881     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
4882     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
4883     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
4884     AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4885     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4886     cpentei2->cd();
4887     pentea->Draw();
4888     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4889     cpentei3->cd();
4890     pente->Draw();
4891   }
4892   else {
4893     delete pentea;
4894     delete pente;
4895     if(polynome) delete polynome;
4896     if(polynomea) delete polynomea;
4897     if(polynomeb) delete polynomeb;
4898     //if(x) delete [] x;
4899     //if(y) delete [] y;
4900     if(line) delete line;
4901
4902   }
4903
4904   //Provisoire
4905   //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4906   //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4907
4908   projPH->SetDirectory(0);
4909
4910 }
4911
4912 //_____________________________________________________________________________
4913 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4914 {
4915   //
4916   // Fit methode for the drift velocity
4917   //
4918
4919   // Constants
4920   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4921
4922   // Some variables
4923   TAxis   *xpph   = projPH->GetXaxis();
4924   Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4925
4926   TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4927   fPH->SetParameter(0,0.469);     // Scaling
4928   fPH->SetParameter(1,0.18);      // Start
4929   fPH->SetParameter(2,0.0857325); // AR
4930   fPH->SetParameter(3,1.89);      // DR
4931   fPH->SetParameter(4,0.08);      // QA/QD
4932   fPH->SetParameter(5,0.0);       // Baseline
4933
4934   TLine *line = new TLine();
4935
4936   fCurrentCoef[0]     = 0.0;
4937   fCurrentCoef2[0]    = 0.0;
4938   fCurrentCoefE       = 0.0;
4939   fCurrentCoefE2      = 0.0;
4940
4941   if (idect%fFitPHPeriode == 0) {
4942
4943     AliInfo(Form("The detector %d will be fitted",idect));
4944     fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4945     fPH->SetParameter(1,fPhd[0] - 0.1);                                                                 // Start
4946     fPH->SetParameter(2,fPhd[1] - fPhd[0]);                                                             // AR
4947     fPH->SetParameter(3,fPhd[2] - fPhd[1]);                                                             // DR
4948     fPH->SetParameter(4,0.225);                                                                         // QA/QD
4949     fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4950
4951     if (fDebugLevel != 1) {
4952       projPH->Fit(fPH,"0M","",0.0,upedge);
4953     }
4954     else {
4955       TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4956       cpente->cd();
4957       projPH->Fit(fPH,"M+","",0.0,upedge);
4958       projPH->Draw("E0");
4959       line->SetLineColor(4);
4960       line->DrawLine(fPH->GetParameter(1)
4961                     ,0
4962                     ,fPH->GetParameter(1)
4963                     ,projPH->GetMaximum());
4964       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4965                     ,0
4966                     ,fPH->GetParameter(1)+fPH->GetParameter(2)
4967                     ,projPH->GetMaximum());
4968       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4969                     ,0
4970                     ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4971                     ,projPH->GetMaximum());
4972     }
4973
4974     if (fPH->GetParameter(3) != 0) {
4975       fNumberFitSuccess++;
4976       fCurrentCoef[0]    = kDrWidth / (fPH->GetParameter(3));
4977       fCurrentCoefE      = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4978       fCurrentCoef2[0]   = fPH->GetParameter(1);
4979       fCurrentCoefE2     = fPH->GetParError(1);
4980     }
4981     else {
4982       fCurrentCoef[0]     = -TMath::Abs(fCurrentCoef[1]);
4983       fCurrentCoef2[0]    = fCurrentCoef2[1] + 100.0;
4984     }
4985
4986   }
4987   else {
4988
4989     // Put the default value
4990     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
4991     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4992   }
4993
4994   if (fDebugLevel != 1) {
4995     delete fPH;
4996   }
4997
4998 }
4999 //_____________________________________________________________________________
5000 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
5001 {
5002   //
5003   // Fit methode for the sigma of the pad response function
5004   //
5005
5006   TVectorD param(3);
5007
5008   fCurrentCoef[0]  = 0.0;
5009   fCurrentCoefE = 0.0;
5010
5011   Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
5012
5013   if(TMath::Abs(ret+4) <= 0.000000001){
5014     fCurrentCoef[0] = -fCurrentCoef[1];
5015     return kFALSE;
5016   }
5017   else {
5018     fNumberFitSuccess++;
5019     fCurrentCoef[0] = param[2];
5020     fCurrentCoefE   = ret;
5021     return kTRUE;
5022   }
5023 }
5024 //_____________________________________________________________________________
5025 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)
5026 {
5027   //
5028   // Fit methode for the sigma of the pad response function
5029   //
5030
5031   //We should have at least 3 points
5032   if(nBins <=3) return -4.0;
5033
5034   TLinearFitter fitter(3,"pol2");
5035   fitter.StoreData(kFALSE);
5036   fitter.ClearPoints();
5037   TVectorD  par(3);
5038   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5039   Float_t entries = 0;
5040   Int_t   nbbinwithentries = 0;
5041   for (Int_t i=0; i<nBins; i++){
5042     entries+=arraye[i];
5043     if(arraye[i] > 15) nbbinwithentries++;
5044     //printf("entries for i %d: %f\n",i,arraye[i]);
5045   }
5046   if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5047   //printf("entries %f\n",entries);
5048   //printf("nbbinwithentries %d\n",nbbinwithentries);
5049
5050   Int_t npoints=0;
5051   Float_t errorm = 0.0;
5052   Float_t errorn = 0.0;
5053   Float_t error  = 0.0;
5054
5055   //
5056   for (Int_t ibin=0;ibin<nBins; ibin++){
5057       Float_t entriesI = arraye[ibin];
5058       Float_t valueI   = arraym[ibin];
5059       Double_t xcenter = 0.0;
5060       Float_t  val     = 0.0;
5061       if ((entriesI>15) && (valueI>0.0)){
5062         xcenter = xMin+(ibin+0.5)*binWidth;
5063         errorm   = 0.0;
5064         errorn   = 0.0;
5065         error    = 0.0;
5066         if(!bError){
5067           if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5068           if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5069           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5070         }
5071         else{
5072           if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5073           if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5074           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5075         }
5076         if(TMath::Abs(error) < 0.000000001) continue;
5077         val      = TMath::Log(Float_t(valueI));
5078         fitter.AddPoint(&xcenter,val,error);
5079         npoints++;
5080       }
5081
5082       if(fDebugLevel > 1){
5083
5084       if ( !fDebugStreamer ) {
5085         //debug stream
5086         TDirectory *backup = gDirectory;
5087         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5088         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5089       }
5090
5091       Int_t    detector     = fCountDet;
5092       Int_t    layer        = GetLayer(fCountDet);
5093       Int_t    group        = ibin;
5094
5095       (* fDebugStreamer) << "FitGausMIFill"<<
5096         "detector="<<detector<<
5097         "layer="<<layer<<
5098         "nbins="<<nBins<<
5099         "group="<<group<<
5100         "entriesI="<<entriesI<<
5101         "valueI="<<valueI<<
5102         "val="<<val<<
5103         "xcenter="<<xcenter<<
5104         "errorm="<<errorm<<
5105         "errorn="<<errorn<<
5106         "error="<<error<<
5107         "bError="<<bError<<
5108         "\n";
5109     }
5110
5111   }
5112
5113   if(npoints <=3) return -4.0;
5114
5115   Double_t chi2 = 0;
5116   if (npoints>3){
5117     fitter.Eval();
5118     fitter.GetParameters(par);
5119     chi2 = fitter.GetChisquare()/Float_t(npoints);
5120
5121
5122     if (!param)  param  = new TVectorD(3);
5123     if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5124     Double_t  x      = TMath::Sqrt(TMath::Abs(-2*par[2]));
5125     Double_t deltax = (fitter.GetParError(2))/x;
5126     Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5127     chi2 = errorparam2;
5128
5129     (*param)[1] = par[1]/(-2.*par[2]);
5130     (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5131     Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
5132     if ( lnparam0>307 ) return -4;
5133     (*param)[0] = TMath::Exp(lnparam0);
5134
5135     if(fDebugLevel > 1){
5136
5137       if ( !fDebugStreamer ) {
5138         //debug stream
5139         TDirectory *backup = gDirectory;
5140         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5141         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5142       }
5143
5144       Int_t    detector     = fCountDet;
5145       Int_t    layer        = GetLayer(fCountDet);
5146
5147
5148       (* fDebugStreamer) << "FitGausMIFit"<<
5149         "detector="<<detector<<
5150         "layer="<<layer<<
5151         "nbins="<<nBins<<
5152         "errorsigma="<<chi2<<
5153         "mean="<<(*param)[1]<<
5154         "sigma="<<(*param)[2]<<
5155         "constant="<<(*param)[0]<<
5156         "\n";
5157     }
5158   }
5159
5160   if((chi2/(*param)[2]) > 0.1){
5161     if(bError){
5162       chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5163     }
5164     else return -4.0;
5165   }
5166
5167   if(fDebugLevel == 1){
5168     TString name("PRF");
5169     name += (Int_t)xMin;
5170     name += (Int_t)xMax;
5171     TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5172     c1->cd();
5173     name += "histo";
5174     TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5175     for(Int_t k = 0; k < nBins; k++){
5176       histo->SetBinContent(k+1,arraym[k]);
5177       histo->SetBinError(k+1,arrayme[k]);
5178     }
5179     histo->Draw();
5180     name += "functionf";
5181     TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5182     f1->SetParameter(0, (*param)[0]);
5183     f1->SetParameter(1, (*param)[1]);
5184     f1->SetParameter(2, (*param)[2]);
5185     f1->Draw("same");
5186   }
5187
5188
5189   return chi2;
5190
5191 }
5192 //_____________________________________________________________________________
5193 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5194 {
5195   //
5196   // Fit methode for the sigma of the pad response function
5197   //
5198
5199   fCurrentCoef[0]  = 0.0;
5200   fCurrentCoefE = 0.0;
5201
5202   if (fDebugLevel != 1) {
5203     projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5204   }
5205   else {
5206     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5207     cfit->cd();
5208     projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5209     projPRF->Draw();
5210   }
5211   fCurrentCoef[0]  = projPRF->GetFunction("gaus")->GetParameter(2);
5212   fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5213   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5214   else {
5215     fNumberFitSuccess++;
5216   }
5217 }
5218 //_____________________________________________________________________________
5219 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5220 {
5221   //
5222   // Fit methode for the sigma of the pad response function
5223   //
5224   fCurrentCoef[0]   = 0.0;
5225   fCurrentCoefE  = 0.0;
5226   if (fDebugLevel == 1) {
5227     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5228     cfit->cd();
5229     projPRF->Draw();
5230   }
5231   fCurrentCoef[0] = projPRF->GetRMS();
5232   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5233   else {
5234     fNumberFitSuccess++;
5235   }
5236 }
5237 //_____________________________________________________________________________
5238 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5239 {
5240   //
5241   // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5242   //
5243
5244   TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5245
5246
5247   Int_t   nbins    = (Int_t)(nybins/(2*nbg));
5248   Float_t lowedge  = -3.0*nbg;
5249   Float_t upedge   = lowedge + 3.0;
5250   Int_t   offset   = 0;
5251   Int_t   npoints  = 0;
5252   Double_t xvalues = -0.2*nbg+0.1;
5253   Double_t y       = 0.0;
5254   Int_t   total    = 2*nbg;
5255
5256
5257   for(Int_t k = 0; k < total; k++){
5258     if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5259       npoints++;
5260       y = fCurrentCoef[0]*fCurrentCoef[0];
5261       linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5262     }
5263
5264     if(fDebugLevel > 1){
5265
5266       if ( !fDebugStreamer ) {
5267         //debug stream
5268         TDirectory *backup = gDirectory;
5269         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5270         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5271       }
5272
5273       Int_t    detector     = fCountDet;
5274       Int_t    layer        = GetLayer(fCountDet);
5275       Int_t    nbtotal      = total;
5276       Int_t    group        = k;
5277       Float_t  low          = lowedge;
5278       Float_t  up           = upedge;
5279       Float_t  tnp          = xvalues;
5280       Float_t  wid          = fCurrentCoef[0];
5281       Float_t  widfE        = fCurrentCoefE;
5282
5283       (* fDebugStreamer) << "FitTnpRange0"<<
5284         "detector="<<detector<<
5285         "layer="<<layer<<
5286         "nbtotal="<<nbtotal<<
5287         "group="<<group<<
5288         "low="<<low<<
5289         "up="<<up<<
5290         "offset="<<offset<<
5291         "tnp="<<tnp<<
5292         "wid="<<wid<<
5293         "widfE="<<widfE<<
5294         "\n";
5295     }
5296
5297     offset  += nbins;
5298     lowedge += 3.0;
5299     upedge  += 3.0;
5300     xvalues += 0.2;
5301
5302   }
5303
5304   fCurrentCoefE = 0.0;
5305   fCurrentCoef[0] = 0.0;
5306
5307   //printf("npoints\n",npoints);
5308
5309   if(npoints < 3){
5310     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5311   }
5312   else{
5313
5314     TVectorD pars0;
5315     linearfitter.Eval();
5316     linearfitter.GetParameters(pars0);
5317     Double_t pointError0  =  TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5318     Double_t errorsx0     =  linearfitter.GetParError(2)*pointError0;
5319     Double_t min0         = 0.0;
5320     Double_t ermin0       = 0.0;
5321     //Double_t prfe0      = 0.0;
5322     Double_t prf0         = 0.0;
5323     if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5324       min0 = -pars0[1]/(2*pars0[2]);
5325       ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5326       prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5327       if(prf0 > 0.0) {
5328         /*
5329           prfe0 = linearfitter->GetParError(0)*pointError0
5330           +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5331           +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5332           prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5333           fCurrentCoefE   = (Float_t) prfe0;
5334         */
5335         fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5336       }
5337       else{
5338         fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5339       }
5340     }
5341     else {
5342       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5343     }
5344
5345     if(fDebugLevel > 1){
5346
5347       if ( !fDebugStreamer ) {
5348         //debug stream
5349         TDirectory *backup = gDirectory;
5350         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5351         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5352       }
5353
5354       Int_t    detector     = fCountDet;
5355       Int_t    layer        = GetLayer(fCountDet);
5356       Int_t    nbtotal      = total;
5357       Double_t colsize[6]   = {0.635,0.665,0.695,0.725,0.755,0.785};
5358       Double_t sigmax       = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5359
5360       (* fDebugStreamer) << "FitTnpRange1"<<
5361         "detector="<<detector<<
5362         "layer="<<layer<<
5363         "nbtotal="<<nbtotal<<
5364         "par0="<<pars0[0]<<
5365         "par1="<<pars0[1]<<
5366         "par2="<<pars0[2]<<
5367         "npoints="<<npoints<<
5368         "sigmax="<<sigmax<<
5369         "tan="<<min0<<
5370         "sigmaprf="<<fCurrentCoef[0]<<
5371         "sigprf="<<fCurrentCoef[1]<<
5372         "\n";
5373     }
5374
5375   }
5376
5377 }
5378 //_____________________________________________________________________________
5379 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5380 {
5381   //
5382   // Only mean methode for the gain factor
5383   //
5384
5385   fCurrentCoef[0] = mean;
5386   fCurrentCoefE   = 0.0;
5387   if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5388   if (fDebugLevel == 1) {
5389     TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5390     cpmean->cd();
5391     projch->Draw();
5392   }
5393   CalculChargeCoefMean(kTRUE);
5394   fNumberFitSuccess++;
5395 }
5396 //_____________________________________________________________________________
5397 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5398 {
5399   //
5400   // mean w methode for the gain factor
5401   //
5402
5403   //Number of bins
5404   Int_t nybins = projch->GetNbinsX();
5405
5406   //The weight function
5407   Double_t a = 0.00228515;
5408   Double_t b = -0.00231487;
5409   Double_t c = 0.00044298;
5410   Double_t d = -0.00379239;
5411   Double_t e = 0.00338349;
5412
5413 //   0 |0.00228515
5414 //    1 |-0.00231487
5415 //    2 |0.00044298
5416 //    3 |-0.00379239
5417 //    4 |0.00338349
5418
5419
5420
5421   //A arbitrary error for the moment
5422   fCurrentCoefE = 0.0;
5423   fCurrentCoef[0] = 0.0;
5424
5425   //Calcul
5426   Double_t sumw = 0.0;
5427   Double_t sum = 0.0;
5428   Float_t sumAll   = (Float_t) nentries;
5429   Int_t sumCurrent = 0;
5430   for(Int_t k = 0; k <nybins; k++){
5431     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5432     if (fraction>0.95) break;
5433     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5434       e*fraction*fraction*fraction*fraction;
5435     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5436     sum  += weight*projch->GetBinContent(k+1);
5437     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5438     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5439   }
5440   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5441
5442   if (fDebugLevel == 1) {
5443     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5444     cpmeanw->cd();
5445     projch->Draw();
5446   }
5447   fNumberFitSuccess++;
5448   CalculChargeCoefMean(kTRUE);
5449 }
5450 //_____________________________________________________________________________
5451 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5452 {
5453   //
5454   // mean w methode for the gain factor
5455   //
5456
5457   //Number of bins
5458   Int_t nybins = projch->GetNbinsX();
5459
5460   //The weight function
5461   Double_t a = 0.00228515;
5462   Double_t b = -0.00231487;
5463   Double_t c = 0.00044298;
5464   Double_t d = -0.00379239;
5465   Double_t e = 0.00338349;
5466
5467 //   0 |0.00228515
5468 //    1 |-0.00231487
5469 //    2 |0.00044298
5470 //    3 |-0.00379239
5471 //    4 |0.00338349
5472
5473
5474
5475   //A arbitrary error for the moment
5476   fCurrentCoefE = 0.0;
5477   fCurrentCoef[0] = 0.0;
5478
5479   //Calcul
5480   Double_t sumw = 0.0;
5481   Double_t sum = 0.0;
5482   Int_t sumCurrent = 0;
5483   for(Int_t k = 0; k <nybins; k++){
5484     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5485     if (fraction>0.95) break;
5486     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5487       e*fraction*fraction*fraction*fraction;
5488     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5489     sum  += weight*projch->GetBinContent(k+1);
5490     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5491     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5492   }
5493   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5494
5495   if (fDebugLevel == 1) {
5496     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5497     cpmeanw->cd();
5498     projch->Draw();
5499   }
5500   fNumberFitSuccess++;
5501 }
5502 //_____________________________________________________________________________
5503 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5504 {
5505   //
5506   // Fit methode for the gain factor
5507   //
5508
5509   fCurrentCoef[0]  = 0.0;
5510   fCurrentCoefE    = 0.0;
5511   Double_t chisqrl = 0.0;
5512   Double_t chisqrg = 0.0;
5513   Double_t chisqr  = 0.0;
5514   TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5515
5516   projch->Fit("landau","0",""
5517              ,(Double_t) mean/fBeginFitCharge
5518              ,projch->GetBinCenter(projch->GetNbinsX()));
5519   Double_t l3P0         = projch->GetFunction("landau")->GetParameter(0);
5520   Double_t l3P1         = projch->GetFunction("landau")->GetParameter(1);
5521   Double_t l3P2         = projch->GetFunction("landau")->GetParameter(2);
5522   chisqrl = projch->GetFunction("landau")->GetChisquare();
5523
5524   projch->Fit("gaus","0",""
5525               ,(Double_t) mean/fBeginFitCharge
5526               ,projch->GetBinCenter(projch->GetNbinsX()));
5527   Double_t g3P0         = projch->GetFunction("gaus")->GetParameter(0);
5528   Double_t g3P2         = projch->GetFunction("gaus")->GetParameter(2);
5529   chisqrg = projch->GetFunction("gaus")->GetChisquare();
5530
5531   fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5532   if (fDebugLevel != 1) {
5533     projch->Fit("fLandauGaus","0",""
5534                 ,(Double_t) mean/fBeginFitCharge
5535                 ,projch->GetBinCenter(projch->GetNbinsX()));
5536     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5537   }
5538   else  {
5539     TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5540     cp->cd();
5541     projch->Fit("fLandauGaus","+",""
5542                 ,(Double_t) mean/fBeginFitCharge
5543                 ,projch->GetBinCenter(projch->GetNbinsX()));
5544     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5545     projch->Draw();
5546     fLandauGaus->Draw("same");
5547   }
5548
5549   if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5550     //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5551     fNumberFitSuccess++;
5552     CalculChargeCoefMean(kTRUE);
5553     fCurrentCoef[0]  = projch->GetFunction("fLandauGaus")->GetParameter(1);
5554     fCurrentCoefE    = projch->GetFunction("fLandauGaus")->GetParError(1);
5555   }
5556   else {
5557     CalculChargeCoefMean(kFALSE);
5558     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5559   }
5560
5561   if (fDebugLevel != 1) {
5562     delete fLandauGaus;
5563   }
5564
5565 }
5566 //_____________________________________________________________________________
5567 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5568 {
5569   //
5570   // Fit methode for the gain factor more time consuming
5571   //
5572
5573
5574   //Some parameters to initialise
5575   Double_t widthLandau, widthGaus, mPV, integral;
5576   Double_t chisquarel = 0.0;
5577   Double_t chisquareg = 0.0;
5578   projch->Fit("landau","0M+",""
5579               ,(Double_t) mean/6
5580               ,projch->GetBinCenter(projch->GetNbinsX()));
5581   widthLandau  = projch->GetFunction("landau")->GetParameter(2);
5582   chisquarel = projch->GetFunction("landau")->GetChisquare();
5583   projch->Fit("gaus","0M+",""
5584               ,(Double_t) mean/6
5585               ,projch->GetBinCenter(projch->GetNbinsX()));
5586   widthGaus    = projch->GetFunction("gaus")->GetParameter(2);
5587   chisquareg = projch->GetFunction("gaus")->GetChisquare();
5588
5589   mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5590   integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5591
5592   // Setting fit range and start values
5593   Double_t fr[2];
5594   //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5595   //Double_t sv[4]   = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5596   Double_t sv[4]   = { widthLandau, mPV, integral, widthGaus};
5597   Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5598   Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5599   Double_t fp[4]   = { 1.0, 1.0, 1.0, 1.0 };
5600   Double_t fpe[4]  = { 1.0, 1.0, 1.0, 1.0 };
5601   fr[0]            = 0.3 * mean;
5602   fr[1]            = 3.0 * mean;
5603   fCurrentCoef[0]  = 0.0;
5604   fCurrentCoefE    = 0.0;
5605
5606   Double_t chisqr;
5607   Int_t    ndf;
5608   TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5609                                 ,&pllo[0],&plhi[0]
5610                                 ,&fp[0],&fpe[0]
5611                                 ,&chisqr,&ndf);
5612
5613   Double_t projchPeak;
5614   Double_t projchFWHM;
5615   LanGauPro(fp,projchPeak,projchFWHM);
5616
5617   if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5618     //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5619     fNumberFitSuccess++;
5620     CalculChargeCoefMean(kTRUE);
5621     fCurrentCoef[0]  = fp[1];
5622     fCurrentCoefE = fpe[1];
5623     //chargeCoefE2 = chisqr;
5624   }
5625   else {
5626     CalculChargeCoefMean(kFALSE);
5627     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5628   }
5629   if (fDebugLevel == 1) {
5630     AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5631     TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5632     cpy->cd();
5633     projch->Draw();
5634     fitsnr->Draw("same");
5635   }
5636   else {
5637     delete fitsnr;
5638   }
5639 }
5640 //_____________________________________________________________________________
5641 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
5642 {
5643   //
5644   // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5645   //
5646   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5647   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5648   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5649
5650   c4 = 0.0;
5651   c3 = 0.0;
5652   c2 = x0+x1+x2;
5653   c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5654   c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5655
5656 }
5657
5658 //_____________________________________________________________________________
5659 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
5660 {
5661   //
5662   // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5663   //
5664   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
5665   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
5666   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
5667   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
5668
5669   c4 = 0.0;
5670   c3 = x0+x1+x2+x3;
5671   c2 = -(x0*(x[1]+x[2]+x[3])
5672            +x1*(x[0]+x[2]+x[3])
5673            +x2*(x[0]+x[1]+x[3])
5674            +x3*(x[0]+x[1]+x[2]));
5675   c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
5676           +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
5677           +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
5678           +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
5679
5680   c0 = -(x0*x[1]*x[2]*x[3]
5681           +x1*x[0]*x[2]*x[3]
5682           +x2*x[0]*x[1]*x[3]
5683           +x3*x[0]*x[1]*x[2]);
5684
5685
5686 }
5687
5688 //_____________________________________________________________________________
5689 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
5690 {
5691   //
5692   // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5693   //
5694   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5695   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5696   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5697   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5698   Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5699
5700
5701   c4 = x0+x1+x2+x3+x4;
5702   c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
5703            +x1*(x[0]+x[2]+x[3]+x[4])
5704            +x2*(x[0]+x[1]+x[3]+x[4])
5705            +x3*(x[0]+x[1]+x[2]+x[4])
5706            +x4*(x[0]+x[1]+x[2]+x[3]));
5707   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])
5708           +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])
5709           +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])
5710           +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])
5711           +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]));
5712
5713   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])
5714           +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])
5715           +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])
5716           +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])
5717           +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]));
5718
5719   c0 = (x0*x[1]*x[2]*x[3]*x[4]
5720           +x1*x[0]*x[2]*x[3]*x[4]
5721           +x2*x[0]*x[1]*x[3]*x[4]
5722           +x3*x[0]*x[1]*x[2]*x[4]
5723           +x4*x[0]*x[1]*x[2]*x[3]);
5724
5725 }
5726 //_____________________________________________________________________________
5727 void AliTRDCalibraFit::NormierungCharge()
5728 {
5729   //
5730   // Normalisation of the gain factor resulting for the fits
5731   //
5732
5733   // Calcul of the mean of choosen method by fFitChargeNDB
5734   Double_t sum         = 0.0;
5735   //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5736   for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5737     Int_t    total    = 0;
5738     Int_t    detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5739     Float_t *coef     = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5740     //printf("detector %d coef[0] %f\n",detector,coef[0]);
5741     if (GetStack(detector) == 2) {
5742       total = 1728;
5743     }
5744     if (GetStack(detector) != 2) {
5745       total = 2304;
5746     }
5747     for (Int_t j = 0; j < total; j++) {
5748       if (coef[j] >= 0) {
5749         sum += coef[j];
5750       }
5751     }
5752   }
5753
5754   if (sum > 0) {
5755     fScaleFitFactor = fScaleFitFactor / sum;
5756   }
5757   else {
5758     fScaleFitFactor = 1.0;
5759   }
5760
5761   //methode de boeuf mais bon...
5762   Double_t scalefactor = fScaleFitFactor;
5763
5764   if(fDebugLevel > 1){
5765
5766     if ( !fDebugStreamer ) {
5767       //debug stream
5768       TDirectory *backup = gDirectory;
5769       fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5770       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5771     }
5772     (* fDebugStreamer) << "NormierungCharge"<<
5773       "scalefactor="<<scalefactor<<
5774       "\n";
5775     }
5776 }
5777 //_____________________________________________________________________________
5778 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5779 {
5780   //
5781   // Rebin of the 1D histo for the gain calibration if needed.
5782   // you have to choose fRebin, divider of fNumberBinCharge
5783   //
5784
5785  TAxis *xhist  = hist->GetXaxis();
5786  TH1I  *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5787                                         ,xhist->GetBinLowEdge(1)
5788                                         ,xhist->GetBinUpEdge(xhist->GetNbins()));
5789
5790  AliInfo(Form("fRebin: %d",fRebin));
5791  Int_t i = 1;
5792  for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5793    Double_t sum = 0.0;
5794    for (Int_t ji = i; ji < i+fRebin; ji++) {
5795      sum += hist->GetBinContent(ji);
5796    }
5797    sum = sum / fRebin;
5798    rehist->SetBinContent(k,sum);
5799    i += fRebin;
5800  }
5801
5802  return rehist;
5803
5804 }
5805
5806 //_____________________________________________________________________________
5807 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5808 {
5809   //
5810   // Rebin of the 1D histo for the gain calibration if needed
5811   // you have to choose fRebin divider of fNumberBinCharge
5812   //
5813
5814   TAxis *xhist  = hist->GetXaxis();
5815   TH1F  *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5816                                          ,xhist->GetBinLowEdge(1)
5817                                          ,xhist->GetBinUpEdge(xhist->GetNbins()));
5818
5819   AliInfo(Form("fRebin: %d",fRebin));
5820   Int_t i = 1;
5821   for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5822     Double_t sum = 0.0;
5823     for (Int_t ji = i; ji < i+fRebin; ji++) {
5824       sum += hist->GetBinContent(ji);
5825     }
5826     sum = sum/fRebin;
5827     rehist->SetBinContent(k,sum);
5828     i += fRebin;
5829   }
5830
5831   return rehist;
5832
5833 }
5834 //
5835 //____________Some basic geometry function_____________________________________
5836 //
5837
5838 //_____________________________________________________________________________
5839 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5840 {
5841   //
5842   // Reconstruct the plane number from the detector number
5843   //
5844
5845   return ((Int_t) (d % 6));
5846
5847 }
5848
5849 //_____________________________________________________________________________
5850 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5851 {
5852   //
5853   // Reconstruct the stack number from the detector number
5854   //
5855   const Int_t kNlayer = 6;
5856
5857   return ((Int_t) (d % 30) / kNlayer);
5858
5859 }
5860
5861 //_____________________________________________________________________________
5862 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5863 {
5864   //
5865   // Reconstruct the sector number from the detector number
5866   //
5867   Int_t fg = 30;
5868
5869   return ((Int_t) (d / fg));
5870
5871 }
5872
5873 //
5874 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5875 //
5876 //_______________________________________________________________________________
5877 void AliTRDCalibraFit::ResetVectorFit()
5878 {
5879   //
5880   // Reset the VectorFits
5881   //
5882
5883   fVectorFit.SetOwner();
5884   fVectorFit.Clear();
5885   fVectorFit2.SetOwner();
5886   fVectorFit2.Clear();
5887
5888 }
5889 //
5890 //____________Private Functions________________________________________________
5891 //
5892
5893 //_____________________________________________________________________________
5894 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5895 {
5896   //
5897   // Function for the fit
5898   //
5899
5900   //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5901
5902   //PARAMETERS FOR FIT PH
5903   // PASAv.4
5904   //fAsymmGauss->SetParameter(0,0.113755);
5905   //fAsymmGauss->SetParameter(1,0.350706);
5906   //fAsymmGauss->SetParameter(2,0.0604244);
5907   //fAsymmGauss->SetParameter(3,7.65596);
5908   //fAsymmGauss->SetParameter(4,1.00124);
5909   //fAsymmGauss->SetParameter(5,0.870597);  // No tail cancelation
5910
5911   Double_t xx = x[0];
5912
5913   if (xx < par[1]) {
5914     return par[5];
5915   }
5916
5917   Double_t dx       = 0.005;
5918   Double_t xs       = par[1];
5919   Double_t ss       = 0.0;
5920   Double_t paras[2] = { 0.0, 0.0 };
5921
5922   while (xs < xx) {
5923     if ((xs >= par[1]) &&
5924         (xs < (par[1]+par[2]))) {
5925       //fAsymmGauss->SetParameter(0,par[0]);
5926       //fAsymmGauss->SetParameter(1,xs);
5927       //ss += fAsymmGauss->Eval(xx);
5928       paras[0] = par[0];
5929       paras[1] = xs;
5930       ss += AsymmGauss(&xx,paras);
5931     }
5932     if ((xs >= (par[1]+par[2])) &&
5933         (xs <  (par[1]+par[2]+par[3]))) {
5934       //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5935       //fAsymmGauss->SetParameter(1,xs);
5936       //ss += fAsymmGauss->Eval(xx);
5937       paras[0] = par[0]*par[4];
5938       paras[1] = xs;
5939       ss += AsymmGauss(&xx,paras);
5940     }
5941     xs += dx;
5942   }
5943
5944   return ss + par[5];
5945
5946 }
5947
5948 //_____________________________________________________________________________
5949 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5950 {
5951   //
5952   // Function for the fit
5953   //
5954
5955   //par[0] = normalization
5956   //par[1] = mean
5957   //par[2] = sigma
5958   //norm0  = 1
5959   //par[3] = lambda0
5960   //par[4] = norm1
5961   //par[5] = lambda1
5962
5963   Double_t par1save = par[1];
5964   //Double_t par2save = par[2];
5965   Double_t par2save = 0.0604244;
5966   //Double_t par3save = par[3];
5967   Double_t par3save = 7.65596;
5968   //Double_t par5save = par[5];
5969   Double_t par5save = 0.870597;
5970   Double_t dx       = x[0] - par1save;
5971
5972   Double_t  sigma2  = par2save*par2save;
5973   Double_t  sqrt2   = TMath::Sqrt(2.0);
5974   Double_t  exp1    = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5975                                * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5976   Double_t  exp2    = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5977                                * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5978
5979   //return par[0]*(exp1+par[4]*exp2);
5980   return par[0] * (exp1 + 1.00124 * exp2);
5981
5982 }
5983
5984 //_____________________________________________________________________________
5985 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5986 {
5987   //
5988   // Sum Landau + Gaus with identical mean
5989   //
5990
5991   Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5992   //Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[4],par[5]);
5993   Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[1],par[4]);
5994   Double_t val       = valLandau + valGaus;
5995
5996   return val;
5997
5998 }
5999
6000 //_____________________________________________________________________________
6001 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6002 {
6003   //
6004   // Function for the fit
6005   //
6006   // Fit parameters:
6007   // par[0]=Width (scale) parameter of Landau density
6008   // par[1]=Most Probable (MP, location) parameter of Landau density
6009   // par[2]=Total area (integral -inf to inf, normalization constant)
6010   // par[3]=Width (sigma) of convoluted Gaussian function
6011   //
6012   // In the Landau distribution (represented by the CERNLIB approximation),
6013   // the maximum is located at x=-0.22278298 with the location parameter=0.
6014   // This shift is corrected within this function, so that the actual
6015   // maximum is identical to the MP parameter.
6016   //
6017
6018   // Numeric constants
6019   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
6020   Double_t mpshift  = -0.22278298;       // Landau maximum location
6021
6022   // Control constants
6023   Double_t np       = 100.0;             // Number of convolution steps
6024   Double_t sc       =   5.0;             // Convolution extends to +-sc Gaussian sigmas
6025
6026   // Variables
6027   Double_t xx;
6028   Double_t mpc;
6029   Double_t fland;
6030   Double_t sum = 0.0;
6031   Double_t xlow;
6032   Double_t xupp;
6033   Double_t step;
6034   Double_t i;
6035
6036   // MP shift correction
6037   mpc = par[1] - mpshift * par[0];
6038
6039   // Range of convolution integral
6040   xlow = x[0] - sc * par[3];
6041   xupp = x[0] + sc * par[3];
6042
6043   step = (xupp - xlow) / np;
6044
6045   // Convolution integral of Landau and Gaussian by sum
6046   for (i = 1.0; i <= np/2; i++) {
6047
6048     xx    = xlow + (i-.5) * step;
6049     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6050     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
6051
6052     xx    = xupp - (i-.5) * step;
6053     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6054     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
6055
6056   }
6057
6058   return (par[2] * step * sum * invsq2pi / par[3]);
6059
6060 }
6061 //_____________________________________________________________________________
6062 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6063                                       , const Double_t *parlimitslo, const Double_t *parlimitshi
6064                                       , Double_t *fitparams, Double_t *fiterrors
6065                                       , Double_t *chiSqr, Int_t *ndf) const
6066 {
6067   //
6068   // Function for the fit
6069   //
6070
6071   Int_t i;
6072   Char_t funname[100];
6073
6074   TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6075   if (ffitold) {
6076     delete ffitold;
6077   }
6078
6079   TF1 *ffit    = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6080   ffit->SetParameters(startvalues);
6081   ffit->SetParNames("Width","MP","Area","GSigma");
6082
6083   for (i = 0; i < 4; i++) {
6084     ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6085   }
6086
6087   his->Fit(funname,"RB0");                   // Fit within specified range, use ParLimits, do not plot
6088
6089   ffit->GetParameters(fitparams);            // Obtain fit parameters
6090   for (i = 0; i < 4; i++) {
6091     fiterrors[i] = ffit->GetParError(i);     // Obtain fit parameter errors
6092   }
6093   chiSqr[0] = ffit->GetChisquare();          // Obtain chi^2
6094   ndf[0]    = ffit->GetNDF();                // Obtain ndf
6095
6096   return (ffit);                             // Return fit function
6097
6098 }
6099
6100 //_____________________________________________________________________________
6101 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6102 {
6103   //
6104   // Function for the fit
6105   //
6106
6107   Double_t p;
6108   Double_t x;
6109   Double_t fy;
6110   Double_t fxr;
6111   Double_t fxl;
6112   Double_t step;
6113   Double_t l;
6114   Double_t lold;
6115
6116   Int_t    i        = 0;
6117   Int_t    maxcalls = 10000;
6118
6119   // Search for maximum
6120   p    = params[1] - 0.1 * params[0];
6121   step = 0.05 * params[0];
6122   lold = -2.0;
6123   l    = -1.0;
6124
6125   while ((l != lold) && (i < maxcalls)) {
6126     i++;
6127     lold = l;
6128     x    = p + step;
6129     l    = LanGauFun(&x,params);
6130     if (l < lold) {
6131       step = -step / 10.0;
6132     }
6133     p += step;
6134   }
6135
6136   if (i == maxcalls) {
6137     return (-1);
6138   }
6139   maxx = x;
6140   fy = l / 2.0;
6141
6142   // Search for right x location of fy
6143   p    = maxx + params[0];
6144   step = params[0];
6145   lold = -2.0;
6146   l    = -1e300;
6147   i    = 0;
6148
6149   while ( (l != lold) && (i < maxcalls) ) {
6150     i++;
6151
6152     lold = l;
6153     x = p + step;
6154     l = TMath::Abs(LanGauFun(&x,params) - fy);
6155
6156     if (l > lold)
6157       step = -step/10;
6158
6159     p += step;
6160   }
6161
6162   if (i == maxcalls)
6163     return (-2);
6164
6165   fxr = x;
6166
6167
6168   // Search for left x location of fy
6169
6170   p = maxx - 0.5 * params[0];
6171   step = -params[0];
6172   lold = -2.0;
6173   l    = -1.0e300;
6174   i    = 0;
6175
6176   while ((l != lold) && (i < maxcalls)) {
6177     i++;
6178     lold = l;
6179     x    = p + step;
6180     l    = TMath::Abs(LanGauFun(&x,params) - fy);
6181     if (l > lold) {
6182       step = -step / 10.0;
6183     }
6184     p += step;
6185   }
6186
6187   if (i == maxcalls) {
6188     return (-3);
6189   }
6190
6191   fxl  = x;
6192   fwhm = fxr - fxl;
6193
6194   return (0);
6195 }
6196 //_____________________________________________________________________________
6197 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6198 {
6199   //
6200   // Gaus with identical mean
6201   //
6202
6203   Double_t gauss   = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
6204
6205   return gauss;
6206
6207 }
6208
6209
6210