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