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