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