TRD/AliTRDCalibraFillHisto.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 // AliTRDCalibraFillHisto
  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 fills histos or vectors.
  24 // It can be used for the calibration per chamber but also per group of pads and eventually per pad.
  25 // The user has to choose with the functions SetNz and SetNrphi the precision of the calibration (see AliTRDCalibraMode).
  26 // 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
  27 // from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
  28 // in the function "FollowBackProlongation" (AliTRDtracker)
  29 // Per default the functions to fill are off.
  30 //
  31 // Author:
  32 //   R. Bailhache (R.Bailhache@gsi.de)
  33 //
  34 //////////////////////////////////////////////////////////////////////////////////////
  35
  36 #include <TTree.h>
  37 #include <TProfile2D.h>
  38 #include <TProfile.h>
  39 #include <TFile.h>
  40 #include <TChain.h>
  41 #include <TStyle.h>
  42 #include <TCanvas.h>
  43 #include <TGraphErrors.h>
  44 #include <TObjArray.h>
  45 #include <TH1F.h>
  46 #include <TH2I.h>
  47 #include <TH2.h>
  48 #include <TStopwatch.h>
  49 #include <TMath.h>
  50 #include <TDirectory.h>
  51 #include <TROOT.h>
  52 #include <TTreeStream.h>
  53 #include <TVectorD.h>
  54 #include <TF1.h>
  55
  56 #include "AliLog.h"
  57 #include "AliCDBManager.h"
  58
  59 #include "AliTRDCalibraFillHisto.h"
  60 #include "AliTRDCalibraFit.h"
  61 #include "AliTRDCalibraMode.h"
  62 #include "AliTRDCalibraVector.h"
  63 #include "AliTRDCalibraVdriftLinearFit.h"
  64 #include "AliTRDcalibDB.h"
  65 #include "AliTRDCommonParam.h"
  66 #include "AliTRDmcmTracklet.h"
  67 #include "AliTRDpadPlane.h"
  68 #include "AliTRDcluster.h"
  69 #include "AliTRDtrack.h"
  70 #include "AliTRDRawStream.h"
  71 #include "AliRawReader.h"
  72 #include "AliRawReaderDate.h"
  73 #include "AliTRDgeometry.h"
  74
  75 #ifdef ALI_DATE
  76 #include "event.h"
  77 #endif
  78
  79
  80 ClassImp(AliTRDCalibraFillHisto)
  81
  82 AliTRDCalibraFillHisto* AliTRDCalibraFillHisto::fgInstance = 0;
  83 Bool_t AliTRDCalibraFillHisto::fgTerminated = kFALSE;
  84
  85 //_____________singleton implementation_________________________________________________
  86 AliTRDCalibraFillHisto *AliTRDCalibraFillHisto::Instance()
  87 {
  88   //
  89   // Singleton implementation
  90   //
  91
  92   if (fgTerminated != kFALSE) {
  93     return 0;
  94   }
  95
  96   if (fgInstance == 0) {
  97     fgInstance = new AliTRDCalibraFillHisto();
  98   }
  99
 100   return fgInstance;
 101
 102 }
 103
 104 //______________________________________________________________________________________
 105 void AliTRDCalibraFillHisto::Terminate()
 106 {
 107   //
 108   // Singleton implementation
 109   // Deletes the instance of this class
 110   //
 111
 112   fgTerminated = kTRUE;
 113
 114   if (fgInstance != 0) {
 115     delete fgInstance;
 116     fgInstance = 0;
 117   }
 118
 119 }
 120
 121 //______________________________________________________________________________________
 122 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
 123   :TObject()
 124   ,fGeo(0)
 125   ,fMITracking(kFALSE)
 126   ,fMcmTracking(kFALSE)
 127   ,fMcmCorrectAngle(kFALSE)
 128   ,fCH2dOn(kFALSE)
 129   ,fPH2dOn(kFALSE)
 130   ,fPRF2dOn(kFALSE)
 131   ,fHisto2d(kFALSE)
 132   ,fVector2d(kFALSE)
 133   ,fLinearFitterOn(kFALSE)
 134   ,fLinearFitterDebugOn(kFALSE)
 135   ,fRelativeScale(0)
 136   ,fThresholdClusterPRF2(15.0)
 137   ,fCalibraMode(new AliTRDCalibraMode())
 138   ,fDebugStreamer(0)
 139   ,fDebugLevel(0)
 140   ,fDetectorAliTRDtrack(kFALSE)
 141   ,fDetectorPreviousTrack(-1)
 142   ,fNumberClusters(18)
 143   ,fProcent(6.0)
 144   ,fDifference(17)
 145   ,fNumberTrack(0)
 146   ,fTimeMax(0)
 147   ,fSf(10.0)
 148   ,fNumberBinCharge(100)
 149   ,fNumberBinPRF(40)
 150   ,fNgroupprf(0)
 151   ,fListClusters(new TObjArray())
 152   ,fPar0(0x0)
 153   ,fPar1(0x0)
 154   ,fPar2(0x0)
 155   ,fPar3(0x0)
 156   ,fPar4(0x0)
 157   ,fAmpTotal(0x0)
 158   ,fPHPlace(0x0)
 159   ,fPHValue(0x0)
 160   ,fGoodTracklet(kTRUE)
 161   ,fGoodTrack(kTRUE)
 162   ,fEntriesCH(0x0)
 163   ,fEntriesLinearFitter(0x0)
 164   ,fCalibraVector(0x0)
 165   ,fPH2d(0x0)
 166   ,fPRF2d(0x0)
 167   ,fCH2d(0x0)
 168   ,fLinearFitterArray(0)
 169   ,fLinearVdriftFit(0x0)
 170 {
 171   //
 172   // Default constructor
 173   //
 174
 175   //
 176   // Init some default values
 177   //
 178
 179   fNumberUsedCh[0]       = 0;
 180   fNumberUsedCh[1]       = 0;
 181   fNumberUsedPh[0]       = 0;
 182   fNumberUsedPh[1]       = 0;
 183
 184   fGeo = new AliTRDgeometry();
 185
 186 }
 187
 188 //______________________________________________________________________________________
 189 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
 190   :TObject(c)
 191   ,fGeo(0)
 192   ,fMITracking(c.fMITracking)
 193   ,fMcmTracking(c.fMcmTracking)
 194   ,fMcmCorrectAngle(c.fMcmCorrectAngle)
 195   ,fCH2dOn(c.fCH2dOn)
 196   ,fPH2dOn(c.fPH2dOn)
 197   ,fPRF2dOn(c.fPRF2dOn)
 198   ,fHisto2d(c.fHisto2d)
 199   ,fVector2d(c.fVector2d)
 200   ,fLinearFitterOn(c.fLinearFitterOn)
 201   ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
 202   ,fRelativeScale(c.fRelativeScale)
 203   ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
 204   ,fCalibraMode(0x0)
 205   ,fDebugStreamer(0)
 206   ,fDebugLevel(c.fDebugLevel)
 207   ,fDetectorAliTRDtrack(c.fDetectorAliTRDtrack)
 208   ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
 209   ,fNumberClusters(c.fNumberClusters)
 210   ,fProcent(c.fProcent)
 211   ,fDifference(c.fDifference)
 212   ,fNumberTrack(c.fNumberTrack)
 213   ,fTimeMax(c.fTimeMax)
 214   ,fSf(c.fSf)
 215   ,fNumberBinCharge(c.fNumberBinCharge)
 216   ,fNumberBinPRF(c.fNumberBinPRF)
 217   ,fNgroupprf(c.fNgroupprf)
 218   ,fListClusters(new TObjArray())
 219   ,fPar0(0x0)
 220   ,fPar1(0x0)
 221   ,fPar2(0x0)
 222   ,fPar3(0x0)
 223   ,fPar4(0x0)
 224   ,fAmpTotal(0x0)
 225   ,fPHPlace(0x0)
 226   ,fPHValue(0x0)
 227   ,fGoodTracklet(c.fGoodTracklet)
 228   ,fGoodTrack(c.fGoodTrack)
 229   ,fEntriesCH(0x0)
 230   ,fEntriesLinearFitter(0x0)
 231   ,fCalibraVector(0x0)
 232   ,fPH2d(0x0)
 233   ,fPRF2d(0x0)
 234   ,fCH2d(0x0)
 235   ,fLinearFitterArray(0)
 236   ,fLinearVdriftFit(0x0)
 237 {
 238   //
 239   // Copy constructor
 240   //
 241   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
 242   if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
 243   if(c.fPH2d) {
 244     fPH2d = new TProfile2D(*c.fPH2d);
 245     fPH2d->SetDirectory(0);
 246   }
 247   if(c.fPRF2d) {
 248     fPRF2d = new TProfile2D(*c.fPRF2d);
 249     fPRF2d->SetDirectory(0);
 250   }
 251   if(c.fCH2d) {
 252     fCH2d = new TH2I(*c.fCH2d);
 253     fCH2d->SetDirectory(0);
 254   }
 255   if(c.fLinearVdriftFit){
 256     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
 257   }
 258
 259   if (fGeo) {
 260     delete fGeo;
 261   }
 262   fGeo = new AliTRDgeometry();
 263 }
 264
 265 //____________________________________________________________________________________
 266 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
 267 {
 268   //
 269   // AliTRDCalibraFillHisto destructor
 270   //
 271
 272   ClearHistos();
 273   if ( fDebugStreamer ) delete fDebugStreamer;
 274
 275   if (fGeo) {
 276     delete fGeo;
 277   }
 278
 279 }
 280
 281 //_____________________________________________________________________________
 282 void AliTRDCalibraFillHisto::Destroy()
 283 {
 284   //
 285   // Delete instance
 286   //
 287
 288   if (fgInstance) {
 289     delete fgInstance;
 290     fgInstance = 0x0;
 291   }
 292
 293 }
 294
 295 //_____________________________________________________________________________
 296 void AliTRDCalibraFillHisto::ClearHistos()
 297 {
 298   //
 299   // Delete the histos
 300   //
 301
 302   if (fPH2d) {
 303     delete fPH2d;
 304     fPH2d  = 0x0;
 305   }
 306   if (fCH2d) {
 307     delete fCH2d;
 308     fCH2d  = 0x0;
 309   }
 310   if (fPRF2d) {
 311     delete fPRF2d;
 312     fPRF2d = 0x0;
 313   }
 314
 315 }
 316
 317 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
 318 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
 319 {
 320   //
 321   // For the offline tracking
 322   // This function will be called in the function AliReconstruction::Run()
 323   // Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
 324   //
 325
 326   // DB Setting
 327   // Get cal
 328   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 329   if (!cal) {
 330     AliInfo("Could not get calibDB");
 331     return kFALSE;
 332   }
 333   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
 334   if (!parCom) {
 335     AliInfo("Could not get CommonParam");
 336     return kFALSE;
 337   }
 338
 339   // Some parameters
 340   fTimeMax = cal->GetNumberOfTimeBins();
 341   fSf      = parCom->GetSamplingFrequency();
 342   fRelativeScale = 20;
 343
 344   //Init the tracklet parameters
 345   fPar0 = new Double_t[fTimeMax];
 346   fPar1 = new Double_t[fTimeMax];
 347   fPar2 = new Double_t[fTimeMax];
 348   fPar3 = new Double_t[fTimeMax];
 349   fPar4 = new Double_t[fTimeMax];
 350
 351   for(Int_t k = 0; k < fTimeMax; k++){
 352     fPar0[k] = 0.0;
 353     fPar1[k] = 0.0;
 354     fPar2[k] = 0.0;
 355     fPar3[k] = 0.0;
 356     fPar4[k] = 0.0;
 357   }
 358
 359   // Calcul Xbins Chambd0, Chamb2
 360   Int_t Ntotal0 = CalculateTotalNumberOfBins(0);
 361   Int_t Ntotal1 = CalculateTotalNumberOfBins(1);
 362   Int_t Ntotal2 = CalculateTotalNumberOfBins(2);
 363
 364   // If vector method On initialised all the stuff
 365   if(fVector2d){
 366     fCalibraVector = new AliTRDCalibraVector();
 367     fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
 368     fCalibraVector->SetTimeMax(fTimeMax);
 369     if(fNgroupprf != 0) {
 370       fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
 371       fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
 372     }
 373     else {
 374       fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
 375       fCalibraVector->SetPRFRange(1.5);
 376     }
 377     for(Int_t k = 0; k < 3; k++){
 378       fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
 379       fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
 380     }
 381     TString namech("Nz");
 382     namech += fCalibraMode->GetNz(0);
 383     namech += "Nrphi";
 384     namech += fCalibraMode->GetNrphi(0);
 385     fCalibraVector->SetNameCH((const char* ) namech);
 386     TString nameph("Nz");
 387     nameph += fCalibraMode->GetNz(1);
 388     nameph += "Nrphi";
 389     nameph += fCalibraMode->GetNrphi(1);
 390     fCalibraVector->SetNamePH((const char* ) nameph);
 391     TString nameprf("Nz");
 392     nameprf += fCalibraMode->GetNz(2);
 393     nameprf += "Nrphi";
 394     nameprf += fCalibraMode->GetNrphi(2);
 395     nameprf += "Ngp";
 396     nameprf += fNgroupprf;
 397     fCalibraVector->SetNamePRF((const char* ) nameprf);
 398   }
 399
 400   // Create the 2D histos corresponding to the pad groupCalibration mode
 401   if (fCH2dOn) {
 402
 403     AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
 404                 ,fCalibraMode->GetNz(0)
 405                 ,fCalibraMode->GetNrphi(0)));
 406
 407     // Create the 2D histo
 408     if (fHisto2d) {
 409       CreateCH2d(Ntotal0);
 410     }
 411     // Variable
 412     fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
 413     for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
 414       fAmpTotal[k] = 0.0;
 415     }
 416     //Statistics
 417     fEntriesCH = new Int_t[Ntotal0];
 418     for(Int_t k = 0; k < Ntotal0; k++){
 419       fEntriesCH[k] = 0;
 420     }
 421
 422   }
 423   if (fPH2dOn) {
 424
 425     AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
 426                 ,fCalibraMode->GetNz(1)
 427                 ,fCalibraMode->GetNrphi(1)));
 428
 429     // Create the 2D histo
 430     if (fHisto2d) {
 431       CreatePH2d(Ntotal1);
 432     }
 433     // Variable
 434     fPHPlace = new Short_t[fTimeMax];
 435     for (Int_t k = 0; k < fTimeMax; k++) {
 436       fPHPlace[k] = -1;
 437     }
 438     fPHValue = new Float_t[fTimeMax];
 439     for (Int_t k = 0; k < fTimeMax; k++) {
 440       fPHValue[k] = 0.0;
 441     }
 442   }
 443   if (fLinearFitterOn) {
 444     fLinearFitterArray.Expand(540);
 445     fLinearFitterArray.SetName("ArrayLinearFitters");
 446     fEntriesLinearFitter = new Int_t[540];
 447     for(Int_t k = 0; k < 540; k++){
 448       fEntriesLinearFitter[k] = 0;
 449     }
 450     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
 451   }
 452
 453   if (fPRF2dOn) {
 454
 455     AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
 456                 ,fCalibraMode->GetNz(2)
 457                 ,fCalibraMode->GetNrphi(2)));
 458     // Create the 2D histo
 459     if (fHisto2d) {
 460       CreatePRF2d(Ntotal2);
 461     }
 462   }
 463
 464   return kTRUE;
 465
 466 }
 467
 468 //____________Functions for filling the histos in the code_____________________
 469
 470 //____________Offine tracking in the AliTRDtracker_____________________________
 471 Bool_t AliTRDCalibraFillHisto::ResetTrack()
 472 {
 473   //
 474   // For the offline tracking
 475   // This function will be called in the function
 476   // AliTRDtracker::FollowBackPropagation() at the beginning
 477   // Reset the parameter to know we have a new TRD track
 478   //
 479
 480   fDetectorAliTRDtrack = kFALSE;
 481   fGoodTrack           = kTRUE;
 482   return kTRUE;
 483
 484 }
 485 //____________Offine tracking in the AliTRDtracker_____________________________
 486 void AliTRDCalibraFillHisto::ResetfVariables()
 487 {
 488   //
 489   // Reset values per tracklet
 490   //
 491
 492   // Reset the list of clusters
 493   fListClusters->Clear();
 494
 495   //Reset the tracklet parameters
 496   for(Int_t k = 0; k < fTimeMax; k++){
 497     fPar0[k] = 0.0;
 498     fPar1[k] = 0.0;
 499     fPar2[k] = 0.0;
 500     fPar3[k] = 0.0;
 501     fPar4[k] = 0.0;
 502   }
 503
 504
 505   //Reset good tracklet
 506   fGoodTracklet = kTRUE;
 507
 508   // Reset the fPHValue
 509   if (fPH2dOn) {
 510     //Reset the fPHValue and fPHPlace
 511     for (Int_t k = 0; k < fTimeMax; k++) {
 512       fPHValue[k] = 0.0;
 513       fPHPlace[k] = -1;
 514     }
 515   }
 516
 517   // Reset the fAmpTotal where we put value
 518   if (fCH2dOn) {
 519     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
 520       fAmpTotal[k] = 0.0;
 521     }
 522   }
 523
 524 }
 525 //____________Offline tracking in the AliTRDtracker____________________________
 526 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
 527 {
 528   //
 529   // For the offline tracking
 530   // This function will be called in the function
 531   // AliTRDtracker::FollowBackPropagation() in the loop over the clusters
 532   // of TRD tracks
 533   // Fill the 2D histos or the vectors with the info of the clusters at
 534   // the end of a detectors if the track is "good"
 535   //
 536
 537   // Localisation of the detector
 538   Int_t detector = cl->GetDetector();
 539
 540
 541   // Fill the infos for the previous clusters if not the same
 542   // detector anymore or if not the same track
 543   if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
 544       (fDetectorPreviousTrack != -1)) {
 545
 546     fNumberTrack++;
 547
 548     // If the same track, then look if the previous detector is in
 549     // the same plane, if yes: not a good track
 550     //FollowBack
 551     if (fDetectorAliTRDtrack &&
 552         (GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
 553     //Follow
 554     //if (fDetectorAliTRDtrack &&
 555     //    (GetPlane(detector) >= GetPlane(fDetectorPreviousTrack))) {
 556       fGoodTrack = kFALSE;
 557     }
 558
 559     // Fill only if the track doesn't touch a masked pad or doesn't
 560     // appear in the middle (fGoodTrack)
 561     if (fGoodTrack && fGoodTracklet) {
 562
 563       // drift velocity unables to cut bad tracklets
 564       Bool_t  pass = FindP1TrackPH();
 565
 566       // Gain calibration
 567       if (fCH2dOn) {
 568         FillTheInfoOfTheTrackCH();
 569       }
 570
 571       // PH calibration
 572       if (fPH2dOn) {
 573         FillTheInfoOfTheTrackPH();
 574       }
 575
 576       if(pass && fPRF2dOn) HandlePRF();
 577
 578
 579     } // if a good track
 580
 581     ResetfVariables();
 582
 583   } // Fill at the end the charge
 584
 585   // Calcul the position of the detector
 586   if (detector != fDetectorPreviousTrack) {
 587     LocalisationDetectorXbins(detector);
 588   }
 589
 590   // Reset the detector
 591   fDetectorPreviousTrack = detector;
 592   fDetectorAliTRDtrack   = kTRUE;
 593
 594   // Store the infos of the tracklets
 595   AliTRDcluster *kcl = new AliTRDcluster(*cl);
 596   fListClusters->Add((TObject *)kcl);
 597   Int_t time = cl->GetLocalTimeBin();
 598   fPar0[time] = t->GetY();
 599   fPar1[time] = t->GetZ();
 600   fPar2[time] = t->GetSnp();
 601   fPar3[time] = t->GetTgl();
 602   fPar4[time] = t->Get1Pt();
 603
 604   // Store the info bis of the tracklet
 605   Int_t *rowcol   = CalculateRowCol(cl);
 606   CheckGoodTracklet(detector,rowcol);
 607   Int_t     group[2] = {0,0};
 608   if(fCH2dOn)  group[0]  = CalculateCalibrationGroup(0,rowcol);
 609   if(fPH2dOn)  group[1]  = CalculateCalibrationGroup(1,rowcol);
 610   StoreInfoCHPH(cl,t,group);
 611
 612   return kTRUE;
 613
 614 }
 615 //____________Online trackling in AliTRDtrigger________________________________
 616 Bool_t AliTRDCalibraFillHisto::UpdateHistogramcm(AliTRDmcmTracklet *trk)
 617 {
 618   //
 619   // For the tracking
 620   // This function will be called in the function AliTRDtrigger::TestTracklet
 621   // before applying the pt cut on the tracklets
 622   // Fill the infos for the tracklets fTrkTest if the tracklets is "good"
 623   //
 624
 625   // Localisation of the Xbins involved
 626   Int_t idect = trk->GetDetector();
 627   fDetectorPreviousTrack = idect;
 628   LocalisationDetectorXbins(idect);
 629
 630   // Get the parameter object
 631   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 632   if (!cal) {
 633     AliInfo("Could not get calibDB");
 634     return kFALSE;
 635   }
 636
 637   // Reset
 638   ResetfVariables();
 639
 640   // Row of the tracklet and position in the pad groups
 641   Int_t *rowcol  = new Int_t[2];
 642   rowcol[0]     = trk->GetRow();
 643   Int_t group[3] = {-1,-1,-1};
 644
 645   // Eventuelle correction due to track angle in z direction
 646   Float_t correction = 1.0;
 647   if (fMcmCorrectAngle) {
 648     Float_t z = trk->GetRowz();
 649     Float_t r = trk->GetTime0();
 650     correction = r / TMath::Sqrt((r*r+z*z));
 651   }
 652
 653   // Boucle sur les clusters
 654   // Condition on number of cluster: don't come from the middle of the detector
 655   if (trk->GetNclusters() >= fNumberClusters) {
 656
 657     for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
 658
 659       Float_t amp[3] = { 0.0, 0.0, 0.0 };
 660       Int_t   time   = trk->GetClusterTime(icl);
 661       rowcol[1]      = trk->GetClusterCol(icl);
 662
 663       amp[0] = trk->GetClusterADC(icl)[0] * correction;
 664       amp[1] = trk->GetClusterADC(icl)[1] * correction;
 665       amp[2] = trk->GetClusterADC(icl)[2] * correction;
 666
 667
 668       if ((amp[0] < 0.0) ||
 669           (amp[1] < 0.0) ||
 670           (amp[2] < 0.0)) {
 671         continue;
 672       }
 673
 674       // Col of cluster and position in the pad groups
 675       if(fCH2dOn)  {
 676         group[0] = CalculateCalibrationGroup(0,rowcol);
 677         fAmpTotal[(Int_t) group[0]] += (Float_t) (amp[0]+amp[1]+amp[2]);
 678       }
 679       if(fPH2dOn)  {
 680         group[1] = CalculateCalibrationGroup(1,rowcol);
 681         fPHPlace[time] = group[1];
 682         fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
 683       }
 684       if(fPRF2dOn) group[2] = CalculateCalibrationGroup(2,rowcol);
 685
 686       // See if we are not near a masked pad fGoodTracklet
 687       CheckGoodTracklet(idect,rowcol);
 688
 689       // Fill PRF direct without tnp bins...only for monitoring...
 690       if (fPRF2dOn && fGoodTracklet) {
 691
 692         if ((amp[0] > fThresholdClusterPRF2) &&
 693             (amp[1] > fThresholdClusterPRF2) &&
 694             (amp[2] > fThresholdClusterPRF2) &&
 695             ((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
 696
 697           // Security of the denomiateur is 0
 698           if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
 699              / ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
 700             Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
 701                                   / (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
 702             Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
 703
 704             if (TMath::Abs(xcenter) < 0.5) {
 705               Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
 706               Float_t ymax   = amp[2] / (amp[0]+amp[1]+amp[2]);
 707               // Fill only if it is in the drift region!
 708               if (((Float_t) time / fSf) > 0.3) {
 709                 if (fHisto2d) {
 710                   fPRF2d->Fill(xcenter,(fCalibraMode->GetXbins(2)+group[2]+0.5),ycenter);
 711                   fPRF2d->Fill(-(xcenter+1.0),(fCalibraMode->GetXbins(2)+group[2]+0.5),yminus);
 712                   fPRF2d->Fill((1.0-xcenter),(fCalibraMode->GetXbins(2)+group[2]+0.5),ymax);
 713                 }
 714                 if (fVector2d) {
 715                   fCalibraVector->UpdateVectorPRF(idect,group[2],xcenter,ycenter);
 716                   fCalibraVector->UpdateVectorPRF(idect,group[2],-(xcenter+1.0),yminus);
 717                   fCalibraVector->UpdateVectorPRF(idect,group[2],(1.0-xcenter),ymax);
 718                 }
 719               }//in the drift region
 720             }//in the middle
 721           }//denominateur security
 722         }//cluster shape and thresholds
 723       }//good and PRF On
 724
 725     } // Boucle clusters
 726
 727     // Fill the charge
 728     if(fGoodTracklet){
 729       if (fCH2dOn) FillTheInfoOfTheTrackCH();
 730       if (fPH2dOn) FillTheInfoOfTheTrackPH();
 731     }
 732
 733     fNumberTrack++;
 734
 735   } // Condition on number of clusters
 736
 737   return kTRUE;
 738
 739 }
 740 //_____________________________________________________________________________
 741 Int_t *AliTRDCalibraFillHisto::CalculateRowCol(AliTRDcluster *cl) const
 742 {
 743   //
 744   // Calculate the row and col number of the cluster
 745   //
 746
 747
 748   Int_t *rowcol = new Int_t[2];
 749   rowcol[0] =  0;
 750   rowcol[1] =  0;
 751
 752   // Localisation of the detector
 753   Int_t detector = cl->GetDetector();
 754   Int_t chamber  = GetChamber(detector);
 755   Int_t plane    = GetPlane(detector);
 756
 757   // Localisation of the cluster
 758   Double_t pos[3] = { 0.0, 0.0, 0.0 };
 759   pos[0] = ((AliCluster *)cl)->GetX();
 760   pos[1] = cl->GetY();
 761   pos[2] = cl->GetZ();
 762
 763   // Position of the cluster
 764   AliTRDpadPlane *padplane  = fGeo->GetPadPlane(plane,chamber);
 765   Int_t    row              = padplane->GetPadRowNumber(pos[2]);
 766   //Do not take from here because it was corrected from ExB already....
 767   //Double_t offsetz         = padplane->GetPadRowOffset(row,pos[2]);
 768   //Double_t offsettilt      = padplane->GetTiltOffset(offsetz);
 769   //Int_t    col             = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
 770   //Int_t    col             = padplane->GetPadColNumber(pos[1]+offsettilt);
 771   Int_t    col               = cl->GetPad();
 772
 773   //return
 774   rowcol[0]     = row;
 775   rowcol[1]     = col;
 776   return rowcol;
 777
 778 }
 779 //_____________________________________________________________________________
 780 void AliTRDCalibraFillHisto::CheckGoodTracklet(Int_t detector, Int_t *rowcol)
 781 {
 782   //
 783   // See if we are not near a masked pad
 784   //
 785
 786   Int_t row = rowcol[0];
 787   Int_t col = rowcol[1];
 788
 789   if (!IsPadOn(detector, col, row)) {
 790     fGoodTracklet = kFALSE;
 791   }
 792
 793   if (col > 0) {
 794     if (!IsPadOn(detector, col-1, row)) {
 795       fGoodTracklet = kFALSE;
 796     }
 797   }
 798
 799   if (col < 143) {
 800     if (!IsPadOn(detector, col+1, row)) {
 801       fGoodTracklet = kFALSE;
 802     }
 803   }
 804
 805 }
 806 //_____________________________________________________________________________
 807 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t col, Int_t row) const
 808 {
 809   //
 810   // Look in the choosen database if the pad is On.
 811   // If no the track will be "not good"
 812   //
 813
 814   // Get the parameter object
 815   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 816   if (!cal) {
 817     AliInfo("Could not get calibDB");
 818     return kFALSE;
 819   }
 820
 821   if (!cal->IsChamberInstalled(detector)     ||
 822        cal->IsChamberMasked(detector)        ||
 823        cal->IsPadMasked(detector,col,row)) {
 824     return kFALSE;
 825   }
 826   else {
 827     return kTRUE;
 828   }
 829
 830 }
 831 //_____________________________________________________________________________
 832 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t *rowcol) const
 833 {
 834   //
 835   // Calculate the calibration group number for i
 836   //
 837
 838   // Row of the cluster and position in the pad groups
 839   Int_t posr = 0;
 840   if (fCalibraMode->GetNnZ(i) != 0) {
 841     posr = (Int_t) rowcol[0] / fCalibraMode->GetNnZ(i);
 842   }
 843
 844
 845   // Col of the cluster and position in the pad groups
 846   Int_t posc = 0;
 847   if (fCalibraMode->GetNnRphi(i) != 0) {
 848     posc = (Int_t) rowcol[1] / fCalibraMode->GetNnRphi(i);
 849   }
 850
 851   return posc*fCalibraMode->GetNfragZ(i)+posr;
 852
 853 }
 854 //____________________________________________________________________________________
 855 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
 856 {
 857   //
 858   // Calculate the total number of calibration groups
 859   //
 860
 861   Int_t Ntotal = 0;
 862   fCalibraMode->ModePadCalibration(2,i);
 863   fCalibraMode->ModePadFragmentation(0,2,0,i);
 864   fCalibraMode->SetDetChamb2(i);
 865   Ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
 866   fCalibraMode->ModePadCalibration(0,i);
 867   fCalibraMode->ModePadFragmentation(0,0,0,i);
 868   fCalibraMode->SetDetChamb0(i);
 869   Ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
 870   AliInfo(Form("Total number of Xbins: %d for i %d",Ntotal,i));
 871   return Ntotal;
 872
 873 }
 874 //____________Set the pad calibration variables for the detector_______________
 875 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
 876 {
 877   //
 878   // For the detector calcul the first Xbins and set the number of row
 879   // and col pads per calibration groups, the number of calibration
 880   // groups in the detector.
 881   //
 882
 883   // first Xbins of the detector
 884   if (fCH2dOn) {
 885     fCalibraMode->CalculXBins(detector,0);
 886   }
 887   if (fPH2dOn) {
 888     fCalibraMode->CalculXBins(detector,1);
 889   }
 890   if (fPRF2dOn) {
 891     fCalibraMode->CalculXBins(detector,2);
 892   }
 893
 894   // fragmentation of idect
 895   for (Int_t i = 0; i < 3; i++) {
 896     fCalibraMode->ModePadCalibration((Int_t) GetChamber(detector),i);
 897     fCalibraMode->ModePadFragmentation((Int_t) GetPlane(detector)
 898                        , (Int_t) GetChamber(detector)
 899                        , (Int_t) GetSector(detector),i);
 900   }
 901
 902   return kTRUE;
 903
 904 }
 905 //_____________________________________________________________________________
 906 void AliTRDCalibraFillHisto::StoreInfoCHPH(AliTRDcluster *cl, AliTRDtrack *t, Int_t *group)
 907 {
 908   //
 909   // Store the infos in fAmpTotal, fPHPlace and fPHValue
 910   //
 911
 912   // Charge in the cluster
 913   Float_t  q        = TMath::Abs(cl->GetQ());
 914   Int_t    time     = cl->GetLocalTimeBin();
 915
 916   // Correction due to the track angle
 917   Float_t correction    = 1.0;
 918   Float_t normalisation = 6.67;
 919   if ((q >0) && (t->GetNdedx() > 0)) {
 920     correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (normalisation);
 921   }
 922
 923   // Fill the fAmpTotal with the charge
 924   if (fCH2dOn) {
 925     fAmpTotal[(Int_t) group[0]] += correction;
 926   }
 927
 928   // Fill the fPHPlace and value
 929   if (fPH2dOn) {
 930     fPHPlace[time] = group[1];
 931     fPHValue[time] = correction;
 932   }
 933
 934 }
 935 //_____________________________________________________________________
 936 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDRawStream *rawStream, Bool_t nocheck)
 937 {
 938   //
 939   // Event Processing loop - AliTRDRawStream
 940   // 0 timebin problem
 941   // 1 no input
 942   // 2 input
 943   //
 944
 945   Int_t withInput = 1;
 946
 947   Int_t phvalue[36];
 948   //Int_t row[36];
 949   //Int_t col[36];
 950   for(Int_t k = 0; k < 36; k++){
 951     phvalue[k] = 10;
 952     //row[k]     = -1;
 953     //col[36]    = -1;
 954   }
 955   fDetectorPreviousTrack = -1;
 956   Int_t nbtimebin = 0;
 957
 958   if(!nocheck){
 959
 960     fTimeMax = 0;
 961
 962     while (rawStream->Next()) {
 963
 964       Int_t idetector = rawStream->GetDet();                            //  current detector
 965       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
 966         if(TMath::Mean(fTimeMax,phvalue)>0.0){
 967           withInput = 2;
 968           for(Int_t k = 0; k < fTimeMax; k++){
 969             UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax);
 970             phvalue[k] = 10;
 971             //row[k]     = -1;
 972             //col[k]     = -1;
 973           }
 974         }
 975       }
 976       fDetectorPreviousTrack = idetector;
 977       nbtimebin         = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
 978       if(nbtimebin == 0) return 0;
 979       if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
 980       fTimeMax          = nbtimebin;
 981       Int_t iTimeBin    = rawStream->GetTimeBin();                       //  current time bin
 982       //row[iTimeBin]   = rawStream->GetRow();                           //  current row
 983       //col[iTimeBin]   = rawStream->GetCol();                           //  current col
 984       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
 985
 986       Int_t fin     = TMath::Min(fTimeMax,(iTimeBin+3));
 987       Int_t n       = 0;
 988       for(Int_t itime = iTimeBin; itime < fin; itime++){
 989         // should extract baseline here!
 990         if(signal[n]>5.0) phvalue[itime] = signal[n];
 991         n++;
 992       }
 993     }
 994
 995     // fill the last one
 996     if(fDetectorPreviousTrack != -1){
 997       if(TMath::Mean(fTimeMax,phvalue)>0.0){
 998         withInput = 2;
 999         for(Int_t k = 0; k < fTimeMax; k++){
1000           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax);
1001           phvalue[k] = 10;
1002           //row[k]     = -1;
1003           //col[k]     = -1;
1004         }
1005       }
1006     }
1007
1008   }
1009   else{
1010
1011     while (rawStream->Next()) {
1012
1013       Int_t idetector = rawStream->GetDet();                            //  current detector
1014       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
1015         if(TMath::Mean(nbtimebin,phvalue)>0.0){
1016           withInput = 2;
1017           for(Int_t k = 0; k < nbtimebin; k++){
1018             UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1019             phvalue[k] = 10;
1020             //row[k]     = -1;
1021             //col[k]     = -1;
1022           }
1023         }
1024       }
1025       fDetectorPreviousTrack = idetector;
1026       nbtimebin         = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
1027       Int_t iTimeBin    = rawStream->GetTimeBin();                       //  current time bin
1028       //row[iTimeBin]   = rawStream->GetRow();                           //  current row
1029       //col[iTimeBin]   = rawStream->GetCol();                           //  current col
1030       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
1031
1032       Int_t fin     = TMath::Min(nbtimebin,(iTimeBin+3));
1033       Int_t n       = 0;
1034       for(Int_t itime = iTimeBin; itime < fin; itime++){
1035         // should extract baseline here!
1036         if(signal[n]>5.0) phvalue[itime] = signal[n];
1037         n++;
1038       }
1039     }
1040
1041     // fill the last one
1042     if(fDetectorPreviousTrack != -1){
1043       if(TMath::Mean(nbtimebin,phvalue)>0.0){
1044         withInput = 2;
1045         for(Int_t k = 0; k < nbtimebin; k++){
1046           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1047           phvalue[k] = 10;
1048           //row[k]     = -1;
1049           //col[k]     = -1;
1050         }
1051       }
1052     }
1053   }
1054
1055   return withInput;
1056
1057 }
1058 //_____________________________________________________________________
1059 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
1060 {
1061   //
1062   //  Event processing loop - AliRawReader
1063   //
1064
1065
1066   AliTRDRawStream rawStream(rawReader);
1067
1068   rawReader->Select("TRD");
1069
1070   return ProcessEventDAQ(&rawStream, nocheck);
1071 }
1072 //_________________________________________________________________________
1073 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
1074 #ifdef ALI_DATE
1075                                                eventHeaderStruct *event,
1076                                                Bool_t nocheck
1077 #else
1078                                                eventHeaderStruct* /*event*/,
1079                                                Bool_t /*nocheck*/
1080
1081 #endif
1082                                    )
1083 {
1084   //
1085   //  process date event
1086   //
1087 #ifdef ALI_DATE
1088     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1089     Int_t result=ProcessEventDAQ(rawReader, nocheck);
1090     delete rawReader;
1091     return result;
1092 #else
1093     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1094     return 0;
1095 #endif
1096
1097 }
1098 //____________Online trackling in AliTRDtrigger________________________________
1099 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1100 {
1101   //
1102   // For the DAQ
1103   // Fill a simple average pulse height
1104   //
1105
1106   // Localisation of the Xbins involved
1107   //LocalisationDetectorXbins(det);
1108
1109   // Row  and position in the pad groups
1110   //Int_t posr = 0;
1111   //if (fCalibraMode->GetNnZ(1) != 0) {
1112   //  posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1113   //}
1114
1115   // Col of cluster and position in the pad groups
1116   //Int_t posc = 0;
1117   //if (fCalibraMode->GetNnRphi(1) != 0) {
1118   //  posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1119   //}
1120
1121   //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1122
1123   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1124
1125   return kTRUE;
1126
1127 }
1128 //____________Write_____________________________________________________
1129 //_____________________________________________________________________
1130 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1131 {
1132   //
1133   //  Write infos to file
1134   //
1135
1136   //For debugging
1137   if ( fDebugStreamer ) {
1138     delete fDebugStreamer;
1139     fDebugStreamer = 0x0;
1140   }
1141
1142   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1143                ,fNumberTrack
1144                ,fNumberUsedCh[0]
1145                ,fNumberUsedCh[1]
1146                ,fNumberUsedPh[0]
1147                ,fNumberUsedPh[1]));
1148
1149   TDirectory *backup = gDirectory;
1150   TString option;
1151
1152   if ( append )
1153     option = "update";
1154   else
1155     option = "recreate";
1156
1157   TFile f(filename,option.Data());
1158
1159   TStopwatch stopwatch;
1160   stopwatch.Start();
1161   if(fVector2d) {
1162     f.WriteTObject(fCalibraVector);
1163   }
1164
1165   if (fCH2dOn ) {
1166     if (fHisto2d) {
1167       f.WriteTObject(fCH2d);
1168     }
1169   }
1170   if (fPH2dOn ) {
1171     if (fHisto2d) {
1172       f.WriteTObject(fPH2d);
1173     }
1174   }
1175   if (fPRF2dOn) {
1176     if (fHisto2d) {
1177         f.WriteTObject(fPRF2d);
1178     }
1179   }
1180   if(fLinearFitterOn){
1181     AnalyseLinearFitter();
1182     f.WriteTObject(fLinearVdriftFit);
1183   }
1184
1185   f.Close();
1186
1187   if ( backup ) backup->cd();
1188
1189   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1190                ,stopwatch.RealTime(),stopwatch.CpuTime()));
1191 }
1192 //___________________________________________probe the histos__________________________________________________
1193 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1194 {
1195   //
1196   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1197   // debug mode with 2 for TH2I and 3 for TProfile2D
1198   // It gives a pointer to a Double_t[7] with the info following...
1199   // [0] : number of calibration groups with entries
1200   // [1] : minimal number of entries found
1201   // [2] : calibration group number of the min
1202   // [3] : maximal number of entries found
1203   // [4] : calibration group number of the max
1204   // [5] : mean number of entries found
1205   // [6] : mean relative error
1206   //
1207
1208   Double_t *info = new Double_t[7];
1209
1210   // Number of Xbins (detectors or groups of pads)
1211   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
1212   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
1213
1214   // Initialise
1215   Double_t nbwe = 0; //number of calibration groups with entries
1216   Double_t minentries = 0; //minimal number of entries found
1217   Double_t maxentries = 0; //maximal number of entries found
1218   Double_t placemin = 0; //calibration group number of the min
1219   Double_t placemax = -1; //calibration group number of the max
1220   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1221   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1222
1223   Double_t counter = 0;
1224
1225   //Debug
1226   TH1F *NbEntries = 0x0;//distribution of the number of entries
1227   TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1228   TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1229
1230   // Beginning of the loop over the calibration groups
1231   for (Int_t idect = 0; idect < nbins; idect++) {
1232
1233     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1234     projch->SetDirectory(0);
1235
1236     // Number of entries for this calibration group
1237     Double_t nentries = 0.0;
1238     if((i%2) == 0){
1239       for (Int_t k = 0; k < nxbins; k++) {
1240         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1241       }
1242     }
1243     else{
1244       for (Int_t k = 0; k < nxbins; k++) {
1245         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1246         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1247           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1248           counter++;
1249         }
1250       }
1251     }
1252
1253     //Debug
1254     if(i > 1){
1255       if((!((Bool_t)NbEntries)) && (nentries > 0)){
1256         NbEntries = new TH1F("Number of entries","Number of entries"
1257                                ,100,(Int_t)nentries/2,nentries*2);
1258         NbEntries->SetDirectory(0);
1259         NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1260                                ,nbins,0,nbins);
1261         NbEntriesPerGroup->SetDirectory(0);
1262         NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1263                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1264         NbEntriesPerSp->SetDirectory(0);
1265       }
1266       if(NbEntries){
1267         if(nentries > 0) NbEntries->Fill(nentries);
1268         NbEntriesPerGroup->Fill(idect+0.5,nentries);
1269         NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1270       }
1271     }
1272
1273     //min amd max
1274     if(nentries > maxentries){
1275       maxentries = nentries;
1276       placemax = idect;
1277     }
1278     if(idect == 0) {
1279       minentries = nentries;
1280     }
1281     if(nentries < minentries){
1282       minentries = nentries;
1283       placemin = idect;
1284     }
1285     //nbwe
1286     if(nentries > 0) {
1287       nbwe++;
1288       meanstats += nentries;
1289     }
1290   }//calibration groups loop
1291
1292   if(nbwe > 0) meanstats /= nbwe;
1293   if(counter > 0) meanrelativerror /= counter;
1294
1295   AliInfo(Form("There are %f calibration groups with entries",nbwe));
1296   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1297   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1298   AliInfo(Form("The mean number of entries is %f",meanstats));
1299   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1300
1301   info[0] = nbwe;
1302   info[1] = minentries;
1303   info[2] = placemin;
1304   info[3] = maxentries;
1305   info[4] = placemax;
1306   info[5] = meanstats;
1307   info[6] = meanrelativerror;
1308
1309   if(i > 1){
1310     gStyle->SetPalette(1);
1311     gStyle->SetOptStat(1111);
1312     gStyle->SetPadBorderMode(0);
1313     gStyle->SetCanvasColor(10);
1314     gStyle->SetPadLeftMargin(0.13);
1315     gStyle->SetPadRightMargin(0.01);
1316     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1317     stat->Divide(2,1);
1318     stat->cd(1);
1319     NbEntries->Draw("");
1320     stat->cd(2);
1321     NbEntriesPerSp->SetStats(0);
1322     NbEntriesPerSp->Draw("");
1323     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1324     stat1->cd();
1325     NbEntriesPerGroup->SetStats(0);
1326     NbEntriesPerGroup->Draw("");
1327   }
1328
1329   return info;
1330
1331 }
1332 //____________________________________________________________________________
1333 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1334 {
1335   //
1336   // Return a Int_t[4] with:
1337   // 0 Mean number of entries
1338   // 1 median of number of entries
1339   // 2 rms of number of entries
1340   // 3 number of group with entries
1341   //
1342
1343   Double_t *stat      = new Double_t[4];
1344   stat[3]             = 0.0;
1345
1346   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
1347   Double_t *weight    = new Double_t[nbofgroups];
1348   Int_t    *nonul     = new Int_t[nbofgroups];
1349
1350   for(Int_t k = 0; k < nbofgroups; k++){
1351     if(fEntriesCH[k] > 0) {
1352       weight[k] = 1.0;
1353       nonul[(Int_t)stat[3]] = fEntriesCH[k];
1354       stat[3]++;
1355     }
1356     else weight[k] = 0.0;
1357   }
1358   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
1359   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
1360   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
1361
1362   return stat;
1363
1364 }
1365 //____________________________________________________________________________
1366 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1367 {
1368   //
1369   // Return a Int_t[4] with:
1370   // 0 Mean number of entries
1371   // 1 median of number of entries
1372   // 2 rms of number of entries
1373   // 3 number of group with entries
1374   //
1375
1376   Double_t *stat      = new Double_t[4];
1377   stat[3]             = 0.0;
1378
1379   Int_t    nbofgroups = 540;
1380   Double_t *weight    = new Double_t[nbofgroups];
1381   Int_t    *nonul     = new Int_t[nbofgroups];
1382
1383   for(Int_t k = 0; k < nbofgroups; k++){
1384     if(fEntriesLinearFitter[k] > 0) {
1385       weight[k] = 1.0;
1386       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1387       stat[3]++;
1388     }
1389     else weight[k] = 0.0;
1390   }
1391   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1392   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1393   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
1394
1395   return stat;
1396
1397 }
1398 //_____________________________________________________________________________
1399 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1400 {
1401   //
1402   // Should be between 0 and 6
1403   //
1404
1405   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1406     AliInfo("The number of groups must be between 0 and 6!");
1407   }
1408   else {
1409     fNgroupprf = numberGroupsPRF;
1410   }
1411
1412 }
1413 //_____________________________________________________________________________
1414 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1415 {
1416   //
1417   // Set the factor that will divide the deposited charge
1418   // to fit in the histo range [0,300]
1419   //
1420
1421   if (RelativeScale > 0.0) {
1422     fRelativeScale = RelativeScale;
1423   }
1424   else {
1425     AliInfo("RelativeScale must be strict positif!");
1426   }
1427
1428 }
1429
1430 //_____________________________________________________________________________
1431 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1432 {
1433   //
1434   // Set the mode of calibration group in the z direction for the parameter i
1435   //
1436
1437   if ((Nz >= 0) &&
1438       (Nz <  5)) {
1439     fCalibraMode->SetNz(i, Nz);
1440   }
1441   else {
1442     AliInfo("You have to choose between 0 and 4");
1443   }
1444
1445 }
1446
1447 //_____________________________________________________________________________
1448 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1449 {
1450   //
1451   // Set the mode of calibration group in the rphi direction for the parameter i
1452   //
1453
1454   if ((Nrphi >= 0) &&
1455       (Nrphi <  7)) {
1456     fCalibraMode->SetNrphi(i ,Nrphi);
1457   }
1458   else {
1459     AliInfo("You have to choose between 0 and 6");
1460   }
1461
1462 }
1463 //____________Protected Functions______________________________________________
1464 //____________Create the 2D histo to be filled online__________________________
1465 //
1466 //_____________________________________________________________________________
1467 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1468 {
1469   //
1470   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1471   // If fNgroupprf is zero then no binning in tnp
1472   //
1473
1474   TString name("Nz");
1475   name += fCalibraMode->GetNz(2);
1476   name += "Nrphi";
1477   name += fCalibraMode->GetNrphi(2);
1478   name += "Ngp";
1479   name += fNgroupprf;
1480
1481   if(fNgroupprf != 0){
1482
1483     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1484                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1485     fPRF2d->SetYTitle("Det/pad groups");
1486     fPRF2d->SetXTitle("Position x/W [pad width units]");
1487     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1488     fPRF2d->SetStats(0);
1489   }
1490   else{
1491     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1492                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1493     fPRF2d->SetYTitle("Det/pad groups");
1494     fPRF2d->SetXTitle("Position x/W [pad width units]");
1495     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1496     fPRF2d->SetStats(0);
1497   }
1498
1499 }
1500
1501 //_____________________________________________________________________________
1502 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1503 {
1504   //
1505   // Create the 2D histos
1506   //
1507
1508   TString name("Nz");
1509   name += fCalibraMode->GetNz(1);
1510   name += "Nrphi";
1511   name += fCalibraMode->GetNrphi(1);
1512
1513   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1514                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1515                          ,nn,0,nn);
1516   fPH2d->SetYTitle("Det/pad groups");
1517   fPH2d->SetXTitle("time [#mus]");
1518   fPH2d->SetZTitle("<PH> [a.u.]");
1519   fPH2d->SetStats(0);
1520
1521 }
1522 //_____________________________________________________________________________
1523 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1524 {
1525   //
1526   // Create the 2D histos
1527   //
1528
1529   TString name("Nz");
1530   name += fCalibraMode->GetNz(0);
1531   name += "Nrphi";
1532   name += fCalibraMode->GetNrphi(0);
1533
1534   fCH2d = new TH2I("CH2d",(const Char_t *) name
1535                    ,fNumberBinCharge,0,300,nn,0,nn);
1536   fCH2d->SetYTitle("Det/pad groups");
1537   fCH2d->SetXTitle("charge deposit [a.u]");
1538   fCH2d->SetZTitle("counts");
1539   fCH2d->SetStats(0);
1540   fCH2d->Sumw2();
1541
1542 }
1543
1544 //____________Offine tracking in the AliTRDtracker_____________________________
1545 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1546 {
1547   //
1548   // For the offline tracking or mcm tracklets
1549   // This function will be called in the functions UpdateHistogram...
1550   // to fill the info of a track for the relativ gain calibration
1551   //
1552
1553   Int_t nb            =  0;   // Nombre de zones traversees
1554   Int_t fd            = -1;   // Premiere zone non nulle
1555   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1556
1557
1558
1559
1560   // See if the track goes through different zones
1561   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1562     if (fAmpTotal[k] > 0.0) {
1563       totalcharge += fAmpTotal[k];
1564       nb++;
1565       if (nb == 1) {
1566         fd = k;
1567       }
1568     }
1569   }
1570
1571
1572   switch (nb)
1573     {
1574     case 1:
1575       fNumberUsedCh[0]++;
1576       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1577       if (fHisto2d) {
1578         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1579         //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1580       }
1581       if (fVector2d) {
1582         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1583       }
1584       break;
1585     case 2:
1586       if ((fAmpTotal[fd]   > 0.0) &&
1587           (fAmpTotal[fd+1] > 0.0)) {
1588         // One of the two very big
1589         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1590           if (fHisto2d) {
1591             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1592             //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1593           }
1594           if (fVector2d) {
1595             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1596           }
1597           fNumberUsedCh[1]++;
1598           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1599         }
1600         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1601           if (fHisto2d) {
1602             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1603             //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1604           }
1605           if (fVector2d) {
1606             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1607           }
1608           fNumberUsedCh[1]++;
1609           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1610         }
1611       }
1612       if (fCalibraMode->GetNfragZ(0) > 1) {
1613         if (fAmpTotal[fd] > 0.0) {
1614           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1615             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1616               // One of the two very big
1617               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1618                 if (fHisto2d) {
1619                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1620                   //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1621                 }
1622                 if (fVector2d) {
1623                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1624                 }
1625                 fNumberUsedCh[1]++;
1626                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1627               }
1628               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1629                 if (fHisto2d) {
1630                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1631                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1632                 }
1633                 fNumberUsedCh[1]++;
1634                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1635                 if (fVector2d) {
1636                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1637                 }
1638               }
1639             }
1640           }
1641         }
1642       }
1643       break;
1644     default: break;
1645     }
1646 }
1647 //____________Offine tracking in the AliTRDtracker_____________________________
1648 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1649 {
1650   //
1651   // For the offline tracking or mcm tracklets
1652   // This function will be called in the functions UpdateHistogram...
1653   // to fill the info of a track for the drift velocity  calibration
1654   //
1655
1656   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1657   Int_t fd1 = -1; // Premiere zone non nulle
1658   Int_t fd2 = -1; // Deuxieme zone non nulle
1659   Int_t k1  = -1; // Debut de la premiere zone
1660   Int_t k2  = -1; // Debut de la seconde zone
1661
1662   // See if the track goes through different zones
1663   for (Int_t k = 0; k < fTimeMax; k++) {
1664     if (fPHValue[k] > 0.0) {
1665       if (fd1 == -1) {
1666         fd1 = fPHPlace[k];
1667         k1  = k;
1668       }
1669       if (fPHPlace[k] != fd1) {
1670         if (fd2 == -1) {
1671           k2  = k;
1672           fd2 = fPHPlace[k];
1673           nb  = 2;
1674         }
1675         if (fPHPlace[k] != fd2) {
1676           nb = 3;
1677         }
1678       }
1679     }
1680   }
1681
1682
1683   switch(nb)
1684     {
1685     case 1:
1686       fNumberUsedPh[0]++;
1687       for (Int_t i = 0; i < fTimeMax; i++) {
1688         if (fHisto2d) {
1689           fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1690         }
1691         if (fVector2d) {
1692           fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1693         }
1694       }
1695       break;
1696     case 2:
1697       if ((fd1 == fd2+1) ||
1698           (fd2 == fd1+1)) {
1699         // One of the two fast all the think
1700         if (k2 > (k1+fDifference)) {
1701           //we choose to fill the fd1 with all the values
1702           fNumberUsedPh[1]++;
1703           for (Int_t i = 0; i < fTimeMax; i++) {
1704             if (fHisto2d) {
1705               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1706             }
1707             if (fVector2d) {
1708               fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1709             }
1710           }
1711         }
1712         if ((k2+fDifference) < fTimeMax) {
1713           //we choose to fill the fd2 with all the values
1714           fNumberUsedPh[1]++;
1715           for (Int_t i = 0; i < fTimeMax; i++) {
1716             if (fHisto2d) {
1717               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1718             }
1719           if (fVector2d) {
1720             fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1721           }
1722           }
1723         }
1724       }
1725       // Two zones voisines sinon rien!
1726       if (fCalibraMode->GetNfragZ(1) > 1) {
1727         // Case 2
1728         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
1729           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
1730             // One of the two fast all the think
1731             if (k2 > (k1+fDifference)) {
1732               //we choose to fill the fd1 with all the values
1733               fNumberUsedPh[1]++;
1734               for (Int_t i = 0; i < fTimeMax; i++) {
1735                 if (fHisto2d) {
1736                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1737                 }
1738                 if (fVector2d) {
1739                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1740                 }
1741               }
1742             }
1743             if ((k2+fDifference) < fTimeMax) {
1744               //we choose to fill the fd2 with all the values
1745               fNumberUsedPh[1]++;
1746               for (Int_t i = 0; i < fTimeMax; i++) {
1747                 if (fHisto2d) {
1748                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1749                 }
1750                 if (fVector2d) {
1751                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1752                 }
1753               }
1754             }
1755           }
1756         }
1757         // Two zones voisines sinon rien!
1758         // Case 3
1759         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
1760           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
1761             // One of the two fast all the think
1762             if (k2 > (k1 + fDifference)) {
1763               //we choose to fill the fd1 with all the values
1764               fNumberUsedPh[1]++;
1765               for (Int_t i = 0; i < fTimeMax; i++) {
1766                 if (fHisto2d) {
1767                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1768                 }
1769                 if (fVector2d) {
1770                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1771                 }
1772               }
1773             }
1774             if ((k2+fDifference) < fTimeMax) {
1775               //we choose to fill the fd2 with all the values
1776               fNumberUsedPh[1]++;
1777               for (Int_t i = 0; i < fTimeMax; i++) {
1778                 if (fHisto2d) {
1779                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1780                 }
1781                 if (fVector2d) {
1782                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1783                 }
1784               }
1785             }
1786           }
1787         }
1788       }
1789       break;
1790     default: break;
1791     }
1792 }
1793 //____________Offine tracking in the AliTRDtracker_____________________________
1794 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
1795 {
1796   //
1797   // For the offline tracking
1798   // This function will be called in the functions UpdateHistogram...
1799   // to fill the find the parameter P1 of a track for the drift velocity  calibration
1800   //
1801
1802
1803   //Number of points: if less than 3 return kFALSE
1804   Int_t Npoints = fListClusters->GetEntriesFast();
1805   if(Npoints <= 2) return kFALSE;
1806
1807   //Variables
1808   TLinearFitter linearFitterTracklet      = TLinearFitter(2,"pol1");        // TLinearFitter per tracklet
1809   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
1810   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
1811   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
1812   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
1813   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
1814   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!
1815   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
1816   Double_t pointError                 = 0.0;                                // error after straight line fit
1817   Int_t    detector                   = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
1818   Int_t    snpright                   = 1;                                  // if we took in the middle snp
1819   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
1820   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
1821   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
1822   Int_t     nbli                      = 0;                                  // number linear fitter points
1823   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
1824
1825   linearFitterTracklet.StoreData(kFALSE);
1826   linearFitterTracklet.ClearPoints();
1827
1828   //if more than one row
1829   Int_t    rowp                       = -1;                              // if it crosses a pad row
1830
1831   //tiltingangle
1832   tiltingangle                        = padplane->GetTiltingAngle();
1833   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
1834
1835   //Fill with points
1836   for(Int_t k = 0; k < Npoints; k++){
1837
1838     AliTRDcluster *cl                 = (AliTRDcluster *) fListClusters->At(k);
1839     Double_t ycluster                 = cl->GetY();
1840     Int_t time                        = cl->GetLocalTimeBin();
1841     Double_t timeis                   = time/fSf;
1842     //See if cross two pad rows
1843     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
1844     if(k==0) rowp                     = row;
1845     if(row != rowp) crossrow          = 1;
1846     //Take in the middle of the chamber
1847     //FollowBack
1848     if(time > (Int_t) 10) {
1849     //Follow
1850     //if(time < (Int_t) 11) {
1851       z   = cl->GetZ();
1852       y   = cl->GetY();
1853       snp = fPar2[time];
1854       tgl = fPar3[time];
1855     }
1856     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
1857     nbli++;
1858   }
1859   //FollowBack
1860   if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
1861   //Follow
1862   //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
1863   if(nbli <= 2) return kFALSE;
1864
1865   // Do the straight line fit now
1866   TVectorD pars;
1867   linearFitterTracklet.Eval();
1868   linearFitterTracklet.GetParameters(pars);
1869   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
1870   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
1871   dydt  = pars[1];
1872
1873   if( TMath::Abs(snp) <  1.){
1874     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
1875   }
1876   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1877
1878   if(fDebugLevel > 0){
1879     if ( !fDebugStreamer ) {
1880       //debug stream
1881       TDirectory *backup = gDirectory;
1882       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1883       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1884     }
1885
1886     (* fDebugStreamer) << "VDRIFT0"<<
1887       "Npoints="<<Npoints<<
1888       "\n";
1889
1890
1891     (* fDebugStreamer) << "VDRIFT"<<
1892       "snpright="<<snpright<<
1893       "Npoints="<<Npoints<<
1894       "nbli="<<nbli<<
1895       "detector="<<detector<<
1896       "snp="<<snp<<
1897       "tnp="<<tnp<<
1898       "tgl="<<tgl<<
1899       "tnt="<<tnt<<
1900       "y="<<y<<
1901       "z="<<z<<
1902       "dydt="<<dydt<<
1903       "dzdx="<<dzdx<<
1904       "crossrow="<<crossrow<<
1905       "errorpar="<<errorpar<<
1906       "pointError="<<pointError<<
1907       "\n";
1908
1909   }
1910
1911   if(Npoints < fNumberClusters) return kFALSE;
1912   if(snpright == 0) return kFALSE;
1913   if(pointError >= 0.1) return kFALSE;
1914   if(crossrow == 1) return kFALSE;
1915
1916   if(fLinearFitterOn){
1917     //Add to the linear fitter of the detector
1918     if( TMath::Abs(snp) <  1.){
1919       Double_t x = tnp-dzdx*tnt;
1920       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
1921       if(fLinearFitterDebugOn) {
1922         fLinearVdriftFit->Update(detector,x,pars[1]);
1923       }
1924       fEntriesLinearFitter[detector]++;
1925     }
1926   }
1927   //AliInfo("End of FindP1TrackPH with success!")
1928   return kTRUE;
1929
1930 }
1931 //____________Offine tracking in the AliTRDtracker_____________________________
1932 Bool_t AliTRDCalibraFillHisto::HandlePRF()
1933 {
1934   //
1935   // For the offline tracking
1936   // Fit the tracklet with a line and take the position as reference for the PRF
1937   //
1938
1939   //Number of points
1940   Int_t Npoints  = fListClusters->GetEntriesFast();                         // number of total points
1941   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit
1942   Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
1943
1944
1945   // To see the difference due to the fit
1946   Double_t *padPositions;
1947   padPositions = new Double_t[Npoints];
1948   for(Int_t k = 0; k < Npoints; k++){
1949     padPositions[k] = 0.0;
1950   }
1951
1952
1953   //Find the position by a fit
1954   TLinearFitter fitter(2,"pol1");
1955   fitter.StoreData(kFALSE);
1956   fitter.ClearPoints();
1957   for(Int_t k = 0;  k < Npoints; k++){
1958     //Take the cluster
1959     AliTRDcluster *cl  = (AliTRDcluster *) fListClusters->At(k);
1960     Short_t  *signals  = cl->GetSignals();
1961     Double_t     time  = cl->GetLocalTimeBin();
1962     //Calculate x if possible
1963     Float_t xcenter    = 0.0;
1964     Bool_t  echec      = kTRUE;
1965     if((time<=7) || (time>=21)) continue;
1966     // Center 3 balanced: position with the center of the pad
1967     if ((((Float_t) signals[3]) > 0.0) &&
1968         (((Float_t) signals[2]) > 0.0) &&
1969         (((Float_t) signals[4]) > 0.0)) {
1970       echec = kFALSE;
1971       // Security if the denomiateur is 0
1972       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1973            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1974         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1975           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1976                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1977       }
1978       else {
1979         echec = kTRUE;
1980       }
1981     }
1982     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1983     if(echec) continue;
1984     //if no echec: calculate with the position of the pad
1985     // Position of the cluster
1986     Double_t       padPosition = xcenter +  cl->GetPad();
1987     padPositions[k]            = padPosition;
1988     Nb3pc++;
1989     fitter.AddPoint(&time, padPosition,1);
1990   }//clusters loop
1991
1992   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
1993   if(Nb3pc < 3) return kFALSE;
1994   fitter.Eval();
1995   TVectorD line(2);
1996   fitter.GetParameters(line);
1997   Float_t  pointError  = -1.0;
1998   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
1999
2000
2001   // Now fill the PRF
2002   for(Int_t k = 0;  k < Npoints; k++){
2003     //Take the cluster
2004     AliTRDcluster *cl      = (AliTRDcluster *) fListClusters->At(k);
2005     Short_t  *signals      = cl->GetSignals();              // signal
2006     Double_t     time      = cl->GetLocalTimeBin();         // time bin
2007     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
2008     Float_t padPos         = cl->GetPad();                  // middle pad
2009     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
2010     Float_t ycenter        = 0.0;                           // relative center charge
2011     Float_t ymin           = 0.0;                           // relative left charge
2012     Float_t ymax           = 0.0;                           // relative right charge
2013     Double_t tgl           = fPar3[(Int_t)time];            // dz/dl and not dz/dx
2014     Double_t pt            = fPar4[(Int_t)time];            // pt
2015     Float_t  dzdx          = 0.0;                           // dzdx
2016
2017
2018     //Requiere simply two pads clusters at least
2019     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2020        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2021       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2022       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2023       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
2024       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
2025     }
2026
2027     //calibration group
2028     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2029     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2030     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
2031     Double_t  snp          = fPar2[(Int_t)time];                        // sin angle in xy plan
2032     Float_t   xcl          = cl->GetY();                                // y cluster
2033     Float_t   qcl          = cl->GetQ();                                // charge cluster
2034     Int_t     plane        = GetPlane(detector);                        // plane
2035     Int_t     chamber      = GetChamber(detector);                      // chamber
2036     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2037     Double_t  x            = dpad;                                      // reconstructed position moved
2038     Float_t   Ep           = pointError;                                // error of fit
2039     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2040     Float_t   signal3      = (Float_t)signals[3];                       // signal
2041     Float_t   signal2      = (Float_t)signals[2];                       // signal
2042     Float_t   signal4      = (Float_t)signals[4];                       // signal
2043     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2044     Float_t tnp            = 0.0;
2045     if(TMath::Abs(snp) < 1.0){
2046       tnp = snp / (TMath::Sqrt(1-snp*snp));
2047       dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2048     }
2049
2050
2051     if(fDebugLevel > 0){
2052       if ( !fDebugStreamer ) {
2053         //debug stream
2054         TDirectory *backup = gDirectory;
2055         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2056         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2057       }
2058
2059       (* fDebugStreamer) << "PRF0"<<
2060         "caligroup="<<caligroup<<
2061         "detector="<<detector<<
2062         "plane="<<plane<<
2063         "chamber="<<chamber<<
2064         "Npoints="<<Npoints<<
2065         "Np="<<Nb3pc<<
2066         "Ep="<<Ep<<
2067         "snp="<<snp<<
2068         "tnp="<<tnp<<
2069         "tgl="<<tgl<<
2070         "dzdx="<<dzdx<<
2071         "pt="<<pt<<
2072         "padPos="<<padPos<<
2073         "padPosition="<<padPositions[k]<<
2074         "padPosTracklet="<<padPosTracklet<<
2075         "x="<<x<<
2076         "ycenter="<<ycenter<<
2077         "ymin="<<ymin<<
2078         "ymax="<<ymax<<
2079         "xcl="<<xcl<<
2080         "qcl="<<qcl<<
2081         "signal1="<<signal1<<
2082         "signal2="<<signal2<<
2083         "signal3="<<signal3<<
2084         "signal4="<<signal4<<
2085         "signal5="<<signal5<<
2086         "time="<<time<<
2087         "\n";
2088       x = xdiff;
2089       Int_t type=0;
2090       Float_t y = ycenter;
2091       (* fDebugStreamer) << "PRFALL"<<
2092         "caligroup="<<caligroup<<
2093         "detector="<<detector<<
2094         "plane="<<plane<<
2095         "chamber="<<chamber<<
2096         "Npoints="<<Npoints<<
2097         "Np="<<Nb3pc<<
2098         "Ep="<<Ep<<
2099         "type="<<type<<
2100         "snp="<<snp<<
2101         "tnp="<<tnp<<
2102         "tgl="<<tgl<<
2103         "dzdx="<<dzdx<<
2104         "pt="<<pt<<
2105         "padPos="<<padPos<<
2106         "padPosition="<<padPositions[k]<<
2107         "padPosTracklet="<<padPosTracklet<<
2108         "x="<<x<<
2109         "y="<<y<<
2110         "xcl="<<xcl<<
2111         "qcl="<<qcl<<
2112         "signal1="<<signal1<<
2113         "signal2="<<signal2<<
2114         "signal3="<<signal3<<
2115         "signal4="<<signal4<<
2116         "signal5="<<signal5<<
2117         "time="<<time<<
2118         "\n";
2119       x=-(xdiff+1);
2120       y = ymin;
2121       type=-1;
2122       (* fDebugStreamer) << "PRFALL"<<
2123         "caligroup="<<caligroup<<
2124         "detector="<<detector<<
2125         "plane="<<plane<<
2126         "chamber="<<chamber<<
2127         "Npoints="<<Npoints<<
2128         "Np="<<Nb3pc<<
2129         "Ep="<<Ep<<
2130         "type="<<type<<
2131         "snp="<<snp<<
2132         "tnp="<<tnp<<
2133         "tgl="<<tgl<<
2134         "dzdx="<<dzdx<<
2135         "pt="<<pt<<
2136         "padPos="<<padPos<<
2137         "padPosition="<<padPositions[k]<<
2138         "padPosTracklet="<<padPosTracklet<<
2139         "x="<<x<<
2140         "y="<<y<<
2141         "xcl="<<xcl<<
2142         "qcl="<<qcl<<
2143         "signal1="<<signal1<<
2144         "signal2="<<signal2<<
2145         "signal3="<<signal3<<
2146         "signal4="<<signal4<<
2147         "signal5="<<signal5<<
2148         "time="<<time<<
2149         "\n";
2150       x=1-xdiff;
2151       y = ymax;
2152       type=1;
2153
2154       (* fDebugStreamer) << "PRFALL"<<
2155         "caligroup="<<caligroup<<
2156         "detector="<<detector<<
2157         "plane="<<plane<<
2158         "chamber="<<chamber<<
2159         "Npoints="<<Npoints<<
2160         "Np="<<Nb3pc<<
2161         "Ep="<<Ep<<
2162         "type="<<type<<
2163         "snp="<<snp<<
2164         "tnp="<<tnp<<
2165         "tgl="<<tgl<<
2166         "dzdx="<<dzdx<<
2167         "pt="<<pt<<
2168         "padPos="<<padPos<<
2169         "padPosition="<<padPositions[k]<<
2170         "padPosTracklet="<<padPosTracklet<<
2171         "x="<<x<<
2172         "y="<<y<<
2173         "xcl="<<xcl<<
2174         "qcl="<<qcl<<
2175         "signal1="<<signal1<<
2176         "signal2="<<signal2<<
2177         "signal3="<<signal3<<
2178         "signal4="<<signal4<<
2179         "signal5="<<signal5<<
2180         "time="<<time<<
2181         "\n";
2182
2183     }
2184
2185     // some cuts
2186     if(Npoints < fNumberClusters) continue;
2187     if(Nb3pc <= 5) continue;
2188     if((time >= 21) || (time < 7)) continue;
2189     if(TMath::Abs(snp) >= 1.0) continue;
2190     if(qcl < 80) continue;
2191
2192     Bool_t echec   = kFALSE;
2193     Double_t shift = 0.0;
2194     //Calculate the shift in x coresponding to this tnp
2195     if(fNgroupprf != 0.0){
2196       shift      = -3.0*(fNgroupprf-1)-1.5;
2197       Double_t limithigh  = -0.2*(fNgroupprf-1);
2198       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2199       else{
2200         while(tnp > limithigh){
2201           limithigh += 0.2;
2202           shift += 3.0;
2203         }
2204       }
2205     }
2206     if (fHisto2d && !echec) {
2207       if(TMath::Abs(dpad) < 1.5) {
2208         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2209         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2210       }
2211       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2212         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2213         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2214       }
2215       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2216         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2217         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2218       }
2219     }
2220     //Not equivalent anymore here!
2221     if (fVector2d && !echec) {
2222       if(TMath::Abs(dpad) < 1.5) {
2223         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2224         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2225       }
2226       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2227         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2228         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2229       }
2230       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
2231         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2232         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2233       }
2234     }
2235   }
2236   return kTRUE;
2237
2238 }
2239 //
2240 //____________Some basic geometry function_____________________________________
2241 //
2242 //_____________________________________________________________________________
2243 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2244 {
2245   //
2246   // Reconstruct the plane number from the detector number
2247   //
2248
2249   return ((Int_t) (d % 6));
2250
2251 }
2252
2253 //_____________________________________________________________________________
2254 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2255 {
2256   //
2257   // Reconstruct the chamber number from the detector number
2258   //
2259   Int_t fgkNplan = 6;
2260
2261   return ((Int_t) (d % 30) / fgkNplan);
2262
2263 }
2264
2265 //_____________________________________________________________________________
2266 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2267 {
2268   //
2269   // Reconstruct the sector number from the detector number
2270   //
2271   Int_t fg = 30;
2272
2273   return ((Int_t) (d / fg));
2274
2275 }
2276 //_____________________________________________________________________
2277 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2278 {
2279     //
2280     // return pointer to fPH2d TProfile2D
2281     // create a new TProfile2D if it doesn't exist allready
2282     //
2283     if ( fPH2d )
2284         return fPH2d;
2285
2286     // Some parameters
2287     fTimeMax = nbtimebin;
2288     fSf      = samplefrequency;
2289
2290     /*
2291       AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
2292       ,fCalibraMode->GetNz(1)
2293       ,fCalibraMode->GetNrphi(1)));
2294
2295       // Calcul the number of Xbins
2296       Int_t Ntotal1 = 0;
2297       fCalibraMode->ModePadCalibration(2,1);
2298       fCalibraMode->ModePadFragmentation(0,2,0,1);
2299       fCalibraMode->SetDetChamb2(1);
2300       Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
2301       fCalibraMode->ModePadCalibration(0,1);
2302       fCalibraMode->ModePadFragmentation(0,0,0,1);
2303       fCalibraMode->SetDetChamb0(1);
2304       Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
2305       AliInfo(Form("Total number of Xbins: %d",Ntotal1));
2306
2307       CreatePH2d(Ntotal1);
2308     */
2309
2310     CreatePH2d(540);
2311
2312     return fPH2d;
2313 }
2314 //_____________________________________________________________________
2315 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
2316 {
2317     //
2318     // return pointer to TLinearFitter Calibration
2319     // if force is true create a new TLinearFitter if it doesn't exist allready
2320     //
2321     if ( !force || fLinearFitterArray.UncheckedAt(detector) )
2322         return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
2323
2324     // if we are forced and TLinearFitter doesn't yet exist create it
2325
2326     // new TLinearFitter
2327     TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
2328     fLinearFitterArray.AddAt(linearfitter,detector);
2329     return linearfitter;
2330 }
2331
2332 //_____________________________________________________________________
2333 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
2334 {
2335   //
2336   // FillCH2d: Marian style
2337   //
2338
2339   //skip simply the value out of range
2340   if((y>=300.0) || (y<0.0)) return;
2341
2342   //Calcul the y place
2343   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
2344   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
2345
2346   //Fill
2347   fCH2d->GetArray()[place]++;
2348
2349 }
2350
2351 //____________________________________________________________________________
2352 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
2353 {
2354   //
2355   // Analyse array of linear fitter because can not be written
2356   // Store two arrays: one with the param the other one with the error param + number of entries
2357   //
2358
2359   for(Int_t k = 0; k < 540; k++){
2360     TLinearFitter *linearfitter = GetLinearFitter(k);
2361     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
2362       TVectorD  *par  = new TVectorD(2);
2363       TVectorD   pare = TVectorD(2);
2364       TVectorD  *parE = new TVectorD(3);
2365       linearfitter->Eval();
2366       linearfitter->GetParameters(*par);
2367       linearfitter->GetErrors(pare);
2368       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
2369       (*parE)[0] = pare[0]*ppointError;
2370       (*parE)[1] = pare[1]*ppointError;
2371       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
2372       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
2373       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
2374     }
2375   }
2376 }