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