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   x[0]              = 0.0;
4320   x[1]              = 0.0;
4321   x[2]              = 0.0;
4322   x[3]              = 0.0;
4323   x[4]              = 0.0;
4324   y[0]              = 0.0;
4325   y[1]              = 0.0;
4326   y[2]              = 0.0;
4327   y[3]              = 0.0;
4328   y[4]              = 0.0;
4329   fPhd[0]           = 0.0;
4330   fPhd[1]           = 0.0;
4331   fPhd[2]           = 0.0;
4332   Int_t ju          = 0;
4333   fCurrentCoefE     = 0.0;
4334   fCurrentCoefE2    = 1.0;
4335   fCurrentCoef[0]   = 0.0;
4336   fCurrentCoef2[0]  = 0.0;
4337   TLine *line = new TLine();
4338   TF1 * polynome = 0x0;
4339   TF1 * polynomea = 0x0;
4340   TF1 * polynomeb = 0x0;
4341   Double_t *c = 0x0;
4342
4343   // Some variables
4344   TAxis   *xpph    = projPH->GetXaxis();
4345   Int_t    nbins   = xpph->GetNbins();
4346   Double_t lowedge = xpph->GetBinLowEdge(1);
4347   Double_t upedge  = xpph->GetBinUpEdge(xpph->GetNbins());
4348   Double_t widbins = (upedge - lowedge) / nbins;
4349   Double_t limit   = upedge + 0.5 * widbins;
4350
4351
4352   Bool_t put = kTRUE;
4353
4354   // Beginning of the signal
4355   TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
4356   for (Int_t k = 1; k <  projPH->GetNbinsX(); k++) {
4357     pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4358   }
4359
4360   binmax = (Int_t) pentea->GetMaximumBin();
4361
4362   Double_t minnn = 0.0;
4363   Double_t maxxx = 0.0;
4364
4365   Int_t kase = nbins-binmax;
4366
4367   switch(kase)
4368     {
4369     case 0:
4370       put = kFALSE;
4371       break;
4372     case 1:
4373       minnn = pentea->GetBinCenter(binmax-2);
4374       maxxx = pentea->GetBinCenter(binmax);
4375       x[0] = pentea->GetBinCenter(binmax-2);
4376       x[1] = pentea->GetBinCenter(binmax-1);
4377       x[2] = pentea->GetBinCenter(binmax);
4378       y[0] = pentea->GetBinContent(binmax-2);
4379       y[1] = pentea->GetBinContent(binmax-1);
4380       y[2] = pentea->GetBinContent(binmax);
4381       c = CalculPolynomeLagrange2(x,y);
4382       AliInfo("At the limit for beginning!");
4383       break;
4384     case 2:
4385       minnn = pentea->GetBinCenter(binmax-2);
4386       maxxx = pentea->GetBinCenter(binmax+1);
4387       x[0] = pentea->GetBinCenter(binmax-2);
4388       x[1] = pentea->GetBinCenter(binmax-1);
4389       x[2] = pentea->GetBinCenter(binmax);
4390       x[3] = pentea->GetBinCenter(binmax+1);
4391       y[0] = pentea->GetBinContent(binmax-2);
4392       y[1] = pentea->GetBinContent(binmax-1);
4393       y[2] = pentea->GetBinContent(binmax);
4394       y[3] = pentea->GetBinContent(binmax+1);
4395       c = CalculPolynomeLagrange3(x,y);
4396       break;
4397     default:
4398       switch(binmax){
4399       case 0:
4400         put = kFALSE;
4401         break;
4402       case 1:
4403         minnn = pentea->GetBinCenter(binmax);
4404         maxxx = pentea->GetBinCenter(binmax+2);
4405         x[0] = pentea->GetBinCenter(binmax);
4406         x[1] = pentea->GetBinCenter(binmax+1);
4407         x[2] = pentea->GetBinCenter(binmax+2);
4408         y[0] = pentea->GetBinContent(binmax);
4409         y[1] = pentea->GetBinContent(binmax+1);
4410         y[2] = pentea->GetBinContent(binmax+2);
4411         c = CalculPolynomeLagrange2(x,y);
4412         break;
4413       case 2:
4414         minnn = pentea->GetBinCenter(binmax-1);
4415         maxxx = pentea->GetBinCenter(binmax+2);
4416         x[0] = pentea->GetBinCenter(binmax-1);
4417         x[1] = pentea->GetBinCenter(binmax);
4418         x[2] = pentea->GetBinCenter(binmax+1);
4419         x[3] = pentea->GetBinCenter(binmax+2);
4420         y[0] = pentea->GetBinContent(binmax-1);
4421         y[1] = pentea->GetBinContent(binmax);
4422         y[2] = pentea->GetBinContent(binmax+1);
4423         y[3] = pentea->GetBinContent(binmax+2);
4424         c = CalculPolynomeLagrange3(x,y);
4425         break;
4426       default:
4427         minnn = pentea->GetBinCenter(binmax-2);
4428         maxxx = pentea->GetBinCenter(binmax+2);
4429         x[0] = pentea->GetBinCenter(binmax-2);
4430         x[1] = pentea->GetBinCenter(binmax-1);
4431         x[2] = pentea->GetBinCenter(binmax);
4432         x[3] = pentea->GetBinCenter(binmax+1);
4433         x[4] = pentea->GetBinCenter(binmax+2);
4434         y[0] = pentea->GetBinContent(binmax-2);
4435         y[1] = pentea->GetBinContent(binmax-1);
4436         y[2] = pentea->GetBinContent(binmax);
4437         y[3] = pentea->GetBinContent(binmax+1);
4438         y[4] = pentea->GetBinContent(binmax+2);
4439         c = CalculPolynomeLagrange4(x,y);
4440         break;
4441       }
4442       break;
4443     }
4444
4445
4446   if(put) {
4447     polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
4448     polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4449
4450     Double_t step = (maxxx-minnn)/10000;
4451     Double_t l = minnn;
4452     Double_t maxvalue = 0.0;
4453     Double_t placemaximum = minnn;
4454     for(Int_t o = 0; o < 10000; o++){
4455       if(o == 0) maxvalue = polynomeb->Eval(l);
4456       if(maxvalue < (polynomeb->Eval(l))){
4457         maxvalue = polynomeb->Eval(l);
4458         placemaximum = l;
4459       }
4460       l += step;
4461     }
4462     fPhd[0] = placemaximum;
4463   }
4464
4465   // Amplification region
4466   binmax = 0;
4467   ju     = 0;
4468   for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
4469     if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
4470       binmax = kbin;
4471       ju     = 1;
4472     }
4473   }
4474
4475   Double_t minn = 0.0;
4476   Double_t maxx = 0.0;
4477   x[0] = 0.0;
4478   x[1] = 0.0;
4479   x[2] = 0.0;
4480   x[3] = 0.0;
4481   x[4] = 0.0;
4482   y[0] = 0.0;
4483   y[1] = 0.0;
4484   y[2] = 0.0;
4485   y[3] = 0.0;
4486   y[4] = 0.0;
4487
4488   Int_t    kase1 = nbins - binmax;
4489
4490   //Determination of minn and maxx
4491   //case binmax = nbins
4492   //pol2
4493   switch(kase1)
4494     {
4495     case 0:
4496       minn = projPH->GetBinCenter(binmax-2);
4497       maxx = projPH->GetBinCenter(binmax);
4498       x[0] = projPH->GetBinCenter(binmax-2);
4499       x[1] = projPH->GetBinCenter(binmax-1);
4500       x[2] = projPH->GetBinCenter(binmax);
4501       y[0] = projPH->GetBinContent(binmax-2);
4502       y[1] = projPH->GetBinContent(binmax-1);
4503       y[2] = projPH->GetBinContent(binmax);
4504       c = CalculPolynomeLagrange2(x,y);
4505       //AliInfo("At the limit for the drift!");
4506       break;
4507     case 1:
4508       minn = projPH->GetBinCenter(binmax-2);
4509       maxx = projPH->GetBinCenter(binmax+1);
4510       x[0] = projPH->GetBinCenter(binmax-2);
4511       x[1] = projPH->GetBinCenter(binmax-1);
4512       x[2] = projPH->GetBinCenter(binmax);
4513       x[3] = projPH->GetBinCenter(binmax+1);
4514       y[0] = projPH->GetBinContent(binmax-2);
4515       y[1] = projPH->GetBinContent(binmax-1);
4516       y[2] = projPH->GetBinContent(binmax);
4517       y[3] = projPH->GetBinContent(binmax+1);
4518       c = CalculPolynomeLagrange3(x,y);
4519       break;
4520     default:
4521       switch(binmax)
4522         {
4523         case 0:
4524           put = kFALSE;
4525           break;
4526         case 1:
4527           minn = projPH->GetBinCenter(binmax);
4528           maxx = projPH->GetBinCenter(binmax+2);
4529           x[0] = projPH->GetBinCenter(binmax);
4530           x[1] = projPH->GetBinCenter(binmax+1);
4531           x[2] = projPH->GetBinCenter(binmax+2);
4532           y[0] = projPH->GetBinContent(binmax);
4533           y[1] = projPH->GetBinContent(binmax+1);
4534           y[2] = projPH->GetBinContent(binmax+2);
4535           c = CalculPolynomeLagrange2(x,y);
4536           break;
4537         case 2:
4538           minn = projPH->GetBinCenter(binmax-1);
4539           maxx = projPH->GetBinCenter(binmax+2);
4540           x[0] = projPH->GetBinCenter(binmax-1);
4541           x[1] = projPH->GetBinCenter(binmax);
4542           x[2] = projPH->GetBinCenter(binmax+1);
4543           x[3] = projPH->GetBinCenter(binmax+2);
4544           y[0] = projPH->GetBinContent(binmax-1);
4545           y[1] = projPH->GetBinContent(binmax);
4546           y[2] = projPH->GetBinContent(binmax+1);
4547           y[3] = projPH->GetBinContent(binmax+2);
4548           c = CalculPolynomeLagrange3(x,y);
4549           break;
4550         default:
4551           minn = projPH->GetBinCenter(binmax-2);
4552           maxx = projPH->GetBinCenter(binmax+2);
4553           x[0] = projPH->GetBinCenter(binmax-2);
4554           x[1] = projPH->GetBinCenter(binmax-1);
4555           x[2] = projPH->GetBinCenter(binmax);
4556           x[3] = projPH->GetBinCenter(binmax+1);
4557           x[4] = projPH->GetBinCenter(binmax+2);
4558           y[0] = projPH->GetBinContent(binmax-2);
4559           y[1] = projPH->GetBinContent(binmax-1);
4560           y[2] = projPH->GetBinContent(binmax);
4561           y[3] = projPH->GetBinContent(binmax+1);
4562           y[4] = projPH->GetBinContent(binmax+2);
4563           c = CalculPolynomeLagrange4(x,y);
4564           break;
4565         }
4566       break;
4567     }
4568
4569   if(put) {
4570     polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
4571     polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4572
4573     Double_t step = (maxx-minn)/1000;
4574     Double_t l = minn;
4575     Double_t maxvalue = 0.0;
4576     Double_t placemaximum = minn;
4577     for(Int_t o = 0; o < 1000; o++){
4578       if(o == 0) maxvalue = polynomea->Eval(l);
4579       if(maxvalue < (polynomea->Eval(l))){
4580         maxvalue = polynomea->Eval(l);
4581         placemaximum = l;
4582       }
4583       l += step;
4584     }
4585     fPhd[1] = placemaximum;
4586   }
4587
4588   // Drift region
4589   TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
4590   for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k <  projPH->GetNbinsX(); k++) {
4591     pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
4592   }
4593   binmin = 0;
4594   if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4595
4596   //should not happen
4597   if (binmin <= 1) {
4598     binmin = 2;
4599     put = 1;
4600     AliInfo("Put the binmax from 1 to 2 to enable the fit");
4601   }
4602
4603   //check
4604   if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
4605     AliInfo("Too many fluctuations at the end!");
4606     put = kFALSE;
4607   }
4608   if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
4609     AliInfo("Too many fluctuations at the end!");
4610     put = kFALSE;
4611   }
4612   if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
4613     AliInfo("No entries for the next bin!");
4614     pente->SetBinContent(binmin,0);
4615     if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
4616   }
4617
4618
4619   x[0] = 0.0;
4620   x[1] = 0.0;
4621   x[2] = 0.0;
4622   x[3] = 0.0;
4623   x[4] = 0.0;
4624   y[0] = 0.0;
4625   y[1] = 0.0;
4626   y[2] = 0.0;
4627   y[3] = 0.0;
4628   y[4] = 0.0;
4629   Double_t min = 0.0;
4630   Double_t max = 0.0;
4631   Bool_t case1 = kFALSE;
4632   Bool_t case2 = kFALSE;
4633   Bool_t case4 = kFALSE;
4634
4635   //Determination of min and max
4636   //case binmin <= nbins-3
4637   //pol4 case 3
4638   if((binmin <= (nbins-3)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4639     min = pente->GetBinCenter(binmin-2);
4640     max = pente->GetBinCenter(binmin+2);
4641     x[0] = pente->GetBinCenter(binmin-2);
4642     x[1] = pente->GetBinCenter(binmin-1);
4643     x[2] = pente->GetBinCenter(binmin);
4644     x[3] = pente->GetBinCenter(binmin+1);
4645     x[4] = pente->GetBinCenter(binmin+2);
4646     y[0] = pente->GetBinContent(binmin-2);
4647     y[1] = pente->GetBinContent(binmin-1);
4648     y[2] = pente->GetBinContent(binmin);
4649     y[3] = pente->GetBinContent(binmin+1);
4650     y[4] = pente->GetBinContent(binmin+2);
4651     //Calcul the polynome de Lagrange
4652     c = CalculPolynomeLagrange4(x,y);
4653     //richtung +/-
4654     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4655        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4656       //AliInfo("polynome 4 false 1");
4657       put = kFALSE;
4658     }
4659     if(((binmin+3) <= (nbins-1)) &&
4660        (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
4661        ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
4662        (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
4663       AliInfo("polynome 4 false 2");
4664       put = kFALSE;
4665     }
4666     // poly 3
4667     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
4668        (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) {
4669       //AliInfo("polynome 4 case 1");
4670       case1 = kTRUE;
4671     }
4672     if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) &&
4673        (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4674       //AliInfo("polynome 4 case 4");
4675       case4 = kTRUE;
4676     }
4677
4678   }
4679   //case binmin = nbins-2
4680   //pol3 case 1
4681   if(((binmin == (nbins-2)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4682      (case1)){
4683     min = pente->GetBinCenter(binmin-2);
4684     max = pente->GetBinCenter(binmin+1);
4685     x[0] = pente->GetBinCenter(binmin-2);
4686     x[1] = pente->GetBinCenter(binmin-1);
4687     x[2] = pente->GetBinCenter(binmin);
4688     x[3] = pente->GetBinCenter(binmin+1);
4689     y[0] = pente->GetBinContent(binmin-2);
4690     y[1] = pente->GetBinContent(binmin-1);
4691     y[2] = pente->GetBinContent(binmin);
4692     y[3] = pente->GetBinContent(binmin+1);
4693     //Calcul the polynome de Lagrange
4694     c = CalculPolynomeLagrange3(x,y);
4695     //richtung +: nothing
4696     //richtung -
4697     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4698       //AliInfo("polynome 3- case 2");
4699       case2 = kTRUE;
4700     }
4701   }
4702   //pol3 case 4
4703   if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4704      (case4)){
4705     min = pente->GetBinCenter(binmin-1);
4706     max = pente->GetBinCenter(binmin+2);
4707     x[0] = pente->GetBinCenter(binmin-1);
4708     x[1] = pente->GetBinCenter(binmin);
4709     x[2] = pente->GetBinCenter(binmin+1);
4710     x[3] = pente->GetBinCenter(binmin+2);
4711     y[0] = pente->GetBinContent(binmin-1);
4712     y[1] = pente->GetBinContent(binmin);
4713     y[2] = pente->GetBinContent(binmin+1);
4714     y[3] = pente->GetBinContent(binmin+2);
4715     //Calcul the polynome de Lagrange
4716     c = CalculPolynomeLagrange3(x,y);
4717     //richtung +
4718     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4719       //AliInfo("polynome 3+ case 2");
4720       case2 = kTRUE;
4721     }
4722   }
4723   //pol2 case 5
4724   if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){
4725     min = pente->GetBinCenter(binmin);
4726     max = pente->GetBinCenter(binmin+2);
4727     x[0] = pente->GetBinCenter(binmin);
4728     x[1] = pente->GetBinCenter(binmin+1);
4729     x[2] = pente->GetBinCenter(binmin+2);
4730     y[0] = pente->GetBinContent(binmin);
4731     y[1] = pente->GetBinContent(binmin+1);
4732     y[2] = pente->GetBinContent(binmin+2);
4733     //Calcul the polynome de Lagrange
4734     c = CalculPolynomeLagrange2(x,y);
4735     //richtung +
4736     if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
4737       //AliInfo("polynome 2+ false");
4738       put = kFALSE;
4739     }
4740   }
4741   //pol2 case 2
4742   if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) ||
4743      (case2)){
4744     min = pente->GetBinCenter(binmin-1);
4745     max = pente->GetBinCenter(binmin+1);
4746     x[0] = pente->GetBinCenter(binmin-1);
4747     x[1] = pente->GetBinCenter(binmin);
4748     x[2] = pente->GetBinCenter(binmin+1);
4749     y[0] = pente->GetBinContent(binmin-1);
4750     y[1] = pente->GetBinContent(binmin);
4751     y[2] = pente->GetBinContent(binmin+1);
4752     //Calcul the polynome de Lagrange
4753     c = CalculPolynomeLagrange2(x,y);
4754     //richtung +: nothing
4755     //richtung -: nothing
4756   }
4757   //case binmin = nbins-1
4758   //pol2 case 0
4759   if((binmin == (nbins-1)) && ((binmin-2) >= TMath::Min(binmax+4, projPH->GetNbinsX()))){
4760     min = pente->GetBinCenter(binmin-2);
4761     max = pente->GetBinCenter(binmin);
4762     x[0] = pente->GetBinCenter(binmin-2);
4763     x[1] = pente->GetBinCenter(binmin-1);
4764     x[2] = pente->GetBinCenter(binmin);
4765     y[0] = pente->GetBinContent(binmin-2);
4766     y[1] = pente->GetBinContent(binmin-1);
4767     y[2] = pente->GetBinContent(binmin);
4768     //Calcul the polynome de Lagrange
4769     c = CalculPolynomeLagrange2(x,y);
4770     //AliInfo("At the limit for the drift!");
4771     //fluctuation too big!
4772     //richtung +: nothing
4773     //richtung -
4774     if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
4775       //AliInfo("polynome 2- false ");
4776       put = kFALSE;
4777     }
4778   }
4779   if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
4780     put = kFALSE;
4781     AliInfo("At the limit for the drift and not usable!");
4782   }
4783
4784   //pass
4785   if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
4786     put = kFALSE;
4787     AliInfo("For the drift...problem!");
4788   }
4789   //pass but should not happen
4790   if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
4791     put = kFALSE;
4792     AliInfo("For the drift...problem!");
4793   }
4794
4795   if(put) {
4796     polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
4797     polynome->SetParameters(c[0],c[1],c[2],c[3],c[4]);
4798     //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
4799     Double_t step = (max-min)/1000;
4800     Double_t l = min;
4801     Double_t minvalue = 0.0;
4802     Double_t placeminimum = min;
4803     for(Int_t o = 0; o < 1000; o++){
4804       if(o == 0) minvalue = polynome->Eval(l);
4805       if(minvalue > (polynome->Eval(l))){
4806         minvalue = polynome->Eval(l);
4807         placeminimum = l;
4808       }
4809       l += step;
4810     }
4811     fPhd[2] = placeminimum;
4812   }
4813   //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2);
4814   if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0;
4815   if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0;
4816
4817   Float_t fPhdt0  = 0.0;
4818   Float_t t0Shift = 0.0;
4819   if(fTakeTheMaxPH) {
4820     fPhdt0 = fPhd[1];
4821     t0Shift = fT0Shift1;
4822   }
4823   else {
4824     fPhdt0 = fPhd[0];
4825     t0Shift = fT0Shift0;
4826   }
4827
4828   if ((fPhd[2] > fPhd[0]) &&
4829       (fPhd[2] > fPhd[1]) &&
4830       (fPhd[1] > fPhd[0]) &&
4831       (put)) {
4832     fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
4833     if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] =  -TMath::Abs(fCurrentCoef[1]);
4834     else fNumberFitSuccess++;
4835     if (fPhdt0 >= 0.0) {
4836       fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4837       if (fCurrentCoef2[0] < -1.0) {
4838         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4839       }
4840     }
4841     else {
4842       fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4843     }
4844   }
4845   else {
4846     //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
4847     fCurrentCoef[0]      = -TMath::Abs(fCurrentCoef[1]);
4848
4849     if((fPhd[1] > fPhd[0]) &&
4850        (put)) {
4851       if (fPhdt0 >= 0.0) {
4852         fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
4853         if (fCurrentCoef2[0] < -1.0) {
4854           fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4855         }
4856       }
4857       else {
4858         fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4859       }
4860     }
4861     else{
4862       fCurrentCoef2[0]     = fCurrentCoef2[1] + 100.0;
4863       //printf("Fit failed!\n");
4864     }
4865   }
4866
4867   if (fDebugLevel == 1) {
4868     TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
4869     cpentei->cd();
4870     projPH->Draw();
4871     line->SetLineColor(2);
4872     line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
4873     line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
4874     line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
4875     AliInfo(Form("fPhd[0] (beginning of the signal): %f"                  ,(Float_t) fPhd[0]));
4876     AliInfo(Form("fPhd[1] (end of the amplification region): %f"          ,(Float_t) fPhd[1]));
4877     AliInfo(Form("fPhd[2] (end of the drift region): %f"                  ,(Float_t) fPhd[2]));
4878     AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
4879     TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
4880     cpentei2->cd();
4881     pentea->Draw();
4882     TCanvas *cpentei3 = new TCanvas("cpentei3","cpentei3",50,50,600,800);
4883     cpentei3->cd();
4884     pente->Draw();
4885   }
4886   else {
4887     if(pentea) delete pentea;
4888     if(pente) delete pente;
4889     if(polynome) delete polynome;
4890     if(polynomea) delete polynomea;
4891     if(polynomeb) delete polynomeb;
4892     if(x) delete [] x;
4893     if(y) delete [] y;
4894     if(c) delete [] c;
4895     if(line) delete line;
4896
4897   }
4898
4899   //Provisoire
4900   //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
4901   //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4902
4903   projPH->SetDirectory(0);
4904
4905 }
4906
4907 //_____________________________________________________________________________
4908 void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect)
4909 {
4910   //
4911   // Fit methode for the drift velocity
4912   //
4913
4914   // Constants
4915   const Float_t kDrWidth = AliTRDgeometry::DrThick();
4916
4917   // Some variables
4918   TAxis   *xpph   = projPH->GetXaxis();
4919   Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
4920
4921   TF1 *fPH = new TF1("fPH",AliTRDCalibraFit::PH,-0.05,3.2,6);
4922   fPH->SetParameter(0,0.469);     // Scaling
4923   fPH->SetParameter(1,0.18);      // Start
4924   fPH->SetParameter(2,0.0857325); // AR
4925   fPH->SetParameter(3,1.89);      // DR
4926   fPH->SetParameter(4,0.08);      // QA/QD
4927   fPH->SetParameter(5,0.0);       // Baseline
4928
4929   TLine *line = new TLine();
4930
4931   fCurrentCoef[0]     = 0.0;
4932   fCurrentCoef2[0]    = 0.0;
4933   fCurrentCoefE       = 0.0;
4934   fCurrentCoefE2      = 0.0;
4935
4936   if (idect%fFitPHPeriode == 0) {
4937
4938     AliInfo(Form("The detector %d will be fitted",idect));
4939     fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
4940     fPH->SetParameter(1,fPhd[0] - 0.1);                                                                 // Start
4941     fPH->SetParameter(2,fPhd[1] - fPhd[0]);                                                             // AR
4942     fPH->SetParameter(3,fPhd[2] - fPhd[1]);                                                             // DR
4943     fPH->SetParameter(4,0.225);                                                                         // QA/QD
4944     fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
4945
4946     if (fDebugLevel != 1) {
4947       projPH->Fit(fPH,"0M","",0.0,upedge);
4948     }
4949     else {
4950       TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
4951       cpente->cd();
4952       projPH->Fit(fPH,"M+","",0.0,upedge);
4953       projPH->Draw("E0");
4954       line->SetLineColor(4);
4955       line->DrawLine(fPH->GetParameter(1)
4956                     ,0
4957                     ,fPH->GetParameter(1)
4958                     ,projPH->GetMaximum());
4959       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
4960                     ,0
4961                     ,fPH->GetParameter(1)+fPH->GetParameter(2)
4962                     ,projPH->GetMaximum());
4963       line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4964                     ,0
4965                     ,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
4966                     ,projPH->GetMaximum());
4967     }
4968
4969     if (fPH->GetParameter(3) != 0) {
4970       fNumberFitSuccess++;
4971       fCurrentCoef[0]    = kDrWidth / (fPH->GetParameter(3));
4972       fCurrentCoefE      = (fPH->GetParError(3)/fPH->GetParameter(3))*fCurrentCoef[0];
4973       fCurrentCoef2[0]   = fPH->GetParameter(1);
4974       fCurrentCoefE2     = fPH->GetParError(1);
4975     }
4976     else {
4977       fCurrentCoef[0]     = -TMath::Abs(fCurrentCoef[1]);
4978       fCurrentCoef2[0]    = fCurrentCoef2[1] + 100.0;
4979     }
4980
4981   }
4982   else {
4983
4984     // Put the default value
4985     fCurrentCoef[0]  = -TMath::Abs(fCurrentCoef[1]);
4986     fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
4987   }
4988
4989   if (fDebugLevel != 1) {
4990     delete fPH;
4991   }
4992
4993 }
4994 //_____________________________________________________________________________
4995 Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nBins, Float_t xMin, Float_t xMax)
4996 {
4997   //
4998   // Fit methode for the sigma of the pad response function
4999   //
5000
5001   TVectorD param(3);
5002
5003   fCurrentCoef[0]  = 0.0;
5004   fCurrentCoefE = 0.0;
5005
5006   Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,&param);
5007
5008   if(TMath::Abs(ret+4) <= 0.000000001){
5009     fCurrentCoef[0] = -fCurrentCoef[1];
5010     return kFALSE;
5011   }
5012   else {
5013     fNumberFitSuccess++;
5014     fCurrentCoef[0] = param[2];
5015     fCurrentCoefE   = ret;
5016     return kTRUE;
5017   }
5018 }
5019 //_____________________________________________________________________________
5020 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)
5021 {
5022   //
5023   // Fit methode for the sigma of the pad response function
5024   //
5025
5026   //We should have at least 3 points
5027   if(nBins <=3) return -4.0;
5028
5029   TLinearFitter fitter(3,"pol2");
5030   fitter.StoreData(kFALSE);
5031   fitter.ClearPoints();
5032   TVectorD  par(3);
5033   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
5034   Float_t entries = 0;
5035   Int_t   nbbinwithentries = 0;
5036   for (Int_t i=0; i<nBins; i++){
5037     entries+=arraye[i];
5038     if(arraye[i] > 15) nbbinwithentries++;
5039     //printf("entries for i %d: %f\n",i,arraye[i]);
5040   }
5041   if ((entries<700) || (nbbinwithentries < ((Int_t)(nBins/2)))) return -4;
5042   //printf("entries %f\n",entries);
5043   //printf("nbbinwithentries %d\n",nbbinwithentries);
5044
5045   Int_t npoints=0;
5046   Float_t errorm = 0.0;
5047   Float_t errorn = 0.0;
5048   Float_t error  = 0.0;
5049
5050   //
5051   for (Int_t ibin=0;ibin<nBins; ibin++){
5052       Float_t entriesI = arraye[ibin];
5053       Float_t valueI   = arraym[ibin];
5054       Double_t xcenter = 0.0;
5055       Float_t  val     = 0.0;
5056       if ((entriesI>15) && (valueI>0.0)){
5057         xcenter = xMin+(ibin+0.5)*binWidth;
5058         errorm   = 0.0;
5059         errorn   = 0.0;
5060         error    = 0.0;
5061         if(!bError){
5062           if((valueI + 0.01) > 0.0) errorm = TMath::Log((valueI + 0.01)/valueI);
5063           if((valueI - 0.01) > 0.0) errorn = TMath::Log((valueI - 0.01)/valueI);
5064           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5065         }
5066         else{
5067           if((valueI + arrayme[ibin]) > 0.0) errorm = TMath::Log((valueI + arrayme[ibin])/valueI);
5068           if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
5069           error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
5070         }
5071         if(TMath::Abs(error) < 0.000000001) continue;
5072         val      = TMath::Log(Float_t(valueI));
5073         fitter.AddPoint(&xcenter,val,error);
5074         npoints++;
5075       }
5076
5077       if(fDebugLevel > 1){
5078
5079       if ( !fDebugStreamer ) {
5080         //debug stream
5081         TDirectory *backup = gDirectory;
5082         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5083         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5084       }
5085
5086       Int_t    detector     = fCountDet;
5087       Int_t    layer        = GetLayer(fCountDet);
5088       Int_t    group        = ibin;
5089
5090       (* fDebugStreamer) << "FitGausMIFill"<<
5091         "detector="<<detector<<
5092         "layer="<<layer<<
5093         "nbins="<<nBins<<
5094         "group="<<group<<
5095         "entriesI="<<entriesI<<
5096         "valueI="<<valueI<<
5097         "val="<<val<<
5098         "xcenter="<<xcenter<<
5099         "errorm="<<errorm<<
5100         "errorn="<<errorn<<
5101         "error="<<error<<
5102         "bError="<<bError<<
5103         "\n";
5104     }
5105
5106   }
5107
5108   if(npoints <=3) return -4.0;
5109
5110   Double_t chi2 = 0;
5111   if (npoints>3){
5112     fitter.Eval();
5113     fitter.GetParameters(par);
5114     chi2 = fitter.GetChisquare()/Float_t(npoints);
5115
5116
5117     if (!param)  param  = new TVectorD(3);
5118     if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
5119     Double_t  x      = TMath::Sqrt(TMath::Abs(-2*par[2]));
5120     Double_t deltax = (fitter.GetParError(2))/x;
5121     Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
5122     chi2 = errorparam2;
5123
5124     (*param)[1] = par[1]/(-2.*par[2]);
5125     (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
5126     Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
5127     if ( lnparam0>307 ) return -4;
5128     (*param)[0] = TMath::Exp(lnparam0);
5129
5130     if(fDebugLevel > 1){
5131
5132       if ( !fDebugStreamer ) {
5133         //debug stream
5134         TDirectory *backup = gDirectory;
5135         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5136         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5137       }
5138
5139       Int_t    detector     = fCountDet;
5140       Int_t    layer        = GetLayer(fCountDet);
5141
5142
5143       (* fDebugStreamer) << "FitGausMIFit"<<
5144         "detector="<<detector<<
5145         "layer="<<layer<<
5146         "nbins="<<nBins<<
5147         "errorsigma="<<chi2<<
5148         "mean="<<(*param)[1]<<
5149         "sigma="<<(*param)[2]<<
5150         "constant="<<(*param)[0]<<
5151         "\n";
5152     }
5153   }
5154
5155   if((chi2/(*param)[2]) > 0.1){
5156     if(bError){
5157       chi2 = FitGausMI(arraye,arraym,arrayme,nBins,xMin,xMax,param,kFALSE);
5158     }
5159     else return -4.0;
5160   }
5161
5162   if(fDebugLevel == 1){
5163     TString name("PRF");
5164     name += (Int_t)xMin;
5165     name += (Int_t)xMax;
5166     TCanvas *c1 = new TCanvas((const char *)name,(const char *)name,50,50,600,800);
5167     c1->cd();
5168     name += "histo";
5169     TH1F *histo = new TH1F((const char *)name,(const char *)name,nBins,xMin,xMax);
5170     for(Int_t k = 0; k < nBins; k++){
5171       histo->SetBinContent(k+1,arraym[k]);
5172       histo->SetBinError(k+1,arrayme[k]);
5173     }
5174     histo->Draw();
5175     name += "functionf";
5176     TF1 *f1= new TF1((const char*)name,"[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
5177     f1->SetParameter(0, (*param)[0]);
5178     f1->SetParameter(1, (*param)[1]);
5179     f1->SetParameter(2, (*param)[2]);
5180     f1->Draw("same");
5181   }
5182
5183
5184   return chi2;
5185
5186 }
5187 //_____________________________________________________________________________
5188 void AliTRDCalibraFit::FitPRF(TH1 *projPRF)
5189 {
5190   //
5191   // Fit methode for the sigma of the pad response function
5192   //
5193
5194   fCurrentCoef[0]  = 0.0;
5195   fCurrentCoefE = 0.0;
5196
5197   if (fDebugLevel != 1) {
5198     projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
5199   }
5200   else {
5201     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5202     cfit->cd();
5203     projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
5204     projPRF->Draw();
5205   }
5206   fCurrentCoef[0]  = projPRF->GetFunction("gaus")->GetParameter(2);
5207   fCurrentCoefE = projPRF->GetFunction("gaus")->GetParError(2);
5208   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5209   else {
5210     fNumberFitSuccess++;
5211   }
5212 }
5213 //_____________________________________________________________________________
5214 void AliTRDCalibraFit::RmsPRF(TH1 *projPRF)
5215 {
5216   //
5217   // Fit methode for the sigma of the pad response function
5218   //
5219   fCurrentCoef[0]   = 0.0;
5220   fCurrentCoefE  = 0.0;
5221   if (fDebugLevel == 1) {
5222     TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
5223     cfit->cd();
5224     projPRF->Draw();
5225   }
5226   fCurrentCoef[0] = projPRF->GetRMS();
5227   if(fCurrentCoef[0] <= 0.0) fCurrentCoef[0] = -fCurrentCoef[1];
5228   else {
5229     fNumberFitSuccess++;
5230   }
5231 }
5232 //_____________________________________________________________________________
5233 void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t *arrayme, Int_t nbg, Int_t nybins)
5234 {
5235   //
5236   // Fit methode for the sigma of the pad response function with 2*nbg tan bins
5237   //
5238
5239   TLinearFitter linearfitter = TLinearFitter(3,"pol2");
5240
5241
5242   Int_t   nbins    = (Int_t)(nybins/(2*nbg));
5243   Float_t lowedge  = -3.0*nbg;
5244   Float_t upedge   = lowedge + 3.0;
5245   Int_t   offset   = 0;
5246   Int_t   npoints  = 0;
5247   Double_t xvalues = -0.2*nbg+0.1;
5248   Double_t y       = 0.0;
5249   Int_t   total    = 2*nbg;
5250
5251
5252   for(Int_t k = 0; k < total; k++){
5253     if(FitPRFGausMI(arraye+offset, arraym+offset, arrayme+offset, nbins, lowedge, upedge)){
5254       npoints++;
5255       y = fCurrentCoef[0]*fCurrentCoef[0];
5256       linearfitter.AddPoint(&xvalues,y,2*fCurrentCoefE*fCurrentCoef[0]);
5257     }
5258
5259     if(fDebugLevel > 1){
5260
5261       if ( !fDebugStreamer ) {
5262         //debug stream
5263         TDirectory *backup = gDirectory;
5264         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5265         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5266       }
5267
5268       Int_t    detector     = fCountDet;
5269       Int_t    layer        = GetLayer(fCountDet);
5270       Int_t    nbtotal      = total;
5271       Int_t    group        = k;
5272       Float_t  low          = lowedge;
5273       Float_t  up           = upedge;
5274       Float_t  tnp          = xvalues;
5275       Float_t  wid          = fCurrentCoef[0];
5276       Float_t  widfE        = fCurrentCoefE;
5277
5278       (* fDebugStreamer) << "FitTnpRange0"<<
5279         "detector="<<detector<<
5280         "layer="<<layer<<
5281         "nbtotal="<<nbtotal<<
5282         "group="<<group<<
5283         "low="<<low<<
5284         "up="<<up<<
5285         "offset="<<offset<<
5286         "tnp="<<tnp<<
5287         "wid="<<wid<<
5288         "widfE="<<widfE<<
5289         "\n";
5290     }
5291
5292     offset  += nbins;
5293     lowedge += 3.0;
5294     upedge  += 3.0;
5295     xvalues += 0.2;
5296
5297   }
5298
5299   fCurrentCoefE = 0.0;
5300   fCurrentCoef[0] = 0.0;
5301
5302   //printf("npoints\n",npoints);
5303
5304   if(npoints < 3){
5305     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5306   }
5307   else{
5308
5309     TVectorD pars0;
5310     linearfitter.Eval();
5311     linearfitter.GetParameters(pars0);
5312     Double_t pointError0  =  TMath::Sqrt(linearfitter.GetChisquare()/npoints);
5313     Double_t errorsx0     =  linearfitter.GetParError(2)*pointError0;
5314     Double_t min0         = 0.0;
5315     Double_t ermin0       = 0.0;
5316     //Double_t prfe0      = 0.0;
5317     Double_t prf0         = 0.0;
5318     if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
5319       min0 = -pars0[1]/(2*pars0[2]);
5320       ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
5321       prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
5322       if(prf0 > 0.0) {
5323         /*
5324           prfe0 = linearfitter->GetParError(0)*pointError0
5325           +(linearfitter->GetParError(1)*pointError0/pars0[1]+ermin0/min0)*pars0[1]*min0
5326           +(linearfitter->GetParError(2)*pointError0/pars0[2]+2*ermin0/min0)*pars0[2]*min0*min0;
5327           prfe0 = prfe0/(2*TMath::Sqrt(prf0));
5328           fCurrentCoefE   = (Float_t) prfe0;
5329         */
5330         fCurrentCoef[0] = (Float_t) TMath::Sqrt(TMath::Abs(prf0));
5331       }
5332       else{
5333         fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5334       }
5335     }
5336     else {
5337       fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5338     }
5339
5340     if(fDebugLevel > 1){
5341
5342       if ( !fDebugStreamer ) {
5343         //debug stream
5344         TDirectory *backup = gDirectory;
5345         fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root");
5346         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5347       }
5348
5349       Int_t    detector     = fCountDet;
5350       Int_t    layer        = GetLayer(fCountDet);
5351       Int_t    nbtotal      = total;
5352       Double_t colsize[6]   = {0.635,0.665,0.695,0.725,0.755,0.785};
5353       Double_t sigmax       = TMath::Sqrt(TMath::Abs(pars0[2]))*10000*colsize[layer];
5354
5355       (* fDebugStreamer) << "FitTnpRange1"<<
5356         "detector="<<detector<<
5357         "layer="<<layer<<
5358         "nbtotal="<<nbtotal<<
5359         "par0="<<pars0[0]<<
5360         "par1="<<pars0[1]<<
5361         "par2="<<pars0[2]<<
5362         "npoints="<<npoints<<
5363         "sigmax="<<sigmax<<
5364         "tan="<<min0<<
5365         "sigmaprf="<<fCurrentCoef[0]<<
5366         "sigprf="<<fCurrentCoef[1]<<
5367         "\n";
5368     }
5369
5370   }
5371
5372 }
5373 //_____________________________________________________________________________
5374 void AliTRDCalibraFit::FitMean(TH1 *projch, Double_t nentries, Double_t mean)
5375 {
5376   //
5377   // Only mean methode for the gain factor
5378   //
5379
5380   fCurrentCoef[0] = mean;
5381   fCurrentCoefE   = 0.0;
5382   if(nentries > 0) fCurrentCoefE = projch->GetRMS()/TMath::Sqrt(nentries);
5383   if (fDebugLevel == 1) {
5384     TCanvas *cpmean = new TCanvas("cpmean","cpmean",50,50,600,800);
5385     cpmean->cd();
5386     projch->Draw();
5387   }
5388   CalculChargeCoefMean(kTRUE);
5389   fNumberFitSuccess++;
5390 }
5391 //_____________________________________________________________________________
5392 void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries)
5393 {
5394   //
5395   // mean w methode for the gain factor
5396   //
5397
5398   //Number of bins
5399   Int_t nybins = projch->GetNbinsX();
5400
5401   //The weight function
5402   Double_t a = 0.00228515;
5403   Double_t b = -0.00231487;
5404   Double_t c = 0.00044298;
5405   Double_t d = -0.00379239;
5406   Double_t e = 0.00338349;
5407
5408 //   0 |0.00228515
5409 //    1 |-0.00231487
5410 //    2 |0.00044298
5411 //    3 |-0.00379239
5412 //    4 |0.00338349
5413
5414
5415
5416   //A arbitrary error for the moment
5417   fCurrentCoefE = 0.0;
5418   fCurrentCoef[0] = 0.0;
5419
5420   //Calcul
5421   Double_t sumw = 0.0;
5422   Double_t sum = 0.0;
5423   Float_t sumAll   = (Float_t) nentries;
5424   Int_t sumCurrent = 0;
5425   for(Int_t k = 0; k <nybins; k++){
5426     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5427     if (fraction>0.95) break;
5428     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5429       e*fraction*fraction*fraction*fraction;
5430     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5431     sum  += weight*projch->GetBinContent(k+1);
5432     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5433     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5434   }
5435   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5436
5437   if (fDebugLevel == 1) {
5438     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5439     cpmeanw->cd();
5440     projch->Draw();
5441   }
5442   fNumberFitSuccess++;
5443   CalculChargeCoefMean(kTRUE);
5444 }
5445 //_____________________________________________________________________________
5446 void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll)
5447 {
5448   //
5449   // mean w methode for the gain factor
5450   //
5451
5452   //Number of bins
5453   Int_t nybins = projch->GetNbinsX();
5454
5455   //The weight function
5456   Double_t a = 0.00228515;
5457   Double_t b = -0.00231487;
5458   Double_t c = 0.00044298;
5459   Double_t d = -0.00379239;
5460   Double_t e = 0.00338349;
5461
5462 //   0 |0.00228515
5463 //    1 |-0.00231487
5464 //    2 |0.00044298
5465 //    3 |-0.00379239
5466 //    4 |0.00338349
5467
5468
5469
5470   //A arbitrary error for the moment
5471   fCurrentCoefE = 0.0;
5472   fCurrentCoef[0] = 0.0;
5473
5474   //Calcul
5475   Double_t sumw = 0.0;
5476   Double_t sum = 0.0;
5477   Int_t sumCurrent = 0;
5478   for(Int_t k = 0; k <nybins; k++){
5479     Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
5480     if (fraction>0.95) break;
5481     Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
5482       e*fraction*fraction*fraction*fraction;
5483     sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
5484     sum  += weight*projch->GetBinContent(k+1);
5485     sumCurrent += (Int_t) projch->GetBinContent(k+1);
5486     //printf("fraction %f, weight %f, bincontent %f\n",fraction,weight,projch->GetBinContent(k+1));
5487   }
5488   if(sum > 0.0) fCurrentCoef[0] = (sumw/sum);
5489
5490   if (fDebugLevel == 1) {
5491     TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
5492     cpmeanw->cd();
5493     projch->Draw();
5494   }
5495   fNumberFitSuccess++;
5496 }
5497 //_____________________________________________________________________________
5498 void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
5499 {
5500   //
5501   // Fit methode for the gain factor
5502   //
5503
5504   fCurrentCoef[0]  = 0.0;
5505   fCurrentCoefE    = 0.0;
5506   Double_t chisqrl = 0.0;
5507   Double_t chisqrg = 0.0;
5508   Double_t chisqr  = 0.0;
5509   TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
5510
5511   projch->Fit("landau","0",""
5512              ,(Double_t) mean/fBeginFitCharge
5513              ,projch->GetBinCenter(projch->GetNbinsX()));
5514   Double_t l3P0         = projch->GetFunction("landau")->GetParameter(0);
5515   Double_t l3P1         = projch->GetFunction("landau")->GetParameter(1);
5516   Double_t l3P2         = projch->GetFunction("landau")->GetParameter(2);
5517   chisqrl = projch->GetFunction("landau")->GetChisquare();
5518
5519   projch->Fit("gaus","0",""
5520               ,(Double_t) mean/fBeginFitCharge
5521               ,projch->GetBinCenter(projch->GetNbinsX()));
5522   Double_t g3P0         = projch->GetFunction("gaus")->GetParameter(0);
5523   Double_t g3P2         = projch->GetFunction("gaus")->GetParameter(2);
5524   chisqrg = projch->GetFunction("gaus")->GetChisquare();
5525
5526   fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
5527   if (fDebugLevel != 1) {
5528     projch->Fit("fLandauGaus","0",""
5529                 ,(Double_t) mean/fBeginFitCharge
5530                 ,projch->GetBinCenter(projch->GetNbinsX()));
5531     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5532   }
5533   else  {
5534     TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
5535     cp->cd();
5536     projch->Fit("fLandauGaus","+",""
5537                 ,(Double_t) mean/fBeginFitCharge
5538                 ,projch->GetBinCenter(projch->GetNbinsX()));
5539     chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
5540     projch->Draw();
5541     fLandauGaus->Draw("same");
5542   }
5543
5544   if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5545     //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
5546     fNumberFitSuccess++;
5547     CalculChargeCoefMean(kTRUE);
5548     fCurrentCoef[0]  = projch->GetFunction("fLandauGaus")->GetParameter(1);
5549     fCurrentCoefE    = projch->GetFunction("fLandauGaus")->GetParError(1);
5550   }
5551   else {
5552     CalculChargeCoefMean(kFALSE);
5553     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5554   }
5555
5556   if (fDebugLevel != 1) {
5557     delete fLandauGaus;
5558   }
5559
5560 }
5561 //_____________________________________________________________________________
5562 void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
5563 {
5564   //
5565   // Fit methode for the gain factor more time consuming
5566   //
5567
5568
5569   //Some parameters to initialise
5570   Double_t widthLandau, widthGaus, mPV, integral;
5571   Double_t chisquarel = 0.0;
5572   Double_t chisquareg = 0.0;
5573   projch->Fit("landau","0M+",""
5574               ,(Double_t) mean/6
5575               ,projch->GetBinCenter(projch->GetNbinsX()));
5576   widthLandau  = projch->GetFunction("landau")->GetParameter(2);
5577   chisquarel = projch->GetFunction("landau")->GetChisquare();
5578   projch->Fit("gaus","0M+",""
5579               ,(Double_t) mean/6
5580               ,projch->GetBinCenter(projch->GetNbinsX()));
5581   widthGaus    = projch->GetFunction("gaus")->GetParameter(2);
5582   chisquareg = projch->GetFunction("gaus")->GetChisquare();
5583
5584   mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
5585   integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
5586
5587   // Setting fit range and start values
5588   Double_t fr[2];
5589   //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
5590   //Double_t sv[4]   = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
5591   Double_t sv[4]   = { widthLandau, mPV, integral, widthGaus};
5592   Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
5593   Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
5594   Double_t fp[4]   = { 1.0, 1.0, 1.0, 1.0 };
5595   Double_t fpe[4]  = { 1.0, 1.0, 1.0, 1.0 };
5596   fr[0]            = 0.3 * mean;
5597   fr[1]            = 3.0 * mean;
5598   fCurrentCoef[0]  = 0.0;
5599   fCurrentCoefE    = 0.0;
5600
5601   Double_t chisqr;
5602   Int_t    ndf;
5603   TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
5604                                 ,&pllo[0],&plhi[0]
5605                                 ,&fp[0],&fpe[0]
5606                                 ,&chisqr,&ndf);
5607
5608   Double_t projchPeak;
5609   Double_t projchFWHM;
5610   LanGauPro(fp,projchPeak,projchFWHM);
5611
5612   if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
5613     //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
5614     fNumberFitSuccess++;
5615     CalculChargeCoefMean(kTRUE);
5616     fCurrentCoef[0]  = fp[1];
5617     fCurrentCoefE = fpe[1];
5618     //chargeCoefE2 = chisqr;
5619   }
5620   else {
5621     CalculChargeCoefMean(kFALSE);
5622     fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
5623   }
5624   if (fDebugLevel == 1) {
5625     AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
5626     TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
5627     cpy->cd();
5628     projch->Draw();
5629     fitsnr->Draw("same");
5630   }
5631   else {
5632     delete fitsnr;
5633   }
5634 }
5635 //_____________________________________________________________________________
5636 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
5637 {
5638   //
5639   // Calcul the coefficients of the polynome passant par ces trois points de degre 2
5640   //
5641   Double_t *c = new Double_t[5];
5642   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
5643   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
5644   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
5645
5646   c[4] = 0.0;
5647   c[3] = 0.0;
5648   c[2] = x0+x1+x2;
5649   c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
5650   c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
5651
5652   return c;
5653
5654
5655 }
5656
5657 //_____________________________________________________________________________
5658 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
5659 {
5660   //
5661   // Calcul the coefficients of the polynome passant par ces quatre points de degre 3
5662   //
5663   Double_t *c = new Double_t[5];
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   c[4] = 0.0;
5670   c[3] = x0+x1+x2+x3;
5671   c[2] = -(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   c[1] = (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   c[0] = -(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   return c;
5687
5688
5689 }
5690
5691 //_____________________________________________________________________________
5692 Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
5693 {
5694   //
5695   // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
5696   //
5697   Double_t *c = new Double_t[5];
5698   Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
5699   Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
5700   Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
5701   Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])*(x[3]-x[4]));
5702   Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
5703
5704
5705   c[4] = x0+x1+x2+x3+x4;
5706   c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
5707            +x1*(x[0]+x[2]+x[3]+x[4])
5708            +x2*(x[0]+x[1]+x[3]+x[4])
5709            +x3*(x[0]+x[1]+x[2]+x[4])
5710            +x4*(x[0]+x[1]+x[2]+x[3]));
5711   c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
5712           +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])
5713           +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])
5714           +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])
5715           +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]));
5716
5717   c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
5718           +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])
5719           +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])
5720           +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])
5721           +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]));
5722
5723   c[0] = (x0*x[1]*x[2]*x[3]*x[4]
5724           +x1*x[0]*x[2]*x[3]*x[4]
5725           +x2*x[0]*x[1]*x[3]*x[4]
5726           +x3*x[0]*x[1]*x[2]*x[4]
5727           +x4*x[0]*x[1]*x[2]*x[3]);
5728
5729   return c;
5730
5731
5732 }
5733 //_____________________________________________________________________________
5734 void AliTRDCalibraFit::NormierungCharge()
5735 {
5736   //
5737   // Normalisation of the gain factor resulting for the fits
5738   //
5739
5740   // Calcul of the mean of choosen method by fFitChargeNDB
5741   Double_t sum         = 0.0;
5742   //printf("total number of entries %d\n",fVectorFitCH->GetEntriesFast());
5743   for (Int_t k = 0; k < (Int_t) fVectorFit.GetEntriesFast(); k++) {
5744     Int_t    total    = 0;
5745     Int_t    detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
5746     Float_t *coef     = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
5747     //printf("detector %d coef[0] %f\n",detector,coef[0]);
5748     if (GetStack(detector) == 2) {
5749       total = 1728;
5750     }
5751     if (GetStack(detector) != 2) {
5752       total = 2304;
5753     }
5754     for (Int_t j = 0; j < total; j++) {
5755       if (coef[j] >= 0) {
5756         sum += coef[j];
5757       }
5758     }
5759   }
5760
5761   if (sum > 0) {
5762     fScaleFitFactor = fScaleFitFactor / sum;
5763   }
5764   else {
5765     fScaleFitFactor = 1.0;
5766   }
5767
5768   //methode de boeuf mais bon...
5769   Double_t scalefactor = fScaleFitFactor;
5770
5771   if(fDebugLevel > 1){
5772
5773     if ( !fDebugStreamer ) {
5774       //debug stream
5775       TDirectory *backup = gDirectory;
5776       fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root");
5777       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
5778     }
5779     (* fDebugStreamer) << "NormierungCharge"<<
5780       "scalefactor="<<scalefactor<<
5781       "\n";
5782     }
5783 }
5784 //_____________________________________________________________________________
5785 TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const
5786 {
5787   //
5788   // Rebin of the 1D histo for the gain calibration if needed.
5789   // you have to choose fRebin, divider of fNumberBinCharge
5790   //
5791
5792  TAxis *xhist  = hist->GetXaxis();
5793  TH1I  *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5794                                         ,xhist->GetBinLowEdge(1)
5795                                         ,xhist->GetBinUpEdge(xhist->GetNbins()));
5796
5797  AliInfo(Form("fRebin: %d",fRebin));
5798  Int_t i = 1;
5799  for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5800    Double_t sum = 0.0;
5801    for (Int_t ji = i; ji < i+fRebin; ji++) {
5802      sum += hist->GetBinContent(ji);
5803    }
5804    sum = sum / fRebin;
5805    rehist->SetBinContent(k,sum);
5806    i += fRebin;
5807  }
5808
5809  return rehist;
5810
5811 }
5812
5813 //_____________________________________________________________________________
5814 TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const
5815 {
5816   //
5817   // Rebin of the 1D histo for the gain calibration if needed
5818   // you have to choose fRebin divider of fNumberBinCharge
5819   //
5820
5821   TAxis *xhist  = hist->GetXaxis();
5822   TH1F  *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
5823                                          ,xhist->GetBinLowEdge(1)
5824                                          ,xhist->GetBinUpEdge(xhist->GetNbins()));
5825
5826   AliInfo(Form("fRebin: %d",fRebin));
5827   Int_t i = 1;
5828   for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
5829     Double_t sum = 0.0;
5830     for (Int_t ji = i; ji < i+fRebin; ji++) {
5831       sum += hist->GetBinContent(ji);
5832     }
5833     sum = sum/fRebin;
5834     rehist->SetBinContent(k,sum);
5835     i += fRebin;
5836   }
5837
5838   return rehist;
5839
5840 }
5841 //
5842 //____________Some basic geometry function_____________________________________
5843 //
5844
5845 //_____________________________________________________________________________
5846 Int_t AliTRDCalibraFit::GetLayer(Int_t d) const
5847 {
5848   //
5849   // Reconstruct the plane number from the detector number
5850   //
5851
5852   return ((Int_t) (d % 6));
5853
5854 }
5855
5856 //_____________________________________________________________________________
5857 Int_t AliTRDCalibraFit::GetStack(Int_t d) const
5858 {
5859   //
5860   // Reconstruct the stack number from the detector number
5861   //
5862   const Int_t kNlayer = 6;
5863
5864   return ((Int_t) (d % 30) / kNlayer);
5865
5866 }
5867
5868 //_____________________________________________________________________________
5869 Int_t AliTRDCalibraFit::GetSector(Int_t d) const
5870 {
5871   //
5872   // Reconstruct the sector number from the detector number
5873   //
5874   Int_t fg = 30;
5875
5876   return ((Int_t) (d / fg));
5877
5878 }
5879
5880 //
5881 //____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
5882 //
5883 //_______________________________________________________________________________
5884 void AliTRDCalibraFit::ResetVectorFit()
5885 {
5886   //
5887   // Reset the VectorFits
5888   //
5889
5890   fVectorFit.SetOwner();
5891   fVectorFit.Clear();
5892   fVectorFit2.SetOwner();
5893   fVectorFit2.Clear();
5894
5895 }
5896 //
5897 //____________Private Functions________________________________________________
5898 //
5899
5900 //_____________________________________________________________________________
5901 Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par)
5902 {
5903   //
5904   // Function for the fit
5905   //
5906
5907   //TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
5908
5909   //PARAMETERS FOR FIT PH
5910   // PASAv.4
5911   //fAsymmGauss->SetParameter(0,0.113755);
5912   //fAsymmGauss->SetParameter(1,0.350706);
5913   //fAsymmGauss->SetParameter(2,0.0604244);
5914   //fAsymmGauss->SetParameter(3,7.65596);
5915   //fAsymmGauss->SetParameter(4,1.00124);
5916   //fAsymmGauss->SetParameter(5,0.870597);  // No tail cancelation
5917
5918   Double_t xx = x[0];
5919
5920   if (xx < par[1]) {
5921     return par[5];
5922   }
5923
5924   Double_t dx       = 0.005;
5925   Double_t xs       = par[1];
5926   Double_t ss       = 0.0;
5927   Double_t paras[2] = { 0.0, 0.0 };
5928
5929   while (xs < xx) {
5930     if ((xs >= par[1]) &&
5931         (xs < (par[1]+par[2]))) {
5932       //fAsymmGauss->SetParameter(0,par[0]);
5933       //fAsymmGauss->SetParameter(1,xs);
5934       //ss += fAsymmGauss->Eval(xx);
5935       paras[0] = par[0];
5936       paras[1] = xs;
5937       ss += AsymmGauss(&xx,paras);
5938     }
5939     if ((xs >= (par[1]+par[2])) &&
5940         (xs <  (par[1]+par[2]+par[3]))) {
5941       //fAsymmGauss->SetParameter(0,par[0]*par[4]);
5942       //fAsymmGauss->SetParameter(1,xs);
5943       //ss += fAsymmGauss->Eval(xx);
5944       paras[0] = par[0]*par[4];
5945       paras[1] = xs;
5946       ss += AsymmGauss(&xx,paras);
5947     }
5948     xs += dx;
5949   }
5950
5951   return ss + par[5];
5952
5953 }
5954
5955 //_____________________________________________________________________________
5956 Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par)
5957 {
5958   //
5959   // Function for the fit
5960   //
5961
5962   //par[0] = normalization
5963   //par[1] = mean
5964   //par[2] = sigma
5965   //norm0  = 1
5966   //par[3] = lambda0
5967   //par[4] = norm1
5968   //par[5] = lambda1
5969
5970   Double_t par1save = par[1];
5971   //Double_t par2save = par[2];
5972   Double_t par2save = 0.0604244;
5973   //Double_t par3save = par[3];
5974   Double_t par3save = 7.65596;
5975   //Double_t par5save = par[5];
5976   Double_t par5save = 0.870597;
5977   Double_t dx       = x[0] - par1save;
5978
5979   Double_t  sigma2  = par2save*par2save;
5980   Double_t  sqrt2   = TMath::Sqrt(2.0);
5981   Double_t  exp1    = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
5982                                * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
5983   Double_t  exp2    = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
5984                                * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
5985
5986   //return par[0]*(exp1+par[4]*exp2);
5987   return par[0] * (exp1 + 1.00124 * exp2);
5988
5989 }
5990
5991 //_____________________________________________________________________________
5992 Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par)
5993 {
5994   //
5995   // Sum Landau + Gaus with identical mean
5996   //
5997
5998   Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
5999   //Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[4],par[5]);
6000   Double_t valGaus   = par[3] * TMath::Gaus(x[0],par[1],par[4]);
6001   Double_t val       = valLandau + valGaus;
6002
6003   return val;
6004
6005 }
6006
6007 //_____________________________________________________________________________
6008 Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par)
6009 {
6010   //
6011   // Function for the fit
6012   //
6013   // Fit parameters:
6014   // par[0]=Width (scale) parameter of Landau density
6015   // par[1]=Most Probable (MP, location) parameter of Landau density
6016   // par[2]=Total area (integral -inf to inf, normalization constant)
6017   // par[3]=Width (sigma) of convoluted Gaussian function
6018   //
6019   // In the Landau distribution (represented by the CERNLIB approximation),
6020   // the maximum is located at x=-0.22278298 with the location parameter=0.
6021   // This shift is corrected within this function, so that the actual
6022   // maximum is identical to the MP parameter.
6023   //
6024
6025   // Numeric constants
6026   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
6027   Double_t mpshift  = -0.22278298;       // Landau maximum location
6028
6029   // Control constants
6030   Double_t np       = 100.0;             // Number of convolution steps
6031   Double_t sc       =   5.0;             // Convolution extends to +-sc Gaussian sigmas
6032
6033   // Variables
6034   Double_t xx;
6035   Double_t mpc;
6036   Double_t fland;
6037   Double_t sum = 0.0;
6038   Double_t xlow;
6039   Double_t xupp;
6040   Double_t step;
6041   Double_t i;
6042
6043   // MP shift correction
6044   mpc = par[1] - mpshift * par[0];
6045
6046   // Range of convolution integral
6047   xlow = x[0] - sc * par[3];
6048   xupp = x[0] + sc * par[3];
6049
6050   step = (xupp - xlow) / np;
6051
6052   // Convolution integral of Landau and Gaussian by sum
6053   for (i = 1.0; i <= np/2; i++) {
6054
6055     xx    = xlow + (i-.5) * step;
6056     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6057     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
6058
6059     xx    = xupp - (i-.5) * step;
6060     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
6061     sum  += fland * TMath::Gaus(x[0],xx,par[3]);
6062
6063   }
6064
6065   return (par[2] * step * sum * invsq2pi / par[3]);
6066
6067 }
6068 //_____________________________________________________________________________
6069 TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
6070                                       , const Double_t *parlimitslo, const Double_t *parlimitshi
6071                                       , Double_t *fitparams, Double_t *fiterrors
6072                                       , Double_t *chiSqr, Int_t *ndf) const
6073 {
6074   //
6075   // Function for the fit
6076   //
6077
6078   Int_t i;
6079   Char_t funname[100];
6080
6081   TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
6082   if (ffitold) {
6083     delete ffitold;
6084   }
6085
6086   TF1 *ffit    = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
6087   ffit->SetParameters(startvalues);
6088   ffit->SetParNames("Width","MP","Area","GSigma");
6089
6090   for (i = 0; i < 4; i++) {
6091     ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
6092   }
6093
6094   his->Fit(funname,"RB0");                   // Fit within specified range, use ParLimits, do not plot
6095
6096   ffit->GetParameters(fitparams);            // Obtain fit parameters
6097   for (i = 0; i < 4; i++) {
6098     fiterrors[i] = ffit->GetParError(i);     // Obtain fit parameter errors
6099   }
6100   chiSqr[0] = ffit->GetChisquare();          // Obtain chi^2
6101   ndf[0]    = ffit->GetNDF();                // Obtain ndf
6102
6103   return (ffit);                             // Return fit function
6104
6105 }
6106
6107 //_____________________________________________________________________________
6108 Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
6109 {
6110   //
6111   // Function for the fit
6112   //
6113
6114   Double_t p;
6115   Double_t x;
6116   Double_t fy;
6117   Double_t fxr;
6118   Double_t fxl;
6119   Double_t step;
6120   Double_t l;
6121   Double_t lold;
6122
6123   Int_t    i        = 0;
6124   Int_t    maxcalls = 10000;
6125
6126   // Search for maximum
6127   p    = params[1] - 0.1 * params[0];
6128   step = 0.05 * params[0];
6129   lold = -2.0;
6130   l    = -1.0;
6131
6132   while ((l != lold) && (i < maxcalls)) {
6133     i++;
6134     lold = l;
6135     x    = p + step;
6136     l    = LanGauFun(&x,params);
6137     if (l < lold) {
6138       step = -step / 10.0;
6139     }
6140     p += step;
6141   }
6142
6143   if (i == maxcalls) {
6144     return (-1);
6145   }
6146   maxx = x;
6147   fy = l / 2.0;
6148
6149   // Search for right x location of fy
6150   p    = maxx + params[0];
6151   step = params[0];
6152   lold = -2.0;
6153   l    = -1e300;
6154   i    = 0;
6155
6156   while ( (l != lold) && (i < maxcalls) ) {
6157     i++;
6158
6159     lold = l;
6160     x = p + step;
6161     l = TMath::Abs(LanGauFun(&x,params) - fy);
6162
6163     if (l > lold)
6164       step = -step/10;
6165
6166     p += step;
6167   }
6168
6169   if (i == maxcalls)
6170     return (-2);
6171
6172   fxr = x;
6173
6174
6175   // Search for left x location of fy
6176
6177   p = maxx - 0.5 * params[0];
6178   step = -params[0];
6179   lold = -2.0;
6180   l    = -1.0e300;
6181   i    = 0;
6182
6183   while ((l != lold) && (i < maxcalls)) {
6184     i++;
6185     lold = l;
6186     x    = p + step;
6187     l    = TMath::Abs(LanGauFun(&x,params) - fy);
6188     if (l > lold) {
6189       step = -step / 10.0;
6190     }
6191     p += step;
6192   }
6193
6194   if (i == maxcalls) {
6195     return (-3);
6196   }
6197
6198   fxl  = x;
6199   fwhm = fxr - fxl;
6200
6201   return (0);
6202 }
6203 //_____________________________________________________________________________
6204 Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
6205 {
6206   //
6207   // Gaus with identical mean
6208   //
6209
6210   Double_t gauss   = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
6211
6212   return gauss;
6213
6214 }
6215