TRD/AliTRDCalibraFillHisto.cxx

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