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 <TProfile2D.h>
  37 #include <TProfile.h>
  38 #include <TFile.h>
  39 #include <TStyle.h>
  40 #include <TCanvas.h>
  41 #include <TObjArray.h>
  42 #include <TObject.h>
  43 #include <TH1F.h>
  44 #include <TH2I.h>
  45 #include <TH2.h>
  46 #include <TStopwatch.h>
  47 #include <TMath.h>
  48 #include <TDirectory.h>
  49 #include <TTreeStream.h>
  50 #include <TVectorD.h>
  51
  52 #include "AliLog.h"
  53
  54 #include "AliTRDCalibraFillHisto.h"
  55 #include "AliTRDCalibraMode.h"
  56 #include "AliTRDCalibraVector.h"
  57 #include "AliTRDCalibraVdriftLinearFit.h"
  58 #include "AliTRDcalibDB.h"
  59 #include "AliTRDCommonParam.h"
  60 #include "AliTRDpadPlane.h"
  61 #include "AliTRDcluster.h"
  62 #include "AliTRDtrack.h"
  63 #include "AliTRDtrackV1.h"
  64 #include "AliTRDrawStreamBase.h"
  65 #include "AliRawReader.h"
  66 #include "AliRawReaderDate.h"
  67 #include "AliTRDgeometry.h"
  68 #include "./Cal/AliTRDCalROC.h"
  69 #include "./Cal/AliTRDCalDet.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 //______________________________________________________________________________________
 118 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto()
 119   :TObject()
 120   ,fGeo(0)
 121   ,fIsHLT(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   ,fLimitChargeIntegration(kFALSE)
 133   ,fFillWithZero(kFALSE)
 134   ,fNormalizeNbOfCluster(kFALSE)
 135   ,fMaxCluster(0)
 136   ,fNbMaxCluster(0)
 137   ,fCalibraMode(new AliTRDCalibraMode())
 138   ,fDebugStreamer(0)
 139   ,fDebugLevel(0)
 140   ,fDetectorPreviousTrack(-1)
 141   ,fMCMPrevious(-1)
 142   ,fROBPrevious(-1)
 143   ,fNumberClusters(1)
 144   ,fNumberClustersf(30)
 145   ,fProcent(6.0)
 146   ,fDifference(17)
 147   ,fNumberTrack(0)
 148   ,fTimeMax(0)
 149   ,fSf(10.0)
 150   ,fNumberBinCharge(100)
 151   ,fNumberBinPRF(40)
 152   ,fNgroupprf(0)
 153   ,fAmpTotal(0x0)
 154   ,fPHPlace(0x0)
 155   ,fPHValue(0x0)
 156   ,fGoodTracklet(kTRUE)
 157   ,fLinearFitterTracklet(0x0)
 158   ,fEntriesCH(0x0)
 159   ,fEntriesLinearFitter(0x0)
 160   ,fCalibraVector(0x0)
 161   ,fPH2d(0x0)
 162   ,fPRF2d(0x0)
 163   ,fCH2d(0x0)
 164   ,fLinearFitterArray(540)
 165   ,fLinearVdriftFit(0x0)
 166   ,fCalDetGain(0x0)
 167   ,fCalROCGain(0x0)
 168 {
 169   //
 170   // Default constructor
 171   //
 172
 173   //
 174   // Init some default values
 175   //
 176
 177   fNumberUsedCh[0]       = 0;
 178   fNumberUsedCh[1]       = 0;
 179   fNumberUsedPh[0]       = 0;
 180   fNumberUsedPh[1]       = 0;
 181
 182   fGeo = new AliTRDgeometry();
 183
 184 }
 185
 186 //______________________________________________________________________________________
 187 AliTRDCalibraFillHisto::AliTRDCalibraFillHisto(const AliTRDCalibraFillHisto &c)
 188   :TObject(c)
 189   ,fGeo(0)
 190   ,fIsHLT(c.fIsHLT)
 191   ,fMcmCorrectAngle(c.fMcmCorrectAngle)
 192   ,fCH2dOn(c.fCH2dOn)
 193   ,fPH2dOn(c.fPH2dOn)
 194   ,fPRF2dOn(c.fPRF2dOn)
 195   ,fHisto2d(c.fHisto2d)
 196   ,fVector2d(c.fVector2d)
 197   ,fLinearFitterOn(c.fLinearFitterOn)
 198   ,fLinearFitterDebugOn(c.fLinearFitterDebugOn)
 199   ,fRelativeScale(c.fRelativeScale)
 200   ,fThresholdClusterPRF2(c.fThresholdClusterPRF2)
 201   ,fLimitChargeIntegration(c.fLimitChargeIntegration)
 202   ,fFillWithZero(c.fFillWithZero)
 203   ,fNormalizeNbOfCluster(c.fNormalizeNbOfCluster)
 204   ,fMaxCluster(c.fMaxCluster)
 205   ,fNbMaxCluster(c.fNbMaxCluster)
 206   ,fCalibraMode(0x0)
 207   ,fDebugStreamer(0)
 208   ,fDebugLevel(c.fDebugLevel)
 209   ,fDetectorPreviousTrack(c.fDetectorPreviousTrack)
 210   ,fMCMPrevious(c.fMCMPrevious)
 211   ,fROBPrevious(c.fROBPrevious)
 212   ,fNumberClusters(c.fNumberClusters)
 213   ,fNumberClustersf(c.fNumberClustersf)
 214   ,fProcent(c.fProcent)
 215   ,fDifference(c.fDifference)
 216   ,fNumberTrack(c.fNumberTrack)
 217   ,fTimeMax(c.fTimeMax)
 218   ,fSf(c.fSf)
 219   ,fNumberBinCharge(c.fNumberBinCharge)
 220   ,fNumberBinPRF(c.fNumberBinPRF)
 221   ,fNgroupprf(c.fNgroupprf)
 222   ,fAmpTotal(0x0)
 223   ,fPHPlace(0x0)
 224   ,fPHValue(0x0)
 225   ,fGoodTracklet(c.fGoodTracklet)
 226   ,fLinearFitterTracklet(0x0)
 227   ,fEntriesCH(0x0)
 228   ,fEntriesLinearFitter(0x0)
 229   ,fCalibraVector(0x0)
 230   ,fPH2d(0x0)
 231   ,fPRF2d(0x0)
 232   ,fCH2d(0x0)
 233   ,fLinearFitterArray(540)
 234   ,fLinearVdriftFit(0x0)
 235   ,fCalDetGain(0x0)
 236   ,fCalROCGain(0x0)
 237 {
 238   //
 239   // Copy constructor
 240   //
 241   if(c.fCalibraMode)   fCalibraMode = new AliTRDCalibraMode(*c.fCalibraMode);
 242   if(c.fCalibraVector) fCalibraVector = new AliTRDCalibraVector(*c.fCalibraVector);
 243   if(c.fPH2d) {
 244     fPH2d = new TProfile2D(*c.fPH2d);
 245     fPH2d->SetDirectory(0);
 246   }
 247   if(c.fPRF2d) {
 248     fPRF2d = new TProfile2D(*c.fPRF2d);
 249     fPRF2d->SetDirectory(0);
 250   }
 251   if(c.fCH2d) {
 252     fCH2d = new TH2I(*c.fCH2d);
 253     fCH2d->SetDirectory(0);
 254   }
 255   if(c.fLinearVdriftFit){
 256     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit(*c.fLinearVdriftFit);
 257   }
 258
 259   if(c.fCalDetGain)  fCalDetGain   = new AliTRDCalDet(*c.fCalDetGain);
 260   if(c.fCalROCGain)  fCalROCGain   = new AliTRDCalROC(*c.fCalROCGain);
 261
 262   if (fGeo) {
 263     delete fGeo;
 264   }
 265   fGeo = new AliTRDgeometry();
 266 }
 267
 268 //____________________________________________________________________________________
 269 AliTRDCalibraFillHisto::~AliTRDCalibraFillHisto()
 270 {
 271   //
 272   // AliTRDCalibraFillHisto destructor
 273   //
 274
 275   ClearHistos();
 276   if ( fDebugStreamer ) delete fDebugStreamer;
 277
 278   if ( fCalDetGain )  delete fCalDetGain;
 279   if ( fCalROCGain )  delete fCalROCGain;
 280
 281   if( fLinearFitterTracklet ) { delete fLinearFitterTracklet; }
 282
 283   delete [] fPHPlace;
 284   delete [] fPHValue;
 285   delete [] fEntriesCH;
 286   delete [] fEntriesLinearFitter;
 287   delete [] fAmpTotal;
 288
 289   for(Int_t idet=0; idet<AliTRDgeometry::kNdet; idet++){
 290     TLinearFitter *f = (TLinearFitter*)fLinearFitterArray.At(idet);
 291     if(f) { delete f;}
 292   }
 293   if (fGeo) {
 294     delete fGeo;
 295   }
 296
 297 }
 298 //_____________________________________________________________________________
 299 void AliTRDCalibraFillHisto::Destroy()
 300 {
 301   //
 302   // Delete instance
 303   //
 304
 305   if (fgInstance) {
 306     delete fgInstance;
 307     fgInstance = 0x0;
 308   }
 309 }
 310 //_____________________________________________________________________________
 311 void AliTRDCalibraFillHisto::DestroyDebugStreamer()
 312 {
 313   //
 314   // Delete DebugStreamer
 315   //
 316
 317   if ( fDebugStreamer ) delete fDebugStreamer;
 318   fDebugStreamer = 0x0;
 319
 320 }
 321 //_____________________________________________________________________________
 322 void AliTRDCalibraFillHisto::ClearHistos()
 323 {
 324   //
 325   // Delete the histos
 326   //
 327
 328   if (fPH2d) {
 329     delete fPH2d;
 330     fPH2d  = 0x0;
 331   }
 332   if (fCH2d) {
 333     delete fCH2d;
 334     fCH2d  = 0x0;
 335   }
 336   if (fPRF2d) {
 337     delete fPRF2d;
 338     fPRF2d = 0x0;
 339   }
 340
 341 }
 342 //////////////////////////////////////////////////////////////////////////////////
 343 // calibration with AliTRDtrackV1: Init, Update
 344 //////////////////////////////////////////////////////////////////////////////////
 345 //____________Functions for initialising the AliTRDCalibraFillHisto in the code_________
 346 Bool_t AliTRDCalibraFillHisto::Init2Dhistos()
 347 {
 348   //
 349   // Init the histograms and stuff to be filled
 350   //
 351
 352   // DB Setting
 353   // Get cal
 354   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 355   if (!cal) {
 356     AliInfo("Could not get calibDB");
 357     return kFALSE;
 358   }
 359   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
 360   if (!parCom) {
 361     AliInfo("Could not get CommonParam");
 362     return kFALSE;
 363   }
 364
 365   // Some parameters
 366   fTimeMax            = cal->GetNumberOfTimeBins();
 367   fSf                 = parCom->GetSamplingFrequency();
 368   if(!fNormalizeNbOfCluster) fRelativeScale = 20.0;
 369   else fRelativeScale = 1.18;
 370   fNumberClustersf    = fTimeMax;
 371   fNumberClusters     = (Int_t)(0.5*fTimeMax);
 372
 373   // Init linear fitter
 374   if(!fLinearFitterTracklet) {
 375     fLinearFitterTracklet = new TLinearFitter(2,"pol1");
 376     fLinearFitterTracklet->StoreData(kTRUE);
 377   }
 378
 379   //calib object from database used for reconstruction
 380   if( fCalDetGain ){
 381     fCalDetGain->~AliTRDCalDet();
 382     new(fCalDetGain) AliTRDCalDet(*(cal->GetGainFactorDet()));
 383   }else fCalDetGain = new AliTRDCalDet(*(cal->GetGainFactorDet()));
 384
 385   // Calcul Xbins Chambd0, Chamb2
 386   Int_t ntotal0 = CalculateTotalNumberOfBins(0);
 387   Int_t ntotal1 = CalculateTotalNumberOfBins(1);
 388   Int_t ntotal2 = CalculateTotalNumberOfBins(2);
 389
 390   // If vector method On initialised all the stuff
 391   if(fVector2d){
 392     fCalibraVector = new AliTRDCalibraVector();
 393     fCalibraVector->SetNumberBinCharge(fNumberBinCharge);
 394     fCalibraVector->SetTimeMax(fTimeMax);
 395     if(fNgroupprf != 0) {
 396       fCalibraVector->SetNumberBinPRF(2*fNgroupprf*fNumberBinPRF);
 397       fCalibraVector->SetPRFRange((Float_t)(3.0*fNgroupprf));
 398     }
 399     else {
 400       fCalibraVector->SetNumberBinPRF(fNumberBinPRF);
 401       fCalibraVector->SetPRFRange(1.5);
 402     }
 403     for(Int_t k = 0; k < 3; k++){
 404       fCalibraVector->SetDetCha0(k,fCalibraMode->GetDetChamb0(k));
 405       fCalibraVector->SetDetCha2(k,fCalibraMode->GetDetChamb2(k));
 406     }
 407     TString namech("Nz");
 408     namech += fCalibraMode->GetNz(0);
 409     namech += "Nrphi";
 410     namech += fCalibraMode->GetNrphi(0);
 411     fCalibraVector->SetNameCH((const char* ) namech);
 412     TString nameph("Nz");
 413     nameph += fCalibraMode->GetNz(1);
 414     nameph += "Nrphi";
 415     nameph += fCalibraMode->GetNrphi(1);
 416     fCalibraVector->SetNamePH((const char* ) nameph);
 417     TString nameprf("Nz");
 418     nameprf += fCalibraMode->GetNz(2);
 419     nameprf += "Nrphi";
 420     nameprf += fCalibraMode->GetNrphi(2);
 421     nameprf += "Ngp";
 422     nameprf += fNgroupprf;
 423     fCalibraVector->SetNamePRF((const char* ) nameprf);
 424   }
 425
 426   // Create the 2D histos corresponding to the pad groupCalibration mode
 427   if (fCH2dOn) {
 428
 429     AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
 430                 ,fCalibraMode->GetNz(0)
 431                 ,fCalibraMode->GetNrphi(0)));
 432
 433     // Create the 2D histo
 434     if (fHisto2d) {
 435       CreateCH2d(ntotal0);
 436     }
 437     // Variable
 438     fAmpTotal = new Float_t[TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0))];
 439     for (Int_t k = 0; k < TMath::Max(fCalibraMode->GetDetChamb2(0),fCalibraMode->GetDetChamb0(0)); k++) {
 440       fAmpTotal[k] = 0.0;
 441     }
 442     //Statistics
 443     fEntriesCH = new Int_t[ntotal0];
 444     for(Int_t k = 0; k < ntotal0; k++){
 445       fEntriesCH[k] = 0;
 446     }
 447
 448   }
 449   if (fPH2dOn) {
 450
 451     AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
 452                 ,fCalibraMode->GetNz(1)
 453                 ,fCalibraMode->GetNrphi(1)));
 454
 455     // Create the 2D histo
 456     if (fHisto2d) {
 457       CreatePH2d(ntotal1);
 458     }
 459     // Variable
 460     fPHPlace = new Short_t[fTimeMax];
 461     for (Int_t k = 0; k < fTimeMax; k++) {
 462       fPHPlace[k] = -1;
 463     }
 464     fPHValue = new Float_t[fTimeMax];
 465     for (Int_t k = 0; k < fTimeMax; k++) {
 466       fPHValue[k] = 0.0;
 467     }
 468   }
 469   if (fLinearFitterOn) {
 470     //fLinearFitterArray.Expand(540);
 471     fLinearFitterArray.SetName("ArrayLinearFitters");
 472     fEntriesLinearFitter = new Int_t[540];
 473     for(Int_t k = 0; k < 540; k++){
 474       fEntriesLinearFitter[k] = 0;
 475     }
 476     fLinearVdriftFit = new AliTRDCalibraVdriftLinearFit();
 477   }
 478
 479   if (fPRF2dOn) {
 480
 481     AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
 482                 ,fCalibraMode->GetNz(2)
 483                 ,fCalibraMode->GetNrphi(2)));
 484     // Create the 2D histo
 485     if (fHisto2d) {
 486       CreatePRF2d(ntotal2);
 487     }
 488   }
 489
 490   return kTRUE;
 491
 492 }
 493 //____________Offline tracking in the AliTRDtracker____________________________
 494 Bool_t AliTRDCalibraFillHisto::UpdateHistograms(AliTRDtrack *t)
 495 {
 496   //
 497   // Use AliTRDtrack for the calibration
 498   //
 499
 500
 501   AliTRDcluster *cl = 0x0;       // pointeur to cluster
 502   Int_t index0 = 0;              // index of the first cluster in the current chamber
 503   Int_t index1 = 0;              // index of the current cluster in the current chamber
 504
 505   //////////////////////////////
 506   // loop over the clusters
 507   ///////////////////////////////
 508   while((cl = t->GetCluster(index1))){
 509
 510     /////////////////////////////////////////////////////////////////////////
 511     // Fill the infos for the previous clusters if not the same detector
 512     ////////////////////////////////////////////////////////////////////////
 513     Int_t detector = cl->GetDetector();
 514     if ((detector != fDetectorPreviousTrack) &&
 515         (index0 != index1)) {
 516
 517       fNumberTrack++;
 518
 519       //If the same track, then look if the previous detector is in
 520       //the same plane, if yes: not a good track
 521       if ((GetLayer(detector) == GetLayer(fDetectorPreviousTrack))) {
 522         return kFALSE;
 523       }
 524
 525       // Fill only if the track doesn't touch a masked pad or doesn't
 526       if (fGoodTracklet) {
 527
 528         // drift velocity unables to cut bad tracklets
 529         Bool_t  pass = FindP1TrackPHtracklet(t,index0,index1);
 530
 531         // Gain calibration
 532         if (fCH2dOn) {
 533           FillTheInfoOfTheTrackCH(index1-index0);
 534         }
 535
 536         // PH calibration
 537         if (fPH2dOn) {
 538           FillTheInfoOfTheTrackPH();
 539         }
 540
 541         if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
 542
 543
 544       } // if a good tracklet
 545
 546       // reset stuff
 547       ResetfVariablestracklet();
 548       index0 = index1;
 549
 550     } // Fill at the end the charge
 551
 552
 553     /////////////////////////////////////////////////////////////
 554     // Localise and take the calib gain object for the detector
 555     ////////////////////////////////////////////////////////////
 556     if (detector != fDetectorPreviousTrack) {
 557
 558       //Localise the detector
 559       LocalisationDetectorXbins(detector);
 560
 561       // Get cal
 562       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 563       if (!cal) {
 564         AliInfo("Could not get calibDB");
 565         return kFALSE;
 566       }
 567
 568       // Get calib objects
 569       if( fCalROCGain ){
 570         if(!fIsHLT){
 571           fCalROCGain->~AliTRDCalROC();
 572           new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
 573         }
 574       }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
 575
 576     }
 577
 578     // Reset the detectbjobsor
 579     fDetectorPreviousTrack = detector;
 580
 581     ///////////////////////////////////////
 582     // Store the info of the cluster
 583     ///////////////////////////////////////
 584     Int_t row = cl->GetPadRow();
 585     Int_t col = cl->GetPadCol();
 586     CheckGoodTrackletV1(cl);
 587     Int_t     group[2] = {0,0};
 588     if(fCH2dOn)  group[0]  = CalculateCalibrationGroup(0,row,col);
 589     if(fPH2dOn)  group[1]  = CalculateCalibrationGroup(1,row,col);
 590     StoreInfoCHPHtrack(cl,t->GetClusterdQdl(index1),group,row,col);
 591
 592     index1++;
 593
 594   } // while on clusters
 595
 596   ///////////////////////////
 597   // Fill the last plane
 598   //////////////////////////
 599   if( index0 != index1 ){
 600
 601     fNumberTrack++;
 602
 603     if (fGoodTracklet) {
 604
 605       // drift velocity unables to cut bad tracklets
 606       Bool_t  pass = FindP1TrackPHtracklet(t,index0,index1);
 607
 608       // Gain calibration
 609       if (fCH2dOn) {
 610         FillTheInfoOfTheTrackCH(index1-index0);
 611       }
 612
 613       // PH calibration
 614       if (fPH2dOn) {
 615         FillTheInfoOfTheTrackPH();
 616       }
 617
 618       if(pass && fPRF2dOn) HandlePRFtracklet(t,index0,index1);
 619
 620     } // if a good tracklet
 621
 622   }
 623
 624   // reset stuff
 625   ResetfVariablestracklet();
 626
 627   return kTRUE;
 628
 629 }
 630 //____________Offline tracking in the AliTRDtracker____________________________
 631 Bool_t AliTRDCalibraFillHisto::UpdateHistogramsV1(AliTRDtrackV1 *t)
 632 {
 633   //
 634   // Use AliTRDtrackV1 for the calibration
 635   //
 636
 637
 638   const AliTRDseedV1 *tracklet = 0x0;          // tracklet per plane
 639   AliTRDcluster *cl      = 0x0;                // cluster attached now to the tracklet
 640   Bool_t         newtr   = kTRUE;              // new track
 641
 642   // Get cal
 643   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
 644   if (!cal) {
 645     AliInfo("Could not get calibDB");
 646     return kFALSE;
 647   }
 648
 649   ///////////////////////////
 650   // loop over the tracklet
 651   ///////////////////////////
 652   for(Int_t itr = 0; itr < 6; itr++){
 653
 654     if(!(tracklet = t->GetTracklet(itr))) continue;
 655     if(!tracklet->IsOK()) continue;
 656     fNumberTrack++;
 657     ResetfVariablestracklet();
 658
 659     //////////////////////////////////////////
 660     // localisation of the tracklet and dqdl
 661     //////////////////////////////////////////
 662     Int_t layer    = tracklet->GetPlane();
 663     Int_t ic = 0;
 664     while(!(cl = tracklet->GetClusters(ic++))) continue;
 665     Int_t detector = cl->GetDetector();
 666     if (detector != fDetectorPreviousTrack) {
 667       // if not a new track
 668       if(!newtr){
 669         // don't use the rest of this track if in the same plane
 670         if (layer == GetLayer(fDetectorPreviousTrack)) {
 671           //printf("bad tracklet, same layer for detector %d\n",detector);
 672           break;
 673         }
 674       }
 675       //Localise the detector bin
 676       LocalisationDetectorXbins(detector);
 677       // Get calib objects
 678       if( fCalROCGain ){
 679         if(!fIsHLT){
 680           fCalROCGain->~AliTRDCalROC();
 681           new(fCalROCGain) AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
 682         }
 683       }else fCalROCGain = new AliTRDCalROC(*(cal->GetGainFactorROC(detector)));
 684
 685       // reset
 686       fDetectorPreviousTrack = detector;
 687     }
 688     newtr = kFALSE;
 689
 690     ////////////////////////////
 691     // loop over the clusters
 692     ////////////////////////////
 693     Int_t nbclusters = 0;
 694     for(int jc=0; jc<AliTRDseed::knTimebins; jc++){
 695       if(!(cl = tracklet->GetClusters(jc))) continue;
 696       nbclusters++;
 697
 698       // Store the info bis of the tracklet
 699       Int_t row = cl->GetPadRow();
 700       Int_t col = cl->GetPadCol();
 701       CheckGoodTrackletV1(cl);
 702       Int_t     group[2] = {0,0};
 703       if(fCH2dOn)  group[0]  = CalculateCalibrationGroup(0,row,col);
 704       if(fPH2dOn)  group[1]  = CalculateCalibrationGroup(1,row,col);
 705       StoreInfoCHPHtrack(cl, tracklet->GetdQdl(jc),group,row,col);
 706     }
 707
 708     ////////////////////////////////////////
 709     // Fill the stuffs if a good tracklet
 710     ////////////////////////////////////////
 711     if (fGoodTracklet) {
 712
 713       // drift velocity unables to cut bad tracklets
 714       Bool_t  pass = FindP1TrackPHtrackletV1(tracklet, nbclusters);
 715
 716       // Gain calibration
 717       if (fCH2dOn) {
 718         FillTheInfoOfTheTrackCH(nbclusters);
 719       }
 720
 721       // PH calibration
 722       if (fPH2dOn) {
 723         FillTheInfoOfTheTrackPH();
 724       }
 725
 726       if(pass && fPRF2dOn) HandlePRFtrackletV1(tracklet,nbclusters);
 727
 728     } // if a good tracklet
 729   }
 730
 731   return kTRUE;
 732
 733 }
 734 ///////////////////////////////////////////////////////////////////////////////////
 735 // Routine inside the update with AliTRDtrack
 736 ///////////////////////////////////////////////////////////////////////////////////
 737 //____________Offine tracking in the AliTRDtracker_____________________________
 738 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
 739 {
 740   //
 741   // Drift velocity calibration:
 742   // Fit the clusters with a straight line
 743   // From the slope find the drift velocity
 744   //
 745
 746   //Number of points: if less than 3 return kFALSE
 747   Int_t npoints = index1-index0;
 748   if(npoints <= 2) return kFALSE;
 749
 750   /////////////////
 751   //Variables
 752   ////////////////
 753   Int_t    detector                   = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
 754   // parameters of the track
 755   Double_t snp                        = t->GetSnpPlane(GetLayer(detector)); // sin angle in the plan yx track
 756   Double_t tgl                        = t->GetTglPlane(GetLayer(detector)); // dz/dl and not dz/dx!
 757   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
 758   if( TMath::Abs(snp) <  1.){
 759     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
 760   }
 761   Float_t dzdx                        = tgl*TMath::Sqrt(1+tnp*tnp);         // dz/dx now from dz/dl
 762   // tilting pad and cross row
 763   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
 764   Int_t    rowp                       = -1;                                 // if it crosses a pad row
 765   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetLayer(detector),GetStack(detector));
 766   Double_t tiltingangle               = padplane->GetTiltingAngle();        // tiltingangle of the pad
 767   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
 768   // linear fit
 769   fLinearFitterTracklet->ClearPoints();
 770   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
 771   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
 772   Double_t pointError                 = 0.0;                                // error after straight line fit
 773   Int_t     nbli                      = 0;                                  // number linear fitter points
 774
 775   //////////////////////////////
 776   // loop over clusters
 777   ////////////////////////////
 778   for(Int_t k = 0; k < npoints; k++){
 779
 780     AliTRDcluster *cl                 = (AliTRDcluster *) t->GetCluster(k+index0);
 781     if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
 782     Double_t ycluster                 = cl->GetY();
 783     Int_t time                        = cl->GetPadTime();
 784     Double_t timeis                   = time/fSf;
 785     //Double_t q                        = cl->GetQ();
 786     //See if cross two pad rows
 787     Int_t    row                      = cl->GetPadRow();
 788     if(k==0) rowp                     = row;
 789     if(row != rowp) crossrow          = 1;
 790
 791     fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
 792     nbli++;
 793
 794   }
 795
 796   //////////////////////////////
 797   // linear fit
 798   //////////////////////////////
 799   if(nbli <= 2){
 800     fLinearFitterTracklet->ClearPoints();
 801     return kFALSE;
 802   }
 803   TVectorD pars;
 804   fLinearFitterTracklet->Eval();
 805   fLinearFitterTracklet->GetParameters(pars);
 806   pointError  =  TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
 807   errorpar    =  fLinearFitterTracklet->GetParError(1)*pointError;
 808   dydt        = pars[1];
 809   fLinearFitterTracklet->ClearPoints();
 810
 811   /////////////////////////////
 812   // debug
 813   ////////////////////////////
 814   if(fDebugLevel > 0){
 815     if ( !fDebugStreamer ) {
 816       //debug stream
 817       TDirectory *backup = gDirectory;
 818       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
 819       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 820     }
 821
 822
 823     (* fDebugStreamer) << "FindP1TrackPHtracklet0"<<
 824       //"snpright="<<snpright<<
 825       "npoints="<<npoints<<
 826       "nbli="<<nbli<<
 827       "detector="<<detector<<
 828       "snp="<<snp<<
 829       "tnp="<<tnp<<
 830       "tgl="<<tgl<<
 831       "tnt="<<tnt<<
 832       "dydt="<<dydt<<
 833       "dzdx="<<dzdx<<
 834       "crossrow="<<crossrow<<
 835       "errorpar="<<errorpar<<
 836       "pointError="<<pointError<<
 837       "\n";
 838
 839
 840     Int_t nbclusters = index1-index0;
 841     Int_t layer      = GetLayer(fDetectorPreviousTrack);
 842
 843     (* fDebugStreamer) << "FindP1TrackPHtracklet1"<<
 844       //"snpright="<<snpright<<
 845       "nbclusters="<<nbclusters<<
 846       "detector="<<fDetectorPreviousTrack<<
 847       "layer="<<layer<<
 848       "\n";
 849
 850   }
 851
 852   ///////////////////////////
 853   // quality cuts
 854   ///////////////////////////
 855   if(npoints < fNumberClusters) return kFALSE;
 856   if(npoints > fNumberClustersf) return kFALSE;
 857   if(pointError >= 0.1) return kFALSE;
 858   if(crossrow == 1) return kFALSE;
 859
 860   ////////////////////////////
 861   // fill
 862   ////////////////////////////
 863   if(fLinearFitterOn){
 864     //Add to the linear fitter of the detector
 865     if( TMath::Abs(snp) <  1.){
 866       Double_t x = tnp-dzdx*tnt;
 867       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
 868       if(fLinearFitterDebugOn) {
 869         fLinearVdriftFit->Update(detector,x,pars[1]);
 870       }
 871       fEntriesLinearFitter[detector]++;
 872     }
 873   }
 874
 875   return kTRUE;
 876 }
 877 //____________Offine tracking in the AliTRDtracker_____________________________
 878 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
 879 {
 880   //
 881   // Drift velocity calibration:
 882   // Fit the clusters with a straight line
 883   // From the slope find the drift velocity
 884   //
 885
 886   ////////////////////////////////////////////////
 887   //Number of points: if less than 3 return kFALSE
 888   /////////////////////////////////////////////////
 889   if(nbclusters <= 2) return kFALSE;
 890
 891   ////////////
 892   //Variables
 893   ////////////
 894   // results of the linear fit
 895   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
 896   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
 897   Double_t pointError                 = 0.0;                                // error after straight line fit
 898   // pad row problemes: avoid tracklet that cross pad rows, tilting angle in the constant
 899   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
 900   Int_t    rowp                       = -1;                                 // if it crosses a pad row
 901   Float_t  tnt                        = tracklet->GetTilt();                // tan tiltingangle
 902   fLinearFitterTracklet->ClearPoints();
 903
 904
 905   ///////////////////////////////////////////
 906   // Take the parameters of the track
 907   //////////////////////////////////////////
 908   // take now the snp, tnp and tgl from the track
 909   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
 910   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
 911   if( TMath::Abs(snp) <  1.){
 912     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
 913   }
 914   Double_t tgl  = tracklet->GetTgl();           // dz/dl
 915   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
 916   // at the entrance
 917   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
 918   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
 919   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
 920   // at the end with correction due to linear fit
 921   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
 922   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
 923
 924
 925   ////////////////////////////
 926   // loop over the clusters
 927   ////////////////////////////
 928   Int_t  nbli = 0;
 929   AliTRDcluster *cl                   = 0x0;
 930   for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
 931     if(!(cl = tracklet->GetClusters(ic))) continue;
 932     if((fLimitChargeIntegration) && (!cl->IsInChamber())) continue;
 933
 934     Double_t ycluster                 = cl->GetY();
 935     Int_t time                        = cl->GetPadTime();
 936     Double_t timeis                   = time/fSf;
 937     //See if cross two pad rows
 938     Int_t    row                      = cl->GetPadRow();
 939     if(rowp==-1) rowp                 = row;
 940     if(row != rowp) crossrow          = 1;
 941
 942     fLinearFitterTracklet->AddPoint(&timeis,ycluster,1);
 943     nbli++;
 944
 945   }
 946
 947   ////////////////////////////////////
 948   // Do the straight line fit now
 949   ///////////////////////////////////
 950   if(nbli <= 2){
 951     fLinearFitterTracklet->ClearPoints();
 952     return kFALSE;
 953   }
 954   TVectorD pars;
 955   fLinearFitterTracklet->Eval();
 956   fLinearFitterTracklet->GetParameters(pars);
 957   pointError  =  TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
 958   errorpar    =  fLinearFitterTracklet->GetParError(1)*pointError;
 959   dydt        = pars[1];
 960   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
 961   fLinearFitterTracklet->ClearPoints();
 962
 963   ////////////////////////////////
 964   // Debug stuff
 965   ///////////////////////////////
 966
 967
 968   if(fDebugLevel > 0){
 969     if ( !fDebugStreamer ) {
 970       //debug stream
 971       TDirectory *backup = gDirectory;
 972       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
 973       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 974     }
 975
 976
 977     Int_t layer = GetLayer(fDetectorPreviousTrack);
 978
 979     (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
 980       //"snpright="<<snpright<<
 981       "nbli="<<nbli<<
 982       "nbclusters="<<nbclusters<<
 983       "detector="<<fDetectorPreviousTrack<<
 984       "layer="<<layer<<
 985       "snp="<<snp<<
 986       "tnp="<<tnp<<
 987       "tgl="<<tgl<<
 988       "tnt="<<tnt<<
 989       "dydt="<<dydt<<
 990       "dzdx="<<dzdx<<
 991       "crossrow="<<crossrow<<
 992       "errorpar="<<errorpar<<
 993       "pointError="<<pointError<<
 994       "\n";
 995
 996   }
 997
 998   /////////////////////////
 999   // Cuts quality
1000   ////////////////////////
1001
1002   if(nbclusters < fNumberClusters) return kFALSE;
1003   if(nbclusters > fNumberClustersf) return kFALSE;
1004   if(pointError >= 0.3) return kFALSE;
1005   if(crossrow == 1) return kFALSE;
1006
1007   ///////////////////////
1008   // Fill
1009   //////////////////////
1010
1011   if(fLinearFitterOn){
1012     //Add to the linear fitter of the detector
1013     if( TMath::Abs(snp) <  1.){
1014       Double_t x = tnp-dzdx*tnt;
1015       (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1016       if(fLinearFitterDebugOn) {
1017         fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1018       }
1019       fEntriesLinearFitter[fDetectorPreviousTrack]++;
1020     }
1021   }
1022
1023   return kTRUE;
1024 }
1025 //____________Offine tracking in the AliTRDtracker_____________________________
1026 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1027 {
1028   //
1029   // PRF width calibration
1030   // Assume a Gaussian shape: determinate the position of the three pad clusters
1031   // Fit with a straight line
1032   // Take the fitted values for all the clusters (3 or 2 pad clusters)
1033   // Fill the PRF as function of angle of the track
1034   //
1035   //
1036
1037
1038   //////////////////////////
1039   // variables
1040   /////////////////////////
1041   Int_t npoints  = index1-index0;                                           // number of total points
1042   Int_t nb3pc    = 0;                                                       // number of three pads clusters used for fit
1043   Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1044   // To see the difference due to the fit
1045   Double_t *padPositions;
1046   padPositions = new Double_t[npoints];
1047   for(Int_t k = 0; k < npoints; k++){
1048     padPositions[k] = 0.0;
1049   }
1050   // Take the tgl and snp with the track t now
1051   Double_t  tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1052   Double_t  snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1053   Float_t  dzdx = 0.0;                                // dzdx
1054   Float_t  tnp  = 0.0;
1055   if(TMath::Abs(snp) < 1.0){
1056     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1057     dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1058   }
1059   // linear fitter
1060   fLinearFitterTracklet->ClearPoints();
1061
1062   ///////////////////////////
1063   // calcul the tnp group
1064   ///////////////////////////
1065   Bool_t echec   = kFALSE;
1066   Double_t shift = 0.0;
1067   //Calculate the shift in x coresponding to this tnp
1068   if(fNgroupprf != 0.0){
1069     shift      = -3.0*(fNgroupprf-1)-1.5;
1070     Double_t limithigh  = -0.2*(fNgroupprf-1);
1071     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1072     else{
1073       while(tnp > limithigh){
1074         limithigh += 0.2;
1075         shift += 3.0;
1076       }
1077     }
1078   }
1079   if(echec) {
1080     delete [] padPositions;
1081     return kFALSE;
1082   }
1083
1084   //////////////////////
1085   // loop clusters
1086   /////////////////////
1087   for(Int_t k = 0;  k < npoints; k++){
1088     //Take the cluster
1089     AliTRDcluster *cl  = (AliTRDcluster *) t->GetCluster(k+index0);
1090     Short_t  *signals  = cl->GetSignals();
1091     Double_t     time  = cl->GetLocalTimeBin();
1092     //Calculate x if possible
1093     Float_t xcenter    = 0.0;
1094     Bool_t  echec1      = kTRUE;
1095     if((time<=7) || (time>=21)) continue;
1096     // Center 3 balanced: position with the center of the pad
1097     if ((((Float_t) signals[3]) > 0.0) &&
1098         (((Float_t) signals[2]) > 0.0) &&
1099         (((Float_t) signals[4]) > 0.0)) {
1100       echec1 = kFALSE;
1101       // Security if the denomiateur is 0
1102       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1103            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1104         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1105           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1106                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1107       }
1108       else {
1109         echec1 = kTRUE;
1110       }
1111     }
1112     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1113     if(echec) continue;
1114     //if no echec: calculate with the position of the pad
1115     // Position of the cluster
1116     Double_t       padPosition = xcenter +  cl->GetPadCol();
1117     padPositions[k]            = padPosition;
1118     nb3pc++;
1119     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1120   }//clusters loop
1121
1122
1123   /////////////////////////////
1124   // fit
1125   ////////////////////////////
1126   if(nb3pc < 3) {
1127     delete [] padPositions;
1128     fLinearFitterTracklet->ClearPoints();
1129     return kFALSE;
1130   }
1131   fLinearFitterTracklet->Eval();
1132   TVectorD line(2);
1133   fLinearFitterTracklet->GetParameters(line);
1134   Float_t  pointError  = -1.0;
1135   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1136     pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1137   }
1138   fLinearFitterTracklet->ClearPoints();
1139
1140   /////////////////////////////////////////////////////
1141   // Now fill the PRF: second loop over clusters
1142   /////////////////////////////////////////////////////
1143   for(Int_t k = 0;  k < npoints; k++){
1144     //Take the cluster
1145     AliTRDcluster *cl      = (AliTRDcluster *) t->GetCluster(k+index0);
1146     Short_t  *signals      = cl->GetSignals();              // signal
1147     Double_t     time      = cl->GetLocalTimeBin();              // time bin
1148     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1149     Float_t padPos         = cl->GetPadCol();               // middle pad
1150     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1151     Float_t ycenter        = 0.0;                           // relative center charge
1152     Float_t ymin           = 0.0;                           // relative left charge
1153     Float_t ymax           = 0.0;                           // relative right charge
1154
1155     //Requiere simply two pads clusters at least
1156     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1157        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1158       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1159       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1160       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1161       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1162     }
1163
1164     //calibration group
1165     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1166     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1167     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1168     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1169     Float_t   xcl          = cl->GetY();                                // y cluster
1170     Float_t   qcl          = cl->GetQ();                                // charge cluster
1171     Int_t     layer        = GetLayer(detector);                        // layer
1172     Int_t     stack        = GetStack(detector);                        // stack
1173     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1174     Double_t  x            = dpad;                                      // reconstructed position moved
1175     Float_t   ep           = pointError;                                // error of fit
1176     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1177     Float_t   signal3      = (Float_t)signals[3];                       // signal
1178     Float_t   signal2      = (Float_t)signals[2];                       // signal
1179     Float_t   signal4      = (Float_t)signals[4];                       // signal
1180     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1181
1182     //////////////////////////////
1183     // debug
1184     /////////////////////////////
1185     if(fDebugLevel > 0){
1186       if ( !fDebugStreamer ) {
1187         //debug stream
1188         TDirectory *backup = gDirectory;
1189         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1190         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1191       }
1192
1193
1194       x = xdiff;
1195       Int_t type=0;
1196       Float_t y = ycenter;
1197       (* fDebugStreamer) << "HandlePRFtracklet"<<
1198         "caligroup="<<caligroup<<
1199         "detector="<<detector<<
1200         "layer="<<layer<<
1201         "stack="<< stack <<
1202         "npoints="<<npoints<<
1203         "Np="<<nb3pc<<
1204         "ep="<<ep<<
1205         "type="<<type<<
1206         "snp="<<snp<<
1207         "tnp="<<tnp<<
1208         "tgl="<<tgl<<
1209         "dzdx="<<dzdx<<
1210         "padPos="<<padPos<<
1211         "padPosition="<<padPositions[k]<<
1212         "padPosTracklet="<<padPosTracklet<<
1213         "x="<<x<<
1214         "y="<<y<<
1215         "xcl="<<xcl<<
1216         "qcl="<<qcl<<
1217         "signal1="<<signal1<<
1218         "signal2="<<signal2<<
1219         "signal3="<<signal3<<
1220         "signal4="<<signal4<<
1221         "signal5="<<signal5<<
1222         "time="<<time<<
1223         "\n";
1224       x=-(xdiff+1);
1225       y = ymin;
1226       type=-1;
1227       (* fDebugStreamer) << "HandlePRFtracklet"<<
1228         "caligroup="<<caligroup<<
1229         "detector="<<detector<<
1230         "layer="<<layer<<
1231         "stack="<<stack<<
1232         "npoints="<<npoints<<
1233         "Np="<<nb3pc<<
1234         "ep="<<ep<<
1235         "type="<<type<<
1236         "snp="<<snp<<
1237         "tnp="<<tnp<<
1238         "tgl="<<tgl<<
1239         "dzdx="<<dzdx<<
1240         "padPos="<<padPos<<
1241         "padPosition="<<padPositions[k]<<
1242         "padPosTracklet="<<padPosTracklet<<
1243         "x="<<x<<
1244         "y="<<y<<
1245         "xcl="<<xcl<<
1246         "qcl="<<qcl<<
1247         "signal1="<<signal1<<
1248         "signal2="<<signal2<<
1249         "signal3="<<signal3<<
1250         "signal4="<<signal4<<
1251         "signal5="<<signal5<<
1252         "time="<<time<<
1253         "\n";
1254       x=1-xdiff;
1255       y = ymax;
1256       type=1;
1257       (* fDebugStreamer) << "HandlePRFtracklet"<<
1258         "caligroup="<<caligroup<<
1259         "detector="<<detector<<
1260         "layer="<<layer<<
1261         "stack="<<stack<<
1262         "npoints="<<npoints<<
1263         "Np="<<nb3pc<<
1264         "ep="<<ep<<
1265         "type="<<type<<
1266         "snp="<<snp<<
1267         "tnp="<<tnp<<
1268         "tgl="<<tgl<<
1269         "dzdx="<<dzdx<<
1270         "padPos="<<padPos<<
1271         "padPosition="<<padPositions[k]<<
1272         "padPosTracklet="<<padPosTracklet<<
1273         "x="<<x<<
1274         "y="<<y<<
1275         "xcl="<<xcl<<
1276         "qcl="<<qcl<<
1277         "signal1="<<signal1<<
1278         "signal2="<<signal2<<
1279         "signal3="<<signal3<<
1280         "signal4="<<signal4<<
1281         "signal5="<<signal5<<
1282         "time="<<time<<
1283         "\n";
1284
1285     }
1286
1287     ////////////////////////////
1288     // quality cuts
1289     ///////////////////////////
1290     if(npoints < fNumberClusters) continue;
1291     if(npoints > fNumberClustersf) continue;
1292     if(nb3pc <= 5) continue;
1293     if((time >= 21) || (time < 7)) continue;
1294     if(TMath::Abs(snp) >= 1.0) continue;
1295     if(TMath::Abs(qcl) < 80) continue;
1296
1297     ////////////////////////////
1298     // Fill
1299     ///////////////////////////
1300     if (fHisto2d) {
1301       if(TMath::Abs(dpad) < 1.5) {
1302         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1303         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1304       }
1305       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1306         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1307         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1308       }
1309       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1310         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1311         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1312       }
1313     }
1314     if (fVector2d) {
1315       if(TMath::Abs(dpad) < 1.5) {
1316         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1317         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1318       }
1319       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1320         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1321         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1322       }
1323       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1324         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1325         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1326       }
1327     }
1328   }
1329   delete [] padPositions;
1330   return kTRUE;
1331
1332 }
1333 //____________Offine tracking in the AliTRDtracker_____________________________
1334 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1335 {
1336   //
1337   // PRF width calibration
1338   // Assume a Gaussian shape: determinate the position of the three pad clusters
1339   // Fit with a straight line
1340   // Take the fitted values for all the clusters (3 or 2 pad clusters)
1341   // Fill the PRF as function of angle of the track
1342   //
1343   //
1344
1345   //printf("begin\n");
1346   ///////////////////////////////////////////
1347   // Take the parameters of the track
1348   //////////////////////////////////////////
1349   // take now the snp, tnp and tgl from the track
1350   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1351   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1352   if( TMath::Abs(snp) <  1.){
1353     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1354   }
1355   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1356   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1357   // at the entrance
1358   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1359   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1360   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1361   // at the end with correction due to linear fit
1362   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1363   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1364
1365   ///////////////////////////////
1366   // Calculate tnp group shift
1367   ///////////////////////////////
1368   Bool_t echec   = kFALSE;
1369   Double_t shift = 0.0;
1370   //Calculate the shift in x coresponding to this tnp
1371   if(fNgroupprf != 0.0){
1372     shift      = -3.0*(fNgroupprf-1)-1.5;
1373     Double_t limithigh  = -0.2*(fNgroupprf-1);
1374     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1375     else{
1376       while(tnp > limithigh){
1377         limithigh += 0.2;
1378         shift += 3.0;
1379       }
1380     }
1381   }
1382   // do nothing if out of tnp range
1383   //printf("echec %d\n",(Int_t)echec);
1384   if(echec) return kFALSE;
1385
1386   ///////////////////////
1387   // Variables
1388   //////////////////////
1389
1390   Int_t nb3pc    = 0;              // number of three pads clusters used for fit
1391   Double_t *padPositions;          // to see the difference between the fit and the 3 pad clusters position
1392   padPositions = new Double_t[AliTRDseed::knTimebins];
1393   for(Int_t k = 0; k < AliTRDseed::knTimebins; k++){
1394     padPositions[k] = 0.0;
1395   }
1396   fLinearFitterTracklet->ClearPoints();
1397
1398   //printf("loop clusters \n");
1399   ////////////////////////////
1400   // loop over the clusters
1401   ////////////////////////////
1402   AliTRDcluster *cl                   = 0x0;
1403   for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1404     if(!(cl = tracklet->GetClusters(ic))) continue;
1405
1406     Double_t     time  = cl->GetLocalTimeBin();
1407     if((time<=7) || (time>=21)) continue;
1408     Short_t  *signals  = cl->GetSignals();
1409     Float_t xcenter    = 0.0;
1410     Bool_t  echec1      = kTRUE;
1411
1412     /////////////////////////////////////////////////////////////
1413     // Center 3 balanced: position with the center of the pad
1414     /////////////////////////////////////////////////////////////
1415     if ((((Float_t) signals[3]) > 0.0) &&
1416         (((Float_t) signals[2]) > 0.0) &&
1417         (((Float_t) signals[4]) > 0.0)) {
1418       echec1 = kFALSE;
1419       // Security if the denomiateur is 0
1420       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1421            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1422         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1423           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1424                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1425       }
1426       else {
1427         echec1 = kTRUE;
1428       }
1429     }
1430     if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1431     if(echec1) continue;
1432
1433     ////////////////////////////////////////////////////////
1434     //if no echec1: calculate with the position of the pad
1435     // Position of the cluster
1436     // fill the linear fitter
1437     ///////////////////////////////////////////////////////
1438     Double_t       padPosition = xcenter +  cl->GetPadCol();
1439     padPositions[ic]            = padPosition;
1440     nb3pc++;
1441     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1442
1443
1444   }//clusters loop
1445
1446   //printf("Fin loop clusters \n");
1447   //////////////////////////////
1448   // fit with a straight line
1449   /////////////////////////////
1450   if(nb3pc < 3){
1451     delete [] padPositions;
1452     fLinearFitterTracklet->ClearPoints();
1453     delete [] padPositions;
1454     return kFALSE;
1455   }
1456   fLinearFitterTracklet->Eval();
1457   TVectorD line(2);
1458   fLinearFitterTracklet->GetParameters(line);
1459   Float_t  pointError  = -1.0;
1460   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1461   pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1462   }
1463   fLinearFitterTracklet->ClearPoints();
1464
1465   //printf("PRF second loop \n");
1466   ////////////////////////////////////////////////
1467   // Fill the PRF: Second loop over clusters
1468   //////////////////////////////////////////////
1469   for(int ic=0; ic<AliTRDseed::knTimebins; ic++){
1470     if(!(cl = tracklet->GetClusters(ic))) continue;
1471
1472     Short_t  *signals      = cl->GetSignals();              // signal
1473     Double_t     time      = cl->GetLocalTimeBin();         // time bin
1474     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1475     Float_t padPos         = cl->GetPadCol();               // middle pad
1476     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1477     Float_t ycenter        = 0.0;                           // relative center charge
1478     Float_t ymin           = 0.0;                           // relative left charge
1479     Float_t ymax           = 0.0;                           // relative right charge
1480
1481     ////////////////////////////////////////////////////////////////
1482     // Calculate ycenter, ymin and ymax even for two pad clusters
1483     ////////////////////////////////////////////////////////////////
1484     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1485        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1486       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1487       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1488       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1489       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1490     }
1491
1492     /////////////////////////
1493     // Calibration group
1494     ////////////////////////
1495     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1496     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1497     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1498     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1499     Float_t   xcl          = cl->GetY();                                // y cluster
1500     Float_t   qcl          = cl->GetQ();                                // charge cluster
1501     Int_t     layer        = GetLayer(fDetectorPreviousTrack);          // layer
1502     Int_t     stack        = GetStack(fDetectorPreviousTrack);          // stack
1503     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1504     Double_t  x            = dpad;                                      // reconstructed position moved
1505     Float_t   ep           = pointError;                                // error of fit
1506     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1507     Float_t   signal3      = (Float_t)signals[3];                       // signal
1508     Float_t   signal2      = (Float_t)signals[2];                       // signal
1509     Float_t   signal4      = (Float_t)signals[4];                       // signal
1510     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1511
1512
1513
1514     /////////////////////
1515     // Debug stuff
1516     ////////////////////
1517
1518     if(fDebugLevel > 0){
1519       if ( !fDebugStreamer ) {
1520         //debug stream
1521         TDirectory *backup = gDirectory;
1522         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1523         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1524       }
1525
1526       x = xdiff;
1527       Int_t type=0;
1528       Float_t y = ycenter;
1529       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1530         "caligroup="<<caligroup<<
1531         "detector="<<fDetectorPreviousTrack<<
1532         "layer="<<layer<<
1533         "stack="<<stack<<
1534         "npoints="<<nbclusters<<
1535         "Np="<<nb3pc<<
1536         "ep="<<ep<<
1537         "type="<<type<<
1538         "snp="<<snp<<
1539         "tnp="<<tnp<<
1540         "tgl="<<tgl<<
1541         "dzdx="<<dzdx<<
1542         "padPos="<<padPos<<
1543         "padPosition="<<padPositions[ic]<<
1544         "padPosTracklet="<<padPosTracklet<<
1545         "x="<<x<<
1546         "y="<<y<<
1547         "xcl="<<xcl<<
1548         "qcl="<<qcl<<
1549         "signal1="<<signal1<<
1550         "signal2="<<signal2<<
1551         "signal3="<<signal3<<
1552         "signal4="<<signal4<<
1553         "signal5="<<signal5<<
1554         "time="<<time<<
1555         "\n";
1556       x=-(xdiff+1);
1557       y = ymin;
1558       type=-1;
1559       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1560         "caligroup="<<caligroup<<
1561         "detector="<<fDetectorPreviousTrack<<
1562         "layer="<<layer<<
1563         "stack="<<stack<<
1564         "npoints="<<nbclusters<<
1565         "Np="<<nb3pc<<
1566         "ep="<<ep<<
1567         "type="<<type<<
1568         "snp="<<snp<<
1569         "tnp="<<tnp<<
1570         "tgl="<<tgl<<
1571         "dzdx="<<dzdx<<
1572         "padPos="<<padPos<<
1573         "padPosition="<<padPositions[ic]<<
1574         "padPosTracklet="<<padPosTracklet<<
1575         "x="<<x<<
1576         "y="<<y<<
1577         "xcl="<<xcl<<
1578         "qcl="<<qcl<<
1579         "signal1="<<signal1<<
1580         "signal2="<<signal2<<
1581         "signal3="<<signal3<<
1582         "signal4="<<signal4<<
1583         "signal5="<<signal5<<
1584         "time="<<time<<
1585         "\n";
1586       x=1-xdiff;
1587       y = ymax;
1588       type=1;
1589       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1590         "caligroup="<<caligroup<<
1591         "detector="<<fDetectorPreviousTrack<<
1592         "layer="<<layer<<
1593         "stack="<<stack<<
1594         "npoints="<<nbclusters<<
1595         "Np="<<nb3pc<<
1596         "ep="<<ep<<
1597         "type="<<type<<
1598         "snp="<<snp<<
1599         "tnp="<<tnp<<
1600         "tgl="<<tgl<<
1601         "dzdx="<<dzdx<<
1602         "padPos="<<padPos<<
1603         "padPosition="<<padPositions[ic]<<
1604         "padPosTracklet="<<padPosTracklet<<
1605         "x="<<x<<
1606         "y="<<y<<
1607         "xcl="<<xcl<<
1608         "qcl="<<qcl<<
1609         "signal1="<<signal1<<
1610         "signal2="<<signal2<<
1611         "signal3="<<signal3<<
1612         "signal4="<<signal4<<
1613         "signal5="<<signal5<<
1614         "time="<<time<<
1615         "\n";
1616
1617     }
1618
1619     /////////////////////
1620     // Cuts quality
1621     /////////////////////
1622     if(nbclusters < fNumberClusters) continue;
1623     if(nbclusters > fNumberClustersf) continue;
1624     if(nb3pc <= 5) continue;
1625     if((time >= 21) || (time < 7)) continue;
1626     if(TMath::Abs(qcl) < 80) continue;
1627     if( TMath::Abs(snp) >  1.) continue;
1628
1629
1630     ////////////////////////
1631     // Fill the histos
1632     ///////////////////////
1633     if (fHisto2d) {
1634       if(TMath::Abs(dpad) < 1.5) {
1635         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1636         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1637         //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1638       }
1639       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1640         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1641         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1642       }
1643       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1644         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1645         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1646       }
1647     }
1648     // vector method
1649     if (fVector2d) {
1650       if(TMath::Abs(dpad) < 1.5) {
1651         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1652         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1653       }
1654       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1655         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1656         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1657       }
1658       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1659         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1660         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1661       }
1662     }
1663   } // second loop over clusters
1664
1665
1666   delete [] padPositions;
1667   return kTRUE;
1668
1669 }
1670 ///////////////////////////////////////////////////////////////////////////////////////
1671 // Pad row col stuff: see if masked or not
1672 ///////////////////////////////////////////////////////////////////////////////////////
1673 //_____________________________________________________________________________
1674 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1675 {
1676   //
1677   // See if we are not near a masked pad
1678   //
1679
1680   if(cl->IsMasked()) fGoodTracklet = kFALSE;
1681
1682
1683 }
1684 //_____________________________________________________________________________
1685 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1686 {
1687   //
1688   // See if we are not near a masked pad
1689   //
1690
1691   if (!IsPadOn(detector, col, row)) {
1692     fGoodTracklet = kFALSE;
1693   }
1694
1695   if (col > 0) {
1696     if (!IsPadOn(detector, col-1, row)) {
1697       fGoodTracklet = kFALSE;
1698     }
1699   }
1700
1701   if (col < 143) {
1702     if (!IsPadOn(detector, col+1, row)) {
1703       fGoodTracklet = kFALSE;
1704     }
1705   }
1706
1707 }
1708 //_____________________________________________________________________________
1709 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1710 {
1711   //
1712   // Look in the choosen database if the pad is On.
1713   // If no the track will be "not good"
1714   //
1715
1716   // Get the parameter object
1717   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1718   if (!cal) {
1719     AliInfo("Could not get calibDB");
1720     return kFALSE;
1721   }
1722
1723   if (!cal->IsChamberInstalled(detector)     ||
1724        cal->IsChamberMasked(detector)        ||
1725        cal->IsPadMasked(detector,col,row)) {
1726     return kFALSE;
1727   }
1728   else {
1729     return kTRUE;
1730   }
1731
1732 }
1733 ///////////////////////////////////////////////////////////////////////////////////////
1734 // Calibration groups: calculate the number of groups, localise...
1735 ////////////////////////////////////////////////////////////////////////////////////////
1736 //_____________________________________________________________________________
1737 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1738 {
1739   //
1740   // Calculate the calibration group number for i
1741   //
1742
1743   // Row of the cluster and position in the pad groups
1744   Int_t posr = 0;
1745   if (fCalibraMode->GetNnZ(i) != 0) {
1746     posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1747   }
1748
1749
1750   // Col of the cluster and position in the pad groups
1751   Int_t posc = 0;
1752   if (fCalibraMode->GetNnRphi(i) != 0) {
1753     posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1754   }
1755
1756   return posc*fCalibraMode->GetNfragZ(i)+posr;
1757
1758 }
1759 //____________________________________________________________________________________
1760 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1761 {
1762   //
1763   // Calculate the total number of calibration groups
1764   //
1765
1766   Int_t ntotal = 0;
1767
1768   // All together
1769   if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1770     ntotal = 1;
1771     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1772     return ntotal;
1773   }
1774
1775   // Per Supermodule
1776   if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1777     ntotal = 18;
1778     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1779     return ntotal;
1780   }
1781
1782   // More
1783   fCalibraMode->ModePadCalibration(2,i);
1784   fCalibraMode->ModePadFragmentation(0,2,0,i);
1785   fCalibraMode->SetDetChamb2(i);
1786   ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1787   fCalibraMode->ModePadCalibration(0,i);
1788   fCalibraMode->ModePadFragmentation(0,0,0,i);
1789   fCalibraMode->SetDetChamb0(i);
1790   ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1791   AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1792   return ntotal;
1793
1794 }
1795 //_____________________________________________________________________________
1796 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1797 {
1798   //
1799   // Set the mode of calibration group in the z direction for the parameter i
1800   //
1801
1802   if ((Nz >= 0) &&
1803       (Nz <  5)) {
1804     fCalibraMode->SetNz(i, Nz);
1805   }
1806   else {
1807     AliInfo("You have to choose between 0 and 4");
1808   }
1809
1810 }
1811 //_____________________________________________________________________________
1812 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1813 {
1814   //
1815   // Set the mode of calibration group in the rphi direction for the parameter i
1816   //
1817
1818   if ((Nrphi >= 0) &&
1819       (Nrphi <  7)) {
1820     fCalibraMode->SetNrphi(i ,Nrphi);
1821   }
1822   else {
1823     AliInfo("You have to choose between 0 and 6");
1824   }
1825
1826 }
1827
1828 //_____________________________________________________________________________
1829 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1830 {
1831   //
1832   // Set the mode of calibration group all together
1833   //
1834   if(fVector2d == kTRUE) {
1835     AliInfo("Can not work with the vector method");
1836     return;
1837   }
1838   fCalibraMode->SetAllTogether(i);
1839
1840 }
1841
1842 //_____________________________________________________________________________
1843 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1844 {
1845   //
1846   // Set the mode of calibration group per supermodule
1847   //
1848   if(fVector2d == kTRUE) {
1849     AliInfo("Can not work with the vector method");
1850     return;
1851   }
1852   fCalibraMode->SetPerSuperModule(i);
1853
1854 }
1855
1856 //____________Set the pad calibration variables for the detector_______________
1857 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1858 {
1859   //
1860   // For the detector calcul the first Xbins and set the number of row
1861   // and col pads per calibration groups, the number of calibration
1862   // groups in the detector.
1863   //
1864
1865   // first Xbins of the detector
1866   if (fCH2dOn) {
1867     fCalibraMode->CalculXBins(detector,0);
1868   }
1869   if (fPH2dOn) {
1870     fCalibraMode->CalculXBins(detector,1);
1871   }
1872   if (fPRF2dOn) {
1873     fCalibraMode->CalculXBins(detector,2);
1874   }
1875
1876   // fragmentation of idect
1877   for (Int_t i = 0; i < 3; i++) {
1878     fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1879     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1880                        , (Int_t) GetStack(detector)
1881                        , (Int_t) GetSector(detector),i);
1882   }
1883
1884   return kTRUE;
1885
1886 }
1887 //_____________________________________________________________________________
1888 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1889 {
1890   //
1891   // Should be between 0 and 6
1892   //
1893
1894   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1895     AliInfo("The number of groups must be between 0 and 6!");
1896   }
1897   else {
1898     fNgroupprf = numberGroupsPRF;
1899   }
1900
1901 }
1902 ///////////////////////////////////////////////////////////////////////////////////////////
1903 // Per tracklet: store or reset the info, fill the histos with the info
1904 //////////////////////////////////////////////////////////////////////////////////////////
1905 //_____________________________________________________________________________
1906 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1907 {
1908   //
1909   // Store the infos in fAmpTotal, fPHPlace and fPHValue
1910   // Correct from the gain correction before
1911   //
1912
1913   // time bin of the cluster not corrected
1914   Int_t    time     = cl->GetPadTime();
1915
1916   //Correct for the gain coefficient used in the database for reconstruction
1917   Float_t correctthegain = 1.0;
1918   if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1919   else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1920   Float_t correction    = 1.0;
1921   Float_t normalisation = 6.67;
1922   // we divide with gain in AliTRDclusterizer::Transform...
1923   if( correctthegain > 0 ) normalisation /= correctthegain;
1924
1925
1926   // take dd/dl corrected from the angle of the track
1927   correction = dqdl / (normalisation);
1928
1929
1930   // Fill the fAmpTotal with the charge
1931   if (fCH2dOn) {
1932     if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1933   }
1934
1935   // Fill the fPHPlace and value
1936   if (fPH2dOn) {
1937     fPHPlace[time] = group[1];
1938     fPHValue[time] = correction;
1939   }
1940
1941 }
1942 //____________Offine tracking in the AliTRDtracker_____________________________
1943 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1944 {
1945   //
1946   // Reset values per tracklet
1947   //
1948
1949   //Reset good tracklet
1950   fGoodTracklet = kTRUE;
1951
1952   // Reset the fPHValue
1953   if (fPH2dOn) {
1954     //Reset the fPHValue and fPHPlace
1955     for (Int_t k = 0; k < fTimeMax; k++) {
1956       fPHValue[k] = 0.0;
1957       fPHPlace[k] = -1;
1958     }
1959   }
1960
1961   // Reset the fAmpTotal where we put value
1962   if (fCH2dOn) {
1963     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1964       fAmpTotal[k] = 0.0;
1965     }
1966   }
1967 }
1968 //____________Offine tracking in the AliTRDtracker_____________________________
1969 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1970 {
1971   //
1972   // For the offline tracking or mcm tracklets
1973   // This function will be called in the functions UpdateHistogram...
1974   // to fill the info of a track for the relativ gain calibration
1975   //
1976
1977   Int_t nb            =  0;   // Nombre de zones traversees
1978   Int_t fd            = -1;   // Premiere zone non nulle
1979   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1980
1981   if(nbclusters < fNumberClusters) return;
1982   if(nbclusters > fNumberClustersf) return;
1983
1984
1985   // Normalize with the number of clusters
1986   Double_t normalizeCst = fRelativeScale;
1987   if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1988
1989   // See if the track goes through different zones
1990   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1991     if (fAmpTotal[k] > 0.0) {
1992       totalcharge += fAmpTotal[k];
1993       nb++;
1994       if (nb == 1) {
1995         fd = k;
1996       }
1997     }
1998   }
1999
2000
2001   switch (nb)
2002     {
2003     case 1:
2004       fNumberUsedCh[0]++;
2005       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2006       if (fHisto2d) {
2007         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2008         //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2009       }
2010       if (fVector2d) {
2011         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2012       }
2013       break;
2014     case 2:
2015       if ((fAmpTotal[fd]   > 0.0) &&
2016           (fAmpTotal[fd+1] > 0.0)) {
2017         // One of the two very big
2018         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2019           if (fHisto2d) {
2020             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2021             //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2022           }
2023           if (fVector2d) {
2024             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2025           }
2026           fNumberUsedCh[1]++;
2027           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2028         }
2029         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2030           if (fHisto2d) {
2031             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2032             //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2033           }
2034           if (fVector2d) {
2035             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2036           }
2037           fNumberUsedCh[1]++;
2038           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2039         }
2040       }
2041       if (fCalibraMode->GetNfragZ(0) > 1) {
2042         if (fAmpTotal[fd] > 0.0) {
2043           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2044             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2045               // One of the two very big
2046               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2047                 if (fHisto2d) {
2048                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2049                   //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2050                 }
2051                 if (fVector2d) {
2052                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2053                 }
2054                 fNumberUsedCh[1]++;
2055                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2056               }
2057               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2058                 if (fHisto2d) {
2059                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2060                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2061                 }
2062                 fNumberUsedCh[1]++;
2063                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2064                 if (fVector2d) {
2065                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2066                 }
2067               }
2068             }
2069           }
2070         }
2071       }
2072       break;
2073     default: break;
2074     }
2075 }
2076 //____________Offine tracking in the AliTRDtracker_____________________________
2077 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2078 {
2079   //
2080   // For the offline tracking or mcm tracklets
2081   // This function will be called in the functions UpdateHistogram...
2082   // to fill the info of a track for the drift velocity  calibration
2083   //
2084
2085   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
2086   Int_t fd1 = -1; // Premiere zone non nulle
2087   Int_t fd2 = -1; // Deuxieme zone non nulle
2088   Int_t k1  = -1; // Debut de la premiere zone
2089   Int_t k2  = -1; // Debut de la seconde zone
2090   Int_t nbclusters = 0; // number of clusters
2091
2092
2093
2094   // See if the track goes through different zones
2095   for (Int_t k = 0; k < fTimeMax; k++) {
2096     if (fPHValue[k] > 0.0) {
2097       nbclusters++;
2098       if (fd1 == -1) {
2099         fd1 = fPHPlace[k];
2100         k1  = k;
2101       }
2102       if (fPHPlace[k] != fd1) {
2103         if (fd2 == -1) {
2104           k2  = k;
2105           fd2 = fPHPlace[k];
2106           nb  = 2;
2107         }
2108         if (fPHPlace[k] != fd2) {
2109           nb = 3;
2110         }
2111       }
2112     }
2113   }
2114
2115   // See if noise before and after
2116   if(fMaxCluster > 0) {
2117     if(fPHValue[0] > fMaxCluster) return;
2118     if(fTimeMax > fNbMaxCluster) {
2119       for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2120         if(fPHValue[k] > fMaxCluster) return;
2121       }
2122     }
2123   }
2124
2125   //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2126
2127   if(nbclusters < fNumberClusters) return;
2128   if(nbclusters > fNumberClustersf) return;
2129
2130   switch(nb)
2131     {
2132     case 1:
2133       fNumberUsedPh[0]++;
2134       for (Int_t i = 0; i < fTimeMax; i++) {
2135         if (fHisto2d) {
2136           if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2137           else {
2138             if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2139               }
2140           //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2141         }
2142         if (fVector2d) {
2143           if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2144           else {
2145             if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2146           }
2147         }
2148       }
2149       break;
2150     case 2:
2151       if ((fd1 == fd2+1) ||
2152           (fd2 == fd1+1)) {
2153         // One of the two fast all the think
2154         if (k2 > (k1+fDifference)) {
2155           //we choose to fill the fd1 with all the values
2156           fNumberUsedPh[1]++;
2157           for (Int_t i = 0; i < fTimeMax; i++) {
2158             if (fHisto2d) {
2159               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2160               else {
2161                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2162                   }
2163             }
2164             if (fVector2d) {
2165               if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2166               else {
2167                 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2168                   }
2169             }
2170           }
2171         }
2172         if ((k2+fDifference) < fTimeMax) {
2173           //we choose to fill the fd2 with all the values
2174           fNumberUsedPh[1]++;
2175           for (Int_t i = 0; i < fTimeMax; i++) {
2176             if (fHisto2d) {
2177               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2178               else {
2179                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2180               }
2181             }
2182           if (fVector2d) {
2183             if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2184             else {
2185               if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2186             }
2187           }
2188           }
2189         }
2190       }
2191       // Two zones voisines sinon rien!
2192       if (fCalibraMode->GetNfragZ(1) > 1) {
2193         // Case 2
2194         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2195           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2196             // One of the two fast all the think
2197             if (k2 > (k1+fDifference)) {
2198               //we choose to fill the fd1 with all the values
2199               fNumberUsedPh[1]++;
2200               for (Int_t i = 0; i < fTimeMax; i++) {
2201                 if (fHisto2d) {
2202                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2203                   else {
2204                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2205                   }
2206                 }
2207                 if (fVector2d) {
2208                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2209                   else {
2210                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2211                   }
2212                 }
2213               }
2214             }
2215             if ((k2+fDifference) < fTimeMax) {
2216               //we choose to fill the fd2 with all the values
2217               fNumberUsedPh[1]++;
2218               for (Int_t i = 0; i < fTimeMax; i++) {
2219                 if (fHisto2d) {
2220                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2221                   else {
2222                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2223                   }
2224                 }
2225                 if (fVector2d) {
2226                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2227                   else {
2228                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2229                   }
2230                 }
2231               }
2232             }
2233           }
2234         }
2235         // Two zones voisines sinon rien!
2236         // Case 3
2237         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2238           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2239             // One of the two fast all the think
2240             if (k2 > (k1 + fDifference)) {
2241               //we choose to fill the fd1 with all the values
2242               fNumberUsedPh[1]++;
2243               for (Int_t i = 0; i < fTimeMax; i++) {
2244                 if (fHisto2d) {
2245                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2246                   else {
2247                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2248                   }
2249                 }
2250                 if (fVector2d) {
2251                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2252                   else {
2253                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2254                   }
2255                 }
2256               }
2257             }
2258             if ((k2+fDifference) < fTimeMax) {
2259               //we choose to fill the fd2 with all the values
2260               fNumberUsedPh[1]++;
2261               for (Int_t i = 0; i < fTimeMax; i++) {
2262                 if (fHisto2d) {
2263                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2264                   else {
2265                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2266                   }
2267                 }
2268                 if (fVector2d) {
2269                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2270                   else {
2271                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2272                   }
2273                 }
2274               }
2275             }
2276           }
2277         }
2278       }
2279       break;
2280     default: break;
2281     }
2282 }
2283 //////////////////////////////////////////////////////////////////////////////////////////
2284 // DAQ process functions
2285 /////////////////////////////////////////////////////////////////////////////////////////
2286 //_____________________________________________________________________
2287 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2288 {
2289   //
2290   // Event Processing loop - AliTRDrawStreamBase
2291   // TestBeam 2007 version
2292   // 0 timebin problem
2293   // 1 no input
2294   // 2 input
2295   //
2296   // Same algorithm as TestBeam but different reader
2297   //
2298
2299   Int_t withInput = 1;
2300
2301   Double_t phvalue[16][144][36];
2302   for(Int_t k = 0; k < 36; k++){
2303     for(Int_t j = 0; j < 16; j++){
2304       for(Int_t c = 0; c < 144; c++){
2305         phvalue[j][c][k] = 0.0;
2306       }
2307     }
2308   }
2309
2310   fDetectorPreviousTrack = -1;
2311   fMCMPrevious           = -1;
2312   fROBPrevious           = -1;
2313   Int_t nbtimebin = 0;
2314   Int_t baseline  = 10;
2315
2316
2317   if(!nocheck){
2318
2319     fTimeMax = 0;
2320
2321     while (rawStream->Next()) {
2322
2323       Int_t idetector = rawStream->GetDet();                            //  current detector
2324       Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2325       Int_t irob      = rawStream->GetROB();                            //  current ROB
2326
2327       //printf("Detector %d\n",idetector);
2328
2329       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2330
2331         // Fill
2332         withInput = TMath::Max(FillDAQ(phvalue),withInput);
2333
2334
2335         // reset
2336         for(Int_t k = 0; k < 36; k++){
2337           for(Int_t j = 0; j < 16; j++){
2338             for(Int_t c = 0; c < 144; c++){
2339               phvalue[j][c][k] = 0.0;
2340             }
2341           }
2342         }
2343       }
2344
2345       fDetectorPreviousTrack = idetector;
2346       fMCMPrevious           = imcm;
2347       fROBPrevious           = irob;
2348
2349       nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2350       if(nbtimebin == 0) return 0;
2351       if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2352       fTimeMax          = nbtimebin;
2353
2354       //baseline          = rawStream->GetCommonAdditive();                // common additive baseline
2355       fNumberClustersf    = fTimeMax;
2356       fNumberClusters     = (Int_t)(0.6*fTimeMax);
2357
2358
2359       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
2360       Int_t  col        = rawStream->GetCol();
2361       Int_t row         = rawStream->GetRow();
2362
2363
2364       //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2365
2366
2367       for(Int_t itime = 0; itime < nbtimebin; itime++){
2368         phvalue[row][col][itime] = signal[itime]-baseline;
2369       }
2370     }
2371
2372     // fill the last one
2373     if(fDetectorPreviousTrack != -1){
2374
2375       // Fill
2376       withInput = TMath::Max(FillDAQ(phvalue),withInput);
2377
2378       // reset
2379       for(Int_t k = 0; k < 36; k++){
2380         for(Int_t j = 0; j < 16; j++){
2381           for(Int_t c = 0; c < 144; c++){
2382             phvalue[j][c][k] = 0.0;
2383           }
2384         }
2385       }
2386     }
2387
2388   }
2389   else{
2390
2391     while (rawStream->Next()) {
2392
2393       Int_t idetector = rawStream->GetDet();                            //  current detector
2394       Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2395       Int_t irob      = rawStream->GetROB();                            //  current ROB
2396
2397       //printf("Detector %d\n",idetector);
2398
2399       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2400
2401         // Fill
2402         withInput = TMath::Max(FillDAQ(phvalue),withInput);
2403
2404         // reset
2405         for(Int_t k = 0; k < 36; k++){
2406           for(Int_t j = 0; j < 16; j++){
2407             for(Int_t c = 0; c < 144; c++){
2408               phvalue[j][c][k] = 0.0;
2409             }
2410           }
2411         }
2412       }
2413
2414       fDetectorPreviousTrack = idetector;
2415       fMCMPrevious           = imcm;
2416       fROBPrevious           = irob;
2417
2418       //baseline          = rawStream->GetCommonAdditive();                //  common baseline
2419
2420       fTimeMax          = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2421       fNumberClustersf    = fTimeMax;
2422       fNumberClusters     = (Int_t)(0.6*fTimeMax);
2423       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
2424       Int_t col         = rawStream->GetCol();
2425       Int_t row         = rawStream->GetRow();
2426
2427
2428       //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2429
2430       for(Int_t itime = 0; itime < fTimeMax; itime++){
2431         phvalue[row][col][itime] = signal[itime]-baseline;
2432       }
2433     }
2434
2435     // fill the last one
2436     if(fDetectorPreviousTrack != -1){
2437
2438       // Fill
2439       withInput = TMath::Max(FillDAQ(phvalue),withInput);
2440
2441       // reset
2442       for(Int_t k = 0; k < 36; k++){
2443         for(Int_t j = 0; j < 16; j++){
2444           for(Int_t c = 0; c < 144; c++){
2445             phvalue[j][c][k] = 0.0;
2446           }
2447         }
2448       }
2449     }
2450   }
2451
2452   return withInput;
2453
2454 }
2455 //_____________________________________________________________________
2456 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2457 {
2458   //
2459   //  Event processing loop - AliRawReader
2460   //  Testbeam 2007 version
2461   //
2462
2463   AliTRDrawStreamBase rawStream(rawReader);
2464
2465   rawReader->Select("TRD");
2466
2467   return ProcessEventDAQ(&rawStream, nocheck);
2468 }
2469
2470 //_________________________________________________________________________
2471 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2472 #ifdef ALI_DATE
2473                                                eventHeaderStruct *event,
2474                                                Bool_t nocheck
2475 #else
2476                                                eventHeaderStruct* /*event*/,
2477                                                Bool_t /*nocheck*/
2478
2479 #endif
2480                                    )
2481 {
2482   //
2483   //  process date event
2484   //  Testbeam 2007 version
2485   //
2486 #ifdef ALI_DATE
2487     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2488     Int_t result=ProcessEventDAQ(rawReader, nocheck);
2489     delete rawReader;
2490     return result;
2491 #else
2492     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2493     return 0;
2494 #endif
2495
2496 }
2497 //////////////////////////////////////////////////////////////////////////////
2498 // Routine inside the DAQ process
2499 /////////////////////////////////////////////////////////////////////////////
2500 //_______________________________________________________________________
2501 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2502
2503   //
2504   // Look for the maximum by collapsing over the time
2505   // Sum over four pad col and two pad row
2506   //
2507
2508   Int_t used = 0;
2509
2510
2511   Int_t idect = fDetectorPreviousTrack;
2512   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2513   Double_t sum[36];
2514   for(Int_t tb = 0; tb < 36; tb++){
2515     sum[tb] = 0.0;
2516   }
2517
2518   //fGoodTracklet = kTRUE;
2519   //fDetectorPreviousTrack = 0;
2520
2521
2522   ///////////////////////////
2523   // look for maximum
2524   /////////////////////////
2525
2526   Int_t imaxRow = 0;
2527   Int_t imaxCol = 0;
2528   Double_t integralMax = -1;
2529
2530   for (Int_t ir = 1; ir <= 15; ir++)
2531     {
2532       for (Int_t ic = 2; ic <= 142; ic++)
2533         {
2534           Double_t integral = 0;
2535           for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2536             {
2537               for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2538                 {
2539                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2540                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2541                     {
2542
2543                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2544                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2545                       }// addtb
2546                     } //addsignal
2547                 } //shiftC
2548             } // shiftR
2549
2550           if (integralMax < integral)
2551             {
2552               imaxRow = ir;
2553               imaxCol = ic;
2554               integralMax = integral;
2555             } // check max integral
2556         } //ic
2557     } // ir
2558
2559   //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2560
2561   if((imaxRow == 0) || (imaxCol == 0)) {
2562     used=1;
2563     return used;
2564   }
2565   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2566   //if(!fGoodTracklet) used = 1;;
2567
2568   //  /////////////////////////////////////////////////////
2569   // sum ober 2 row and 4 pad cols for each time bins
2570   //  ////////////////////////////////////////////////////
2571
2572
2573   for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2574     {
2575       for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2576         {
2577           for(Int_t it = 0; it < fTimeMax; it++){
2578             sum[it] += phvalue[ir][ic][it];
2579           }
2580         }//ic
2581     }//ir
2582
2583   Int_t nbcl = 0;
2584   Double_t sumcharge = 0.0;
2585   for(Int_t it = 0; it < fTimeMax; it++){
2586     sumcharge += sum[it];
2587     if(sum[it] > 20.0) nbcl++;
2588   }
2589
2590
2591   /////////////////////////////////////////////////////////
2592   // Debug
2593   ////////////////////////////////////////////////////////
2594   if(fDebugLevel > 0){
2595     if ( !fDebugStreamer ) {
2596       //debug stream
2597       TDirectory *backup = gDirectory;
2598       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2599       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2600     }
2601
2602     Double_t amph0 = sum[0];
2603     Double_t amphlast = sum[fTimeMax-1];
2604     Double_t rms      = TMath::RMS(fTimeMax,sum);
2605     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2606     for(Int_t it = 0; it < fTimeMax; it++){
2607       Double_t clustera = sum[it];
2608
2609     (* fDebugStreamer) << "FillDAQa"<<
2610       "ampTotal="<<sumcharge<<
2611       "row="<<imaxRow<<
2612       "col="<<imaxCol<<
2613       "detector="<<idect<<
2614       "amph0="<<amph0<<
2615       "amphlast="<<amphlast<<
2616       "goodtracklet="<<goodtracklet<<
2617       "clustera="<<clustera<<
2618       "it="<<it<<
2619       "rms="<<rms<<
2620       "\n";
2621     }
2622   }
2623
2624   ////////////////////////////////////////////////////////
2625   // fill
2626   ///////////////////////////////////////////////////////
2627   if(sum[0] > 100.0) used = 1;
2628   if(nbcl < fNumberClusters) used = 1;
2629   if(nbcl > fNumberClustersf) used = 1;
2630
2631   //if(fDetectorPreviousTrack == 15){
2632   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2633   //}
2634   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2635   if(used == 0){
2636     for(Int_t it = 0; it < fTimeMax; it++){
2637       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2638       else{
2639         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2640       }
2641     }
2642
2643
2644     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2645     used = 2;
2646     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2647
2648   }
2649
2650   return used;
2651
2652 }
2653 //____________Online trackling in AliTRDtrigger________________________________
2654 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2655 {
2656   //
2657   // For the DAQ
2658   // Fill a simple average pulse height
2659   //
2660
2661
2662   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2663
2664
2665   return kTRUE;
2666
2667 }
2668 //____________Write_____________________________________________________
2669 //_____________________________________________________________________
2670 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2671 {
2672   //
2673   //  Write infos to file
2674   //
2675
2676   //For debugging
2677   if ( fDebugStreamer ) {
2678     delete fDebugStreamer;
2679     fDebugStreamer = 0x0;
2680   }
2681
2682   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2683                ,fNumberTrack
2684                ,fNumberUsedCh[0]
2685                ,fNumberUsedCh[1]
2686                ,fNumberUsedPh[0]
2687                ,fNumberUsedPh[1]));
2688
2689   TDirectory *backup = gDirectory;
2690   TString option;
2691
2692   if ( append )
2693     option = "update";
2694   else
2695     option = "recreate";
2696
2697   TFile f(filename,option.Data());
2698
2699   TStopwatch stopwatch;
2700   stopwatch.Start();
2701   if(fVector2d) {
2702     f.WriteTObject(fCalibraVector);
2703   }
2704
2705   if (fCH2dOn ) {
2706     if (fHisto2d) {
2707       f.WriteTObject(fCH2d);
2708     }
2709   }
2710   if (fPH2dOn ) {
2711     if (fHisto2d) {
2712       f.WriteTObject(fPH2d);
2713     }
2714   }
2715   if (fPRF2dOn) {
2716     if (fHisto2d) {
2717         f.WriteTObject(fPRF2d);
2718     }
2719   }
2720   if(fLinearFitterOn){
2721     AnalyseLinearFitter();
2722     f.WriteTObject(fLinearVdriftFit);
2723   }
2724
2725   f.Close();
2726
2727   if ( backup ) backup->cd();
2728
2729   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2730                ,stopwatch.RealTime(),stopwatch.CpuTime()));
2731 }
2732 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2733 // Stats stuff
2734 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2735 //___________________________________________probe the histos__________________________________________________
2736 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2737 {
2738   //
2739   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2740   // debug mode with 2 for TH2I and 3 for TProfile2D
2741   // It gives a pointer to a Double_t[7] with the info following...
2742   // [0] : number of calibration groups with entries
2743   // [1] : minimal number of entries found
2744   // [2] : calibration group number of the min
2745   // [3] : maximal number of entries found
2746   // [4] : calibration group number of the max
2747   // [5] : mean number of entries found
2748   // [6] : mean relative error
2749   //
2750
2751   Double_t *info = new Double_t[7];
2752
2753   // Number of Xbins (detectors or groups of pads)
2754   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
2755   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
2756
2757   // Initialise
2758   Double_t nbwe = 0; //number of calibration groups with entries
2759   Double_t minentries = 0; //minimal number of entries found
2760   Double_t maxentries = 0; //maximal number of entries found
2761   Double_t placemin = 0; //calibration group number of the min
2762   Double_t placemax = -1; //calibration group number of the max
2763   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2764   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2765
2766   Double_t counter = 0;
2767
2768   //Debug
2769   TH1F *nbEntries = 0x0;//distribution of the number of entries
2770   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2771   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2772
2773   // Beginning of the loop over the calibration groups
2774   for (Int_t idect = 0; idect < nbins; idect++) {
2775
2776     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2777     projch->SetDirectory(0);
2778
2779     // Number of entries for this calibration group
2780     Double_t nentries = 0.0;
2781     if((i%2) == 0){
2782       for (Int_t k = 0; k < nxbins; k++) {
2783         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2784       }
2785     }
2786     else{
2787       for (Int_t k = 0; k < nxbins; k++) {
2788         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2789         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2790           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2791           counter++;
2792         }
2793       }
2794     }
2795
2796     //Debug
2797     if(i > 1){
2798       if((!((Bool_t)nbEntries)) && (nentries > 0)){
2799         nbEntries = new TH1F("Number of entries","Number of entries"
2800                                ,100,(Int_t)nentries/2,nentries*2);
2801         nbEntries->SetDirectory(0);
2802         nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
2803                                ,nbins,0,nbins);
2804         nbEntriesPerGroup->SetDirectory(0);
2805         nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
2806                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2807         nbEntriesPerSp->SetDirectory(0);
2808       }
2809       if(nbEntries){
2810         if(nentries > 0) nbEntries->Fill(nentries);
2811         nbEntriesPerGroup->Fill(idect+0.5,nentries);
2812         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2813       }
2814     }
2815
2816     //min amd max
2817     if(nentries > maxentries){
2818       maxentries = nentries;
2819       placemax = idect;
2820     }
2821     if(idect == 0) {
2822       minentries = nentries;
2823     }
2824     if(nentries < minentries){
2825       minentries = nentries;
2826       placemin = idect;
2827     }
2828     //nbwe
2829     if(nentries > 0) {
2830       nbwe++;
2831       meanstats += nentries;
2832     }
2833   }//calibration groups loop
2834
2835   if(nbwe > 0) meanstats /= nbwe;
2836   if(counter > 0) meanrelativerror /= counter;
2837
2838   AliInfo(Form("There are %f calibration groups with entries",nbwe));
2839   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2840   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2841   AliInfo(Form("The mean number of entries is %f",meanstats));
2842   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2843
2844   info[0] = nbwe;
2845   info[1] = minentries;
2846   info[2] = placemin;
2847   info[3] = maxentries;
2848   info[4] = placemax;
2849   info[5] = meanstats;
2850   info[6] = meanrelativerror;
2851
2852   if(i > 1){
2853     gStyle->SetPalette(1);
2854     gStyle->SetOptStat(1111);
2855     gStyle->SetPadBorderMode(0);
2856     gStyle->SetCanvasColor(10);
2857     gStyle->SetPadLeftMargin(0.13);
2858     gStyle->SetPadRightMargin(0.01);
2859     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2860     stat->Divide(2,1);
2861     stat->cd(1);
2862     nbEntries->Draw("");
2863     stat->cd(2);
2864     nbEntriesPerSp->SetStats(0);
2865     nbEntriesPerSp->Draw("");
2866     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2867     stat1->cd();
2868     nbEntriesPerGroup->SetStats(0);
2869     nbEntriesPerGroup->Draw("");
2870   }
2871
2872   return info;
2873
2874 }
2875 //____________________________________________________________________________
2876 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2877 {
2878   //
2879   // Return a Int_t[4] with:
2880   // 0 Mean number of entries
2881   // 1 median of number of entries
2882   // 2 rms of number of entries
2883   // 3 number of group with entries
2884   //
2885
2886   Double_t *stat      = new Double_t[4];
2887   stat[3]             = 0.0;
2888
2889   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
2890   Double_t *weight    = new Double_t[nbofgroups];
2891   Int_t    *nonul     = new Int_t[nbofgroups];
2892
2893   for(Int_t k = 0; k < nbofgroups; k++){
2894     if(fEntriesCH[k] > 0) {
2895       weight[k] = 1.0;
2896       nonul[(Int_t)stat[3]] = fEntriesCH[k];
2897       stat[3]++;
2898     }
2899     else weight[k] = 0.0;
2900   }
2901   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
2902   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
2903   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2904
2905   return stat;
2906
2907 }
2908 //____________________________________________________________________________
2909 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2910 {
2911   //
2912   // Return a Int_t[4] with:
2913   // 0 Mean number of entries
2914   // 1 median of number of entries
2915   // 2 rms of number of entries
2916   // 3 number of group with entries
2917   //
2918
2919   Double_t *stat      = new Double_t[4];
2920   stat[3]             = 0.0;
2921
2922   Int_t    nbofgroups = 540;
2923   Double_t *weight    = new Double_t[nbofgroups];
2924   Int_t    *nonul     = new Int_t[nbofgroups];
2925
2926   for(Int_t k = 0; k < nbofgroups; k++){
2927     if(fEntriesLinearFitter[k] > 0) {
2928       weight[k] = 1.0;
2929       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2930       stat[3]++;
2931     }
2932     else weight[k] = 0.0;
2933   }
2934   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2935   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2936   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2937
2938   return stat;
2939
2940 }
2941 //////////////////////////////////////////////////////////////////////////////////////
2942 // Create Histos
2943 //////////////////////////////////////////////////////////////////////////////////////
2944 //_____________________________________________________________________________
2945 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2946 {
2947   //
2948   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2949   // If fNgroupprf is zero then no binning in tnp
2950   //
2951
2952   TString name("Nz");
2953   name += fCalibraMode->GetNz(2);
2954   name += "Nrphi";
2955   name += fCalibraMode->GetNrphi(2);
2956   name += "Ngp";
2957   name += fNgroupprf;
2958
2959   if(fNgroupprf != 0){
2960
2961     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2962                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2963     fPRF2d->SetYTitle("Det/pad groups");
2964     fPRF2d->SetXTitle("Position x/W [pad width units]");
2965     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2966     fPRF2d->SetStats(0);
2967   }
2968   else{
2969     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2970                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2971     fPRF2d->SetYTitle("Det/pad groups");
2972     fPRF2d->SetXTitle("Position x/W [pad width units]");
2973     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2974     fPRF2d->SetStats(0);
2975   }
2976
2977 }
2978
2979 //_____________________________________________________________________________
2980 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2981 {
2982   //
2983   // Create the 2D histos
2984   //
2985
2986   TString name("Nz");
2987   name += fCalibraMode->GetNz(1);
2988   name += "Nrphi";
2989   name += fCalibraMode->GetNrphi(1);
2990
2991   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2992                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2993                          ,nn,0,nn);
2994   fPH2d->SetYTitle("Det/pad groups");
2995   fPH2d->SetXTitle("time [#mus]");
2996   fPH2d->SetZTitle("<PH> [a.u.]");
2997   fPH2d->SetStats(0);
2998
2999 }
3000 //_____________________________________________________________________________
3001 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3002 {
3003   //
3004   // Create the 2D histos
3005   //
3006
3007   TString name("Nz");
3008   name += fCalibraMode->GetNz(0);
3009   name += "Nrphi";
3010   name += fCalibraMode->GetNrphi(0);
3011
3012   fCH2d = new TH2I("CH2d",(const Char_t *) name
3013                    ,fNumberBinCharge,0,300,nn,0,nn);
3014   fCH2d->SetYTitle("Det/pad groups");
3015   fCH2d->SetXTitle("charge deposit [a.u]");
3016   fCH2d->SetZTitle("counts");
3017   fCH2d->SetStats(0);
3018   fCH2d->Sumw2();
3019
3020 }
3021 //////////////////////////////////////////////////////////////////////////////////
3022 // Set relative scale
3023 /////////////////////////////////////////////////////////////////////////////////
3024 //_____________________________________________________________________________
3025 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3026 {
3027   //
3028   // Set the factor that will divide the deposited charge
3029   // to fit in the histo range [0,300]
3030   //
3031
3032   if (RelativeScale > 0.0) {
3033     fRelativeScale = RelativeScale;
3034   }
3035   else {
3036     AliInfo("RelativeScale must be strict positif!");
3037   }
3038
3039 }
3040 //////////////////////////////////////////////////////////////////////////////////
3041 // Quick way to fill a histo
3042 //////////////////////////////////////////////////////////////////////////////////
3043 //_____________________________________________________________________
3044 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3045 {
3046   //
3047   // FillCH2d: Marian style
3048   //
3049
3050   //skip simply the value out of range
3051   if((y>=300.0) || (y<0.0)) return;
3052
3053   //Calcul the y place
3054   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3055   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3056
3057   //Fill
3058   fCH2d->GetArray()[place]++;
3059
3060 }
3061
3062 //////////////////////////////////////////////////////////////////////////////////
3063 // Geometrical functions
3064 ///////////////////////////////////////////////////////////////////////////////////
3065 //_____________________________________________________________________________
3066 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3067 {
3068   //
3069   // Reconstruct the layer number from the detector number
3070   //
3071
3072   return ((Int_t) (d % 6));
3073
3074 }
3075
3076 //_____________________________________________________________________________
3077 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3078 {
3079   //
3080   // Reconstruct the stack number from the detector number
3081   //
3082   const Int_t kNlayer = 6;
3083
3084   return ((Int_t) (d % 30) / kNlayer);
3085
3086 }
3087
3088 //_____________________________________________________________________________
3089 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3090 {
3091   //
3092   // Reconstruct the sector number from the detector number
3093   //
3094   Int_t fg = 30;
3095
3096   return ((Int_t) (d / fg));
3097
3098 }
3099 ///////////////////////////////////////////////////////////////////////////////////
3100 // Getter functions for DAQ of the CH2d and the PH2d
3101 //////////////////////////////////////////////////////////////////////////////////
3102 //_____________________________________________________________________
3103 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3104 {
3105     //
3106     // return pointer to fPH2d TProfile2D
3107     // create a new TProfile2D if it doesn't exist allready
3108     //
3109     if ( fPH2d )
3110         return fPH2d;
3111
3112     // Some parameters
3113     fTimeMax = nbtimebin;
3114     fSf      = samplefrequency;
3115
3116     CreatePH2d(540);
3117
3118     return fPH2d;
3119 }
3120 //_____________________________________________________________________
3121 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3122 {
3123     //
3124     // return pointer to fCH2d TH2I
3125     // create a new TH2I if it doesn't exist allready
3126     //
3127     if ( fCH2d )
3128         return fCH2d;
3129
3130     CreateCH2d(540);
3131
3132     return fCH2d;
3133 }
3134 ////////////////////////////////////////////////////////////////////////////////////////////
3135 // Drift velocity calibration
3136 ///////////////////////////////////////////////////////////////////////////////////////////
3137 //_____________________________________________________________________
3138 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3139 {
3140     //
3141     // return pointer to TLinearFitter Calibration
3142     // if force is true create a new TLinearFitter if it doesn't exist allready
3143     //
3144
3145   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3146     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3147   }
3148
3149   // if we are forced and TLinearFitter doesn't yet exist create it
3150
3151   // new TLinearFitter
3152   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3153   fLinearFitterArray.AddAt(linearfitter,detector);
3154   return linearfitter;
3155 }
3156
3157 //____________________________________________________________________________
3158 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3159 {
3160   //
3161   // Analyse array of linear fitter because can not be written
3162   // Store two arrays: one with the param the other one with the error param + number of entries
3163   //
3164
3165   for(Int_t k = 0; k < 540; k++){
3166     TLinearFitter *linearfitter = GetLinearFitter(k);
3167     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3168       TVectorD  *par  = new TVectorD(2);
3169       TVectorD   pare = TVectorD(2);
3170       TVectorD  *parE = new TVectorD(3);
3171       linearfitter->Eval();
3172       linearfitter->GetParameters(*par);
3173       linearfitter->GetErrors(pare);
3174       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3175       (*parE)[0] = pare[0]*ppointError;
3176       (*parE)[1] = pare[1]*ppointError;
3177       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3178       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3179       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3180     }
3181   }
3182 }
3183