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