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