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