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