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