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<AliTRDseedV1::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<AliTRDseedV1::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     // Check no shared clusters
 947     //////////////////////////////
 948     for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
 949       if((cl = tracklet->GetClusters(icc)))  crossrow = 1;
 950     }
 951   }
 952
 953   ////////////////////////////////////
 954   // Do the straight line fit now
 955   ///////////////////////////////////
 956   if(nbli <= 2){
 957     fLinearFitterTracklet->ClearPoints();
 958     return kFALSE;
 959   }
 960   TVectorD pars;
 961   fLinearFitterTracklet->Eval();
 962   fLinearFitterTracklet->GetParameters(pars);
 963   pointError  =  TMath::Sqrt(fLinearFitterTracklet->GetChisquare()/(nbli-2));
 964   errorpar    =  fLinearFitterTracklet->GetParError(1)*pointError;
 965   dydt        = pars[1];
 966   //printf("chis %f, nbli %d, pointError %f, parError %f, errorpar %f\n",fLinearFitterTracklet->GetChisquare(),nbli,pointError,fLinearFitterTracklet->GetParError(1),errorpar);
 967   fLinearFitterTracklet->ClearPoints();
 968
 969   ////////////////////////////////
 970   // Debug stuff
 971   ///////////////////////////////
 972
 973
 974   if(fDebugLevel > 0){
 975     if ( !fDebugStreamer ) {
 976       //debug stream
 977       TDirectory *backup = gDirectory;
 978       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
 979       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 980     }
 981
 982
 983     Int_t layer = GetLayer(fDetectorPreviousTrack);
 984
 985     (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
 986       //"snpright="<<snpright<<
 987       "nbli="<<nbli<<
 988       "nbclusters="<<nbclusters<<
 989       "detector="<<fDetectorPreviousTrack<<
 990       "layer="<<layer<<
 991       "snp="<<snp<<
 992       "tnp="<<tnp<<
 993       "tgl="<<tgl<<
 994       "tnt="<<tnt<<
 995       "dydt="<<dydt<<
 996       "dzdx="<<dzdx<<
 997       "crossrow="<<crossrow<<
 998       "errorpar="<<errorpar<<
 999       "pointError="<<pointError<<
1000       "\n";
1001
1002   }
1003
1004   /////////////////////////
1005   // Cuts quality
1006   ////////////////////////
1007
1008   if(nbclusters < fNumberClusters) return kFALSE;
1009   if(nbclusters > fNumberClustersf) return kFALSE;
1010   if(pointError >= 0.3) return kFALSE;
1011   if(crossrow == 1) return kFALSE;
1012
1013   ///////////////////////
1014   // Fill
1015   //////////////////////
1016
1017   if(fLinearFitterOn){
1018     //Add to the linear fitter of the detector
1019     if( TMath::Abs(snp) <  1.){
1020       Double_t x = tnp-dzdx*tnt;
1021       (GetLinearFitter(fDetectorPreviousTrack,kTRUE))->AddPoint(&x,dydt);
1022       if(fLinearFitterDebugOn) {
1023         fLinearVdriftFit->Update(fDetectorPreviousTrack,x,pars[1]);
1024       }
1025       fEntriesLinearFitter[fDetectorPreviousTrack]++;
1026     }
1027   }
1028
1029   return kTRUE;
1030 }
1031 //____________Offine tracking in the AliTRDtracker_____________________________
1032 Bool_t AliTRDCalibraFillHisto::HandlePRFtracklet(AliTRDtrack *t, Int_t index0, Int_t index1)
1033 {
1034   //
1035   // PRF width calibration
1036   // Assume a Gaussian shape: determinate the position of the three pad clusters
1037   // Fit with a straight line
1038   // Take the fitted values for all the clusters (3 or 2 pad clusters)
1039   // Fill the PRF as function of angle of the track
1040   //
1041   //
1042
1043
1044   //////////////////////////
1045   // variables
1046   /////////////////////////
1047   Int_t npoints  = index1-index0;                                           // number of total points
1048   Int_t nb3pc    = 0;                                                       // number of three pads clusters used for fit
1049   Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
1050   // To see the difference due to the fit
1051   Double_t *padPositions;
1052   padPositions = new Double_t[npoints];
1053   for(Int_t k = 0; k < npoints; k++){
1054     padPositions[k] = 0.0;
1055   }
1056   // Take the tgl and snp with the track t now
1057   Double_t  tgl = t->GetTglPlane(GetLayer(detector)); //dz/dl and not dz/dx
1058   Double_t  snp = t->GetSnpPlane(GetLayer(detector)); // sin angle in xy plan
1059   Float_t  dzdx = 0.0;                                // dzdx
1060   Float_t  tnp  = 0.0;
1061   if(TMath::Abs(snp) < 1.0){
1062     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1063     dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
1064   }
1065   // linear fitter
1066   fLinearFitterTracklet->ClearPoints();
1067
1068   ///////////////////////////
1069   // calcul the tnp group
1070   ///////////////////////////
1071   Bool_t echec   = kFALSE;
1072   Double_t shift = 0.0;
1073   //Calculate the shift in x coresponding to this tnp
1074   if(fNgroupprf != 0.0){
1075     shift      = -3.0*(fNgroupprf-1)-1.5;
1076     Double_t limithigh  = -0.2*(fNgroupprf-1);
1077     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1078     else{
1079       while(tnp > limithigh){
1080         limithigh += 0.2;
1081         shift += 3.0;
1082       }
1083     }
1084   }
1085   if(echec) {
1086     delete [] padPositions;
1087     return kFALSE;
1088   }
1089
1090   //////////////////////
1091   // loop clusters
1092   /////////////////////
1093   for(Int_t k = 0;  k < npoints; k++){
1094     //Take the cluster
1095     AliTRDcluster *cl  = (AliTRDcluster *) t->GetCluster(k+index0);
1096     Short_t  *signals  = cl->GetSignals();
1097     Double_t     time  = cl->GetLocalTimeBin();
1098     //Calculate x if possible
1099     Float_t xcenter    = 0.0;
1100     Bool_t  echec1      = kTRUE;
1101     if((time<=7) || (time>=21)) continue;
1102     // Center 3 balanced: position with the center of the pad
1103     if ((((Float_t) signals[3]) > 0.0) &&
1104         (((Float_t) signals[2]) > 0.0) &&
1105         (((Float_t) signals[4]) > 0.0)) {
1106       echec1 = kFALSE;
1107       // Security if the denomiateur is 0
1108       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1109            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1110         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1111           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1112                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1113       }
1114       else {
1115         echec1 = kTRUE;
1116       }
1117     }
1118     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
1119     if(echec) continue;
1120     //if no echec: calculate with the position of the pad
1121     // Position of the cluster
1122     Double_t       padPosition = xcenter +  cl->GetPadCol();
1123     padPositions[k]            = padPosition;
1124     nb3pc++;
1125     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1126   }//clusters loop
1127
1128
1129   /////////////////////////////
1130   // fit
1131   ////////////////////////////
1132   if(nb3pc < 3) {
1133     delete [] padPositions;
1134     fLinearFitterTracklet->ClearPoints();
1135     return kFALSE;
1136   }
1137   fLinearFitterTracklet->Eval();
1138   TVectorD line(2);
1139   fLinearFitterTracklet->GetParameters(line);
1140   Float_t  pointError  = -1.0;
1141   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1142     pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1143   }
1144   fLinearFitterTracklet->ClearPoints();
1145
1146   /////////////////////////////////////////////////////
1147   // Now fill the PRF: second loop over clusters
1148   /////////////////////////////////////////////////////
1149   for(Int_t k = 0;  k < npoints; k++){
1150     //Take the cluster
1151     AliTRDcluster *cl      = (AliTRDcluster *) t->GetCluster(k+index0);
1152     Short_t  *signals      = cl->GetSignals();              // signal
1153     Double_t     time      = cl->GetLocalTimeBin();              // time bin
1154     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1155     Float_t padPos         = cl->GetPadCol();               // middle pad
1156     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1157     Float_t ycenter        = 0.0;                           // relative center charge
1158     Float_t ymin           = 0.0;                           // relative left charge
1159     Float_t ymax           = 0.0;                           // relative right charge
1160
1161     //Requiere simply two pads clusters at least
1162     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1163        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1164       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1165       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1166       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1167       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1168     }
1169
1170     //calibration group
1171     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1172     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1173     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1174     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1175     Float_t   xcl          = cl->GetY();                                // y cluster
1176     Float_t   qcl          = cl->GetQ();                                // charge cluster
1177     Int_t     layer        = GetLayer(detector);                        // layer
1178     Int_t     stack        = GetStack(detector);                        // stack
1179     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1180     Double_t  x            = dpad;                                      // reconstructed position moved
1181     Float_t   ep           = pointError;                                // error of fit
1182     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1183     Float_t   signal3      = (Float_t)signals[3];                       // signal
1184     Float_t   signal2      = (Float_t)signals[2];                       // signal
1185     Float_t   signal4      = (Float_t)signals[4];                       // signal
1186     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1187
1188     //////////////////////////////
1189     // debug
1190     /////////////////////////////
1191     if(fDebugLevel > 0){
1192       if ( !fDebugStreamer ) {
1193         //debug stream
1194         TDirectory *backup = gDirectory;
1195         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1196         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1197       }
1198
1199
1200       x = xdiff;
1201       Int_t type=0;
1202       Float_t y = ycenter;
1203       (* fDebugStreamer) << "HandlePRFtracklet"<<
1204         "caligroup="<<caligroup<<
1205         "detector="<<detector<<
1206         "layer="<<layer<<
1207         "stack="<< stack <<
1208         "npoints="<<npoints<<
1209         "Np="<<nb3pc<<
1210         "ep="<<ep<<
1211         "type="<<type<<
1212         "snp="<<snp<<
1213         "tnp="<<tnp<<
1214         "tgl="<<tgl<<
1215         "dzdx="<<dzdx<<
1216         "padPos="<<padPos<<
1217         "padPosition="<<padPositions[k]<<
1218         "padPosTracklet="<<padPosTracklet<<
1219         "x="<<x<<
1220         "y="<<y<<
1221         "xcl="<<xcl<<
1222         "qcl="<<qcl<<
1223         "signal1="<<signal1<<
1224         "signal2="<<signal2<<
1225         "signal3="<<signal3<<
1226         "signal4="<<signal4<<
1227         "signal5="<<signal5<<
1228         "time="<<time<<
1229         "\n";
1230       x=-(xdiff+1);
1231       y = ymin;
1232       type=-1;
1233       (* fDebugStreamer) << "HandlePRFtracklet"<<
1234         "caligroup="<<caligroup<<
1235         "detector="<<detector<<
1236         "layer="<<layer<<
1237         "stack="<<stack<<
1238         "npoints="<<npoints<<
1239         "Np="<<nb3pc<<
1240         "ep="<<ep<<
1241         "type="<<type<<
1242         "snp="<<snp<<
1243         "tnp="<<tnp<<
1244         "tgl="<<tgl<<
1245         "dzdx="<<dzdx<<
1246         "padPos="<<padPos<<
1247         "padPosition="<<padPositions[k]<<
1248         "padPosTracklet="<<padPosTracklet<<
1249         "x="<<x<<
1250         "y="<<y<<
1251         "xcl="<<xcl<<
1252         "qcl="<<qcl<<
1253         "signal1="<<signal1<<
1254         "signal2="<<signal2<<
1255         "signal3="<<signal3<<
1256         "signal4="<<signal4<<
1257         "signal5="<<signal5<<
1258         "time="<<time<<
1259         "\n";
1260       x=1-xdiff;
1261       y = ymax;
1262       type=1;
1263       (* fDebugStreamer) << "HandlePRFtracklet"<<
1264         "caligroup="<<caligroup<<
1265         "detector="<<detector<<
1266         "layer="<<layer<<
1267         "stack="<<stack<<
1268         "npoints="<<npoints<<
1269         "Np="<<nb3pc<<
1270         "ep="<<ep<<
1271         "type="<<type<<
1272         "snp="<<snp<<
1273         "tnp="<<tnp<<
1274         "tgl="<<tgl<<
1275         "dzdx="<<dzdx<<
1276         "padPos="<<padPos<<
1277         "padPosition="<<padPositions[k]<<
1278         "padPosTracklet="<<padPosTracklet<<
1279         "x="<<x<<
1280         "y="<<y<<
1281         "xcl="<<xcl<<
1282         "qcl="<<qcl<<
1283         "signal1="<<signal1<<
1284         "signal2="<<signal2<<
1285         "signal3="<<signal3<<
1286         "signal4="<<signal4<<
1287         "signal5="<<signal5<<
1288         "time="<<time<<
1289         "\n";
1290
1291     }
1292
1293     ////////////////////////////
1294     // quality cuts
1295     ///////////////////////////
1296     if(npoints < fNumberClusters) continue;
1297     if(npoints > fNumberClustersf) continue;
1298     if(nb3pc <= 5) continue;
1299     if((time >= 21) || (time < 7)) continue;
1300     if(TMath::Abs(snp) >= 1.0) continue;
1301     if(TMath::Abs(qcl) < 80) continue;
1302
1303     ////////////////////////////
1304     // Fill
1305     ///////////////////////////
1306     if (fHisto2d) {
1307       if(TMath::Abs(dpad) < 1.5) {
1308         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1309         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1310       }
1311       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1312         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1313         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1314       }
1315       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1316         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1317         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1318       }
1319     }
1320     if (fVector2d) {
1321       if(TMath::Abs(dpad) < 1.5) {
1322         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1323         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1324       }
1325       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1326         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1327         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1328       }
1329       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1330         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1331         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1332       }
1333     }
1334   }
1335   delete [] padPositions;
1336   return kTRUE;
1337
1338 }
1339 //____________Offine tracking in the AliTRDtracker_____________________________
1340 Bool_t AliTRDCalibraFillHisto::HandlePRFtrackletV1(const AliTRDseedV1 *tracklet, Int_t nbclusters)
1341 {
1342   //
1343   // PRF width calibration
1344   // Assume a Gaussian shape: determinate the position of the three pad clusters
1345   // Fit with a straight line
1346   // Take the fitted values for all the clusters (3 or 2 pad clusters)
1347   // Fill the PRF as function of angle of the track
1348   //
1349   //
1350
1351   //printf("begin\n");
1352   ///////////////////////////////////////////
1353   // Take the parameters of the track
1354   //////////////////////////////////////////
1355   // take now the snp, tnp and tgl from the track
1356   Double_t snp = tracklet->GetSnp();             // sin dy/dx at the end of the chamber
1357   Double_t tnp = 0.0;                            // dy/dx at the end of the chamber
1358   if( TMath::Abs(snp) <  1.){
1359     tnp = snp / TMath::Sqrt((1.-snp)*(1.+snp));
1360   }
1361   Double_t tgl  = tracklet->GetTgl();           // dz/dl
1362   Double_t dzdx = tgl*TMath::Sqrt(1+tnp*tnp);   // dz/dx calculated from dz/dl
1363   // at the entrance
1364   //Double_t tnp = tracklet->GetYref(1);      // dy/dx at the entrance of the chamber
1365   //Double_t tgl = tracklet->GetZref(1);      // dz/dl at the entrance of the chamber
1366   //Double_t dzdx = tgl;                      //*TMath::Sqrt(1+tnp*tnp); // dz/dx from dz/dl
1367   // at the end with correction due to linear fit
1368   //Double_t tnp = tracklet->GetYfit(1);      // dy/dx at the end of the chamber after fit correction
1369   //Double_t tgl = tracklet->GetZfit(1);      // dz/dl at the end of the chamber after fit correction
1370
1371   ///////////////////////////////
1372   // Calculate tnp group shift
1373   ///////////////////////////////
1374   Bool_t echec   = kFALSE;
1375   Double_t shift = 0.0;
1376   //Calculate the shift in x coresponding to this tnp
1377   if(fNgroupprf != 0.0){
1378     shift      = -3.0*(fNgroupprf-1)-1.5;
1379     Double_t limithigh  = -0.2*(fNgroupprf-1);
1380     if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
1381     else{
1382       while(tnp > limithigh){
1383         limithigh += 0.2;
1384         shift += 3.0;
1385       }
1386     }
1387   }
1388   // do nothing if out of tnp range
1389   //printf("echec %d\n",(Int_t)echec);
1390   if(echec) return kFALSE;
1391
1392   ///////////////////////
1393   // Variables
1394   //////////////////////
1395
1396   Int_t nb3pc    = 0;              // number of three pads clusters used for fit
1397   // to see the difference between the fit and the 3 pad clusters position
1398   Double_t padPositions[AliTRDseedV1::kNtb];
1399   memset(padPositions, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
1400   fLinearFitterTracklet->ClearPoints();
1401
1402   //printf("loop clusters \n");
1403   ////////////////////////////
1404   // loop over the clusters
1405   ////////////////////////////
1406   AliTRDcluster *cl                   = 0x0;
1407   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1408     // reject shared clusters on pad row
1409     if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1410
1411     Double_t     time  = cl->GetLocalTimeBin();
1412     if((time<=7) || (time>=21)) continue;
1413     Short_t  *signals  = cl->GetSignals();
1414     Float_t xcenter    = 0.0;
1415     Bool_t  echec1      = kTRUE;
1416
1417     /////////////////////////////////////////////////////////////
1418     // Center 3 balanced: position with the center of the pad
1419     /////////////////////////////////////////////////////////////
1420     if ((((Float_t) signals[3]) > 0.0) &&
1421         (((Float_t) signals[2]) > 0.0) &&
1422         (((Float_t) signals[4]) > 0.0)) {
1423       echec1 = kFALSE;
1424       // Security if the denomiateur is 0
1425       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
1426            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
1427         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
1428           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
1429                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
1430       }
1431       else {
1432         echec1 = kTRUE;
1433       }
1434     }
1435     if(TMath::Abs(xcenter) > 0.5) echec1 = kTRUE;
1436     if(echec1) continue;
1437
1438     ////////////////////////////////////////////////////////
1439     //if no echec1: calculate with the position of the pad
1440     // Position of the cluster
1441     // fill the linear fitter
1442     ///////////////////////////////////////////////////////
1443     Double_t       padPosition = xcenter +  cl->GetPadCol();
1444     padPositions[ic]            = padPosition;
1445     nb3pc++;
1446     fLinearFitterTracklet->AddPoint(&time, padPosition,1);
1447
1448
1449   }//clusters loop
1450
1451   //printf("Fin loop clusters \n");
1452   //////////////////////////////
1453   // fit with a straight line
1454   /////////////////////////////
1455   if(nb3pc < 3){
1456     fLinearFitterTracklet->ClearPoints();
1457     return kFALSE;
1458   }
1459   fLinearFitterTracklet->Eval();
1460   TVectorD line(2);
1461   fLinearFitterTracklet->GetParameters(line);
1462   Float_t  pointError  = -1.0;
1463   if( fLinearFitterTracklet->GetChisquare()>=0.0) {
1464   pointError  =  TMath::Sqrt( fLinearFitterTracklet->GetChisquare()/(nb3pc-2));
1465   }
1466   fLinearFitterTracklet->ClearPoints();
1467
1468   //printf("PRF second loop \n");
1469   ////////////////////////////////////////////////
1470   // Fill the PRF: Second loop over clusters
1471   //////////////////////////////////////////////
1472   for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
1473     // reject shared clusters on pad row
1474     if(((ic+AliTRDseedV1::kNtb) < AliTRDseedV1::kNclusters) && (cl = tracklet->GetClusters(ic+AliTRDseedV1::kNtb))) continue;
1475     //
1476     if(!(cl = tracklet->GetClusters(ic))) continue;
1477
1478     Short_t  *signals      = cl->GetSignals();              // signal
1479     Double_t     time      = cl->GetLocalTimeBin();         // time bin
1480     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
1481     Float_t padPos         = cl->GetPadCol();               // middle pad
1482     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
1483     Float_t ycenter        = 0.0;                           // relative center charge
1484     Float_t ymin           = 0.0;                           // relative left charge
1485     Float_t ymax           = 0.0;                           // relative right charge
1486
1487     ////////////////////////////////////////////////////////////////
1488     // Calculate ycenter, ymin and ymax even for two pad clusters
1489     ////////////////////////////////////////////////////////////////
1490     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
1491        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
1492       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
1493       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
1494       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
1495       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
1496     }
1497
1498     /////////////////////////
1499     // Calibration group
1500     ////////////////////////
1501     Int_t     rowcl        = cl->GetPadRow();                           // row of cluster
1502     Int_t     colcl        = cl->GetPadCol();                           // col of cluster
1503     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcl,colcl);  // calcul the corresponding group
1504     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
1505     Float_t   xcl          = cl->GetY();                                // y cluster
1506     Float_t   qcl          = cl->GetQ();                                // charge cluster
1507     Int_t     layer        = GetLayer(fDetectorPreviousTrack);          // layer
1508     Int_t     stack        = GetStack(fDetectorPreviousTrack);          // stack
1509     Double_t  xdiff        = dpad;                                      // reconstructed position constant
1510     Double_t  x            = dpad;                                      // reconstructed position moved
1511     Float_t   ep           = pointError;                                // error of fit
1512     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
1513     Float_t   signal3      = (Float_t)signals[3];                       // signal
1514     Float_t   signal2      = (Float_t)signals[2];                       // signal
1515     Float_t   signal4      = (Float_t)signals[4];                       // signal
1516     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
1517
1518
1519
1520     /////////////////////
1521     // Debug stuff
1522     ////////////////////
1523
1524     if(fDebugLevel > 0){
1525       if ( !fDebugStreamer ) {
1526         //debug stream
1527         TDirectory *backup = gDirectory;
1528         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
1529         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1530       }
1531
1532       x = xdiff;
1533       Int_t type=0;
1534       Float_t y = ycenter;
1535       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1536         "caligroup="<<caligroup<<
1537         "detector="<<fDetectorPreviousTrack<<
1538         "layer="<<layer<<
1539         "stack="<<stack<<
1540         "npoints="<<nbclusters<<
1541         "Np="<<nb3pc<<
1542         "ep="<<ep<<
1543         "type="<<type<<
1544         "snp="<<snp<<
1545         "tnp="<<tnp<<
1546         "tgl="<<tgl<<
1547         "dzdx="<<dzdx<<
1548         "padPos="<<padPos<<
1549         "padPosition="<<padPositions[ic]<<
1550         "padPosTracklet="<<padPosTracklet<<
1551         "x="<<x<<
1552         "y="<<y<<
1553         "xcl="<<xcl<<
1554         "qcl="<<qcl<<
1555         "signal1="<<signal1<<
1556         "signal2="<<signal2<<
1557         "signal3="<<signal3<<
1558         "signal4="<<signal4<<
1559         "signal5="<<signal5<<
1560         "time="<<time<<
1561         "\n";
1562       x=-(xdiff+1);
1563       y = ymin;
1564       type=-1;
1565       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1566         "caligroup="<<caligroup<<
1567         "detector="<<fDetectorPreviousTrack<<
1568         "layer="<<layer<<
1569         "stack="<<stack<<
1570         "npoints="<<nbclusters<<
1571         "Np="<<nb3pc<<
1572         "ep="<<ep<<
1573         "type="<<type<<
1574         "snp="<<snp<<
1575         "tnp="<<tnp<<
1576         "tgl="<<tgl<<
1577         "dzdx="<<dzdx<<
1578         "padPos="<<padPos<<
1579         "padPosition="<<padPositions[ic]<<
1580         "padPosTracklet="<<padPosTracklet<<
1581         "x="<<x<<
1582         "y="<<y<<
1583         "xcl="<<xcl<<
1584         "qcl="<<qcl<<
1585         "signal1="<<signal1<<
1586         "signal2="<<signal2<<
1587         "signal3="<<signal3<<
1588         "signal4="<<signal4<<
1589         "signal5="<<signal5<<
1590         "time="<<time<<
1591         "\n";
1592       x=1-xdiff;
1593       y = ymax;
1594       type=1;
1595       (* fDebugStreamer) << "HandlePRFtrackletV1"<<
1596         "caligroup="<<caligroup<<
1597         "detector="<<fDetectorPreviousTrack<<
1598         "layer="<<layer<<
1599         "stack="<<stack<<
1600         "npoints="<<nbclusters<<
1601         "Np="<<nb3pc<<
1602         "ep="<<ep<<
1603         "type="<<type<<
1604         "snp="<<snp<<
1605         "tnp="<<tnp<<
1606         "tgl="<<tgl<<
1607         "dzdx="<<dzdx<<
1608         "padPos="<<padPos<<
1609         "padPosition="<<padPositions[ic]<<
1610         "padPosTracklet="<<padPosTracklet<<
1611         "x="<<x<<
1612         "y="<<y<<
1613         "xcl="<<xcl<<
1614         "qcl="<<qcl<<
1615         "signal1="<<signal1<<
1616         "signal2="<<signal2<<
1617         "signal3="<<signal3<<
1618         "signal4="<<signal4<<
1619         "signal5="<<signal5<<
1620         "time="<<time<<
1621         "\n";
1622
1623     }
1624
1625     /////////////////////
1626     // Cuts quality
1627     /////////////////////
1628     if(nbclusters < fNumberClusters) continue;
1629     if(nbclusters > fNumberClustersf) continue;
1630     if(nb3pc <= 5) continue;
1631     if((time >= 21) || (time < 7)) continue;
1632     if(TMath::Abs(qcl) < 80) continue;
1633     if( TMath::Abs(snp) >  1.) continue;
1634
1635
1636     ////////////////////////
1637     // Fill the histos
1638     ///////////////////////
1639     if (fHisto2d) {
1640       if(TMath::Abs(dpad) < 1.5) {
1641         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
1642         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
1643         //printf("place %f, ycenter %f\n",(shift+dpad),ycenter);
1644       }
1645       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1646         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
1647         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
1648       }
1649       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
1650         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
1651         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
1652       }
1653     }
1654     // vector method
1655     if (fVector2d) {
1656       if(TMath::Abs(dpad) < 1.5) {
1657         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
1658         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
1659       }
1660       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
1661         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
1662         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
1663       }
1664       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
1665         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
1666         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
1667       }
1668     }
1669   } // second loop over clusters
1670
1671
1672   return kTRUE;
1673 }
1674 ///////////////////////////////////////////////////////////////////////////////////////
1675 // Pad row col stuff: see if masked or not
1676 ///////////////////////////////////////////////////////////////////////////////////////
1677 //_____________________________________________________________________________
1678 void AliTRDCalibraFillHisto::CheckGoodTrackletV1(AliTRDcluster *cl)
1679 {
1680   //
1681   // See if we are not near a masked pad
1682   //
1683
1684   if(cl->IsMasked()) fGoodTracklet = kFALSE;
1685
1686
1687 }
1688 //_____________________________________________________________________________
1689 void AliTRDCalibraFillHisto::CheckGoodTrackletV0(Int_t detector, Int_t row, Int_t col)
1690 {
1691   //
1692   // See if we are not near a masked pad
1693   //
1694
1695   if (!IsPadOn(detector, col, row)) {
1696     fGoodTracklet = kFALSE;
1697   }
1698
1699   if (col > 0) {
1700     if (!IsPadOn(detector, col-1, row)) {
1701       fGoodTracklet = kFALSE;
1702     }
1703   }
1704
1705   if (col < 143) {
1706     if (!IsPadOn(detector, col+1, row)) {
1707       fGoodTracklet = kFALSE;
1708     }
1709   }
1710
1711 }
1712 //_____________________________________________________________________________
1713 Bool_t AliTRDCalibraFillHisto::IsPadOn(Int_t detector, Int_t row, Int_t col) const
1714 {
1715   //
1716   // Look in the choosen database if the pad is On.
1717   // If no the track will be "not good"
1718   //
1719
1720   // Get the parameter object
1721   AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
1722   if (!cal) {
1723     AliInfo("Could not get calibDB");
1724     return kFALSE;
1725   }
1726
1727   if (!cal->IsChamberInstalled(detector)     ||
1728        cal->IsChamberMasked(detector)        ||
1729        cal->IsPadMasked(detector,col,row)) {
1730     return kFALSE;
1731   }
1732   else {
1733     return kTRUE;
1734   }
1735
1736 }
1737 ///////////////////////////////////////////////////////////////////////////////////////
1738 // Calibration groups: calculate the number of groups, localise...
1739 ////////////////////////////////////////////////////////////////////////////////////////
1740 //_____________________________________________________________________________
1741 Int_t AliTRDCalibraFillHisto::CalculateCalibrationGroup(Int_t i, Int_t row, Int_t col) const
1742 {
1743   //
1744   // Calculate the calibration group number for i
1745   //
1746
1747   // Row of the cluster and position in the pad groups
1748   Int_t posr = 0;
1749   if (fCalibraMode->GetNnZ(i) != 0) {
1750     posr = (Int_t) row / fCalibraMode->GetNnZ(i);
1751   }
1752
1753
1754   // Col of the cluster and position in the pad groups
1755   Int_t posc = 0;
1756   if (fCalibraMode->GetNnRphi(i) != 0) {
1757     posc = (Int_t) col / fCalibraMode->GetNnRphi(i);
1758   }
1759
1760   return posc*fCalibraMode->GetNfragZ(i)+posr;
1761
1762 }
1763 //____________________________________________________________________________________
1764 Int_t AliTRDCalibraFillHisto::CalculateTotalNumberOfBins(Int_t i)
1765 {
1766   //
1767   // Calculate the total number of calibration groups
1768   //
1769
1770   Int_t ntotal = 0;
1771
1772   // All together
1773   if((fCalibraMode->GetNz(i)==100) && (fCalibraMode->GetNrphi(i)==100)){
1774     ntotal = 1;
1775     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1776     return ntotal;
1777   }
1778
1779   // Per Supermodule
1780   if((fCalibraMode->GetNz(i)==10) && (fCalibraMode->GetNrphi(i)==10)){
1781     ntotal = 18;
1782     AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1783     return ntotal;
1784   }
1785
1786   // More
1787   fCalibraMode->ModePadCalibration(2,i);
1788   fCalibraMode->ModePadFragmentation(0,2,0,i);
1789   fCalibraMode->SetDetChamb2(i);
1790   ntotal += 6 * 18 * fCalibraMode->GetDetChamb2(i);
1791   fCalibraMode->ModePadCalibration(0,i);
1792   fCalibraMode->ModePadFragmentation(0,0,0,i);
1793   fCalibraMode->SetDetChamb0(i);
1794   ntotal += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(i);
1795   AliInfo(Form("Total number of Xbins: %d for i %d",ntotal,i));
1796   return ntotal;
1797
1798 }
1799 //_____________________________________________________________________________
1800 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1801 {
1802   //
1803   // Set the mode of calibration group in the z direction for the parameter i
1804   //
1805
1806   if ((Nz >= 0) &&
1807       (Nz <  5)) {
1808     fCalibraMode->SetNz(i, Nz);
1809   }
1810   else {
1811     AliInfo("You have to choose between 0 and 4");
1812   }
1813
1814 }
1815 //_____________________________________________________________________________
1816 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1817 {
1818   //
1819   // Set the mode of calibration group in the rphi direction for the parameter i
1820   //
1821
1822   if ((Nrphi >= 0) &&
1823       (Nrphi <  7)) {
1824     fCalibraMode->SetNrphi(i ,Nrphi);
1825   }
1826   else {
1827     AliInfo("You have to choose between 0 and 6");
1828   }
1829
1830 }
1831
1832 //_____________________________________________________________________________
1833 void AliTRDCalibraFillHisto::SetAllTogether(Int_t i)
1834 {
1835   //
1836   // Set the mode of calibration group all together
1837   //
1838   if(fVector2d == kTRUE) {
1839     AliInfo("Can not work with the vector method");
1840     return;
1841   }
1842   fCalibraMode->SetAllTogether(i);
1843
1844 }
1845
1846 //_____________________________________________________________________________
1847 void AliTRDCalibraFillHisto::SetPerSuperModule(Int_t i)
1848 {
1849   //
1850   // Set the mode of calibration group per supermodule
1851   //
1852   if(fVector2d == kTRUE) {
1853     AliInfo("Can not work with the vector method");
1854     return;
1855   }
1856   fCalibraMode->SetPerSuperModule(i);
1857
1858 }
1859
1860 //____________Set the pad calibration variables for the detector_______________
1861 Bool_t AliTRDCalibraFillHisto::LocalisationDetectorXbins(Int_t detector)
1862 {
1863   //
1864   // For the detector calcul the first Xbins and set the number of row
1865   // and col pads per calibration groups, the number of calibration
1866   // groups in the detector.
1867   //
1868
1869   // first Xbins of the detector
1870   if (fCH2dOn) {
1871     fCalibraMode->CalculXBins(detector,0);
1872   }
1873   if (fPH2dOn) {
1874     fCalibraMode->CalculXBins(detector,1);
1875   }
1876   if (fPRF2dOn) {
1877     fCalibraMode->CalculXBins(detector,2);
1878   }
1879
1880   // fragmentation of idect
1881   for (Int_t i = 0; i < 3; i++) {
1882     fCalibraMode->ModePadCalibration((Int_t) GetStack(detector),i);
1883     fCalibraMode->ModePadFragmentation((Int_t) GetLayer(detector)
1884                        , (Int_t) GetStack(detector)
1885                        , (Int_t) GetSector(detector),i);
1886   }
1887
1888   return kTRUE;
1889
1890 }
1891 //_____________________________________________________________________________
1892 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1893 {
1894   //
1895   // Should be between 0 and 6
1896   //
1897
1898   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1899     AliInfo("The number of groups must be between 0 and 6!");
1900   }
1901   else {
1902     fNgroupprf = numberGroupsPRF;
1903   }
1904
1905 }
1906 ///////////////////////////////////////////////////////////////////////////////////////////
1907 // Per tracklet: store or reset the info, fill the histos with the info
1908 //////////////////////////////////////////////////////////////////////////////////////////
1909 //_____________________________________________________________________________
1910 void AliTRDCalibraFillHisto::StoreInfoCHPHtrack(AliTRDcluster *cl, Double_t dqdl, Int_t *group, Int_t row, Int_t col)
1911 {
1912   //
1913   // Store the infos in fAmpTotal, fPHPlace and fPHValue
1914   // Correct from the gain correction before
1915   //
1916
1917   // time bin of the cluster not corrected
1918   Int_t    time     = cl->GetPadTime();
1919
1920   //Correct for the gain coefficient used in the database for reconstruction
1921   Float_t correctthegain = 1.0;
1922   if(fIsHLT) correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack);
1923   else correctthegain = fCalDetGain->GetValue(fDetectorPreviousTrack)*fCalROCGain->GetValue(col,row);
1924   Float_t correction    = 1.0;
1925   Float_t normalisation = 6.67;
1926   // we divide with gain in AliTRDclusterizer::Transform...
1927   if( correctthegain > 0 ) normalisation /= correctthegain;
1928
1929
1930   // take dd/dl corrected from the angle of the track
1931   correction = dqdl / (normalisation);
1932
1933
1934   // Fill the fAmpTotal with the charge
1935   if (fCH2dOn) {
1936     if((!fLimitChargeIntegration) || (cl->IsInChamber())) fAmpTotal[(Int_t) group[0]] += correction;
1937   }
1938
1939   // Fill the fPHPlace and value
1940   if (fPH2dOn) {
1941     fPHPlace[time] = group[1];
1942     fPHValue[time] = correction;
1943   }
1944
1945 }
1946 //____________Offine tracking in the AliTRDtracker_____________________________
1947 void AliTRDCalibraFillHisto::ResetfVariablestracklet()
1948 {
1949   //
1950   // Reset values per tracklet
1951   //
1952
1953   //Reset good tracklet
1954   fGoodTracklet = kTRUE;
1955
1956   // Reset the fPHValue
1957   if (fPH2dOn) {
1958     //Reset the fPHValue and fPHPlace
1959     for (Int_t k = 0; k < fTimeMax; k++) {
1960       fPHValue[k] = 0.0;
1961       fPHPlace[k] = -1;
1962     }
1963   }
1964
1965   // Reset the fAmpTotal where we put value
1966   if (fCH2dOn) {
1967     for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1968       fAmpTotal[k] = 0.0;
1969     }
1970   }
1971 }
1972 //____________Offine tracking in the AliTRDtracker_____________________________
1973 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH(Int_t nbclusters)
1974 {
1975   //
1976   // For the offline tracking or mcm tracklets
1977   // This function will be called in the functions UpdateHistogram...
1978   // to fill the info of a track for the relativ gain calibration
1979   //
1980
1981   Int_t nb            =  0;   // Nombre de zones traversees
1982   Int_t fd            = -1;   // Premiere zone non nulle
1983   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1984
1985   if(nbclusters < fNumberClusters) return;
1986   if(nbclusters > fNumberClustersf) return;
1987
1988
1989   // Normalize with the number of clusters
1990   Double_t normalizeCst = fRelativeScale;
1991   if(fNormalizeNbOfCluster) normalizeCst = normalizeCst*nbclusters;
1992
1993   // See if the track goes through different zones
1994   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1995     if (fAmpTotal[k] > 0.0) {
1996       totalcharge += fAmpTotal[k];
1997       nb++;
1998       if (nb == 1) {
1999         fd = k;
2000       }
2001     }
2002   }
2003
2004
2005   switch (nb)
2006     {
2007     case 1:
2008       fNumberUsedCh[0]++;
2009       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2010       if (fHisto2d) {
2011         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2012         //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2013       }
2014       if (fVector2d) {
2015         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2016       }
2017       break;
2018     case 2:
2019       if ((fAmpTotal[fd]   > 0.0) &&
2020           (fAmpTotal[fd+1] > 0.0)) {
2021         // One of the two very big
2022         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
2023           if (fHisto2d) {
2024             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2025             //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2026           }
2027           if (fVector2d) {
2028             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2029           }
2030           fNumberUsedCh[1]++;
2031           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2032         }
2033         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
2034           if (fHisto2d) {
2035             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/normalizeCst);
2036             //fCH2d->Fill(fAmpTotal[fd+1]/normalizeCst,fCalibraMode->GetXbins(0)+fd+1.5);
2037           }
2038           if (fVector2d) {
2039             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/normalizeCst);
2040           }
2041           fNumberUsedCh[1]++;
2042           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
2043         }
2044       }
2045       if (fCalibraMode->GetNfragZ(0) > 1) {
2046         if (fAmpTotal[fd] > 0.0) {
2047           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
2048             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
2049               // One of the two very big
2050               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
2051                 if (fHisto2d) {
2052                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/normalizeCst);
2053                   //fCH2d->Fill(fAmpTotal[fd]/normalizeCst,fCalibraMode->GetXbins(0)+fd+0.5);
2054                 }
2055                 if (fVector2d) {
2056                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/normalizeCst);
2057                 }
2058                 fNumberUsedCh[1]++;
2059                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
2060               }
2061               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
2062                 if (fHisto2d) {
2063                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2064                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
2065                 }
2066                 fNumberUsedCh[1]++;
2067                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
2068                 if (fVector2d) {
2069                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/normalizeCst);
2070                 }
2071               }
2072             }
2073           }
2074         }
2075       }
2076       break;
2077     default: break;
2078     }
2079 }
2080 //____________Offine tracking in the AliTRDtracker_____________________________
2081 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
2082 {
2083   //
2084   // For the offline tracking or mcm tracklets
2085   // This function will be called in the functions UpdateHistogram...
2086   // to fill the info of a track for the drift velocity  calibration
2087   //
2088
2089   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
2090   Int_t fd1 = -1; // Premiere zone non nulle
2091   Int_t fd2 = -1; // Deuxieme zone non nulle
2092   Int_t k1  = -1; // Debut de la premiere zone
2093   Int_t k2  = -1; // Debut de la seconde zone
2094   Int_t nbclusters = 0; // number of clusters
2095
2096
2097
2098   // See if the track goes through different zones
2099   for (Int_t k = 0; k < fTimeMax; k++) {
2100     if (fPHValue[k] > 0.0) {
2101       nbclusters++;
2102       if (fd1 == -1) {
2103         fd1 = fPHPlace[k];
2104         k1  = k;
2105       }
2106       if (fPHPlace[k] != fd1) {
2107         if (fd2 == -1) {
2108           k2  = k;
2109           fd2 = fPHPlace[k];
2110           nb  = 2;
2111         }
2112         if (fPHPlace[k] != fd2) {
2113           nb = 3;
2114         }
2115       }
2116     }
2117   }
2118
2119   // See if noise before and after
2120   if(fMaxCluster > 0) {
2121     if(fPHValue[0] > fMaxCluster) return;
2122     if(fTimeMax > fNbMaxCluster) {
2123       for(Int_t k = (fTimeMax-fNbMaxCluster); k < fTimeMax; k++){
2124         if(fPHValue[k] > fMaxCluster) return;
2125       }
2126     }
2127   }
2128
2129   //printf("nbclusters %d, low limit %d, high limit %d\n",nbclusters,fNumberClusters,fNumberClustersf);
2130
2131   if(nbclusters < fNumberClusters) return;
2132   if(nbclusters > fNumberClustersf) return;
2133
2134   switch(nb)
2135     {
2136     case 1:
2137       fNumberUsedPh[0]++;
2138       for (Int_t i = 0; i < fTimeMax; i++) {
2139         if (fHisto2d) {
2140           if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2141           else {
2142             if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2143               }
2144           //printf("Fill the time bin %d with %f\n",i,fPHValue[i]);
2145         }
2146         if (fVector2d) {
2147           if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2148           else {
2149             if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2150           }
2151         }
2152       }
2153       break;
2154     case 2:
2155       if ((fd1 == fd2+1) ||
2156           (fd2 == fd1+1)) {
2157         // One of the two fast all the think
2158         if (k2 > (k1+fDifference)) {
2159           //we choose to fill the fd1 with all the values
2160           fNumberUsedPh[1]++;
2161           for (Int_t i = 0; i < fTimeMax; i++) {
2162             if (fHisto2d) {
2163               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2164               else {
2165                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2166                   }
2167             }
2168             if (fVector2d) {
2169               if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2170               else {
2171                 if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2172                   }
2173             }
2174           }
2175         }
2176         if ((k2+fDifference) < fTimeMax) {
2177           //we choose to fill the fd2 with all the values
2178           fNumberUsedPh[1]++;
2179           for (Int_t i = 0; i < fTimeMax; i++) {
2180             if (fHisto2d) {
2181               if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2182               else {
2183                 if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2184               }
2185             }
2186           if (fVector2d) {
2187             if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2188             else {
2189               if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2190             }
2191           }
2192           }
2193         }
2194       }
2195       // Two zones voisines sinon rien!
2196       if (fCalibraMode->GetNfragZ(1) > 1) {
2197         // Case 2
2198         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2199           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2200             // One of the two fast all the think
2201             if (k2 > (k1+fDifference)) {
2202               //we choose to fill the fd1 with all the values
2203               fNumberUsedPh[1]++;
2204               for (Int_t i = 0; i < fTimeMax; i++) {
2205                 if (fHisto2d) {
2206                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2207                   else {
2208                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2209                   }
2210                 }
2211                 if (fVector2d) {
2212                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2213                   else {
2214                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2215                   }
2216                 }
2217               }
2218             }
2219             if ((k2+fDifference) < fTimeMax) {
2220               //we choose to fill the fd2 with all the values
2221               fNumberUsedPh[1]++;
2222               for (Int_t i = 0; i < fTimeMax; i++) {
2223                 if (fHisto2d) {
2224                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2225                   else {
2226                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2227                   }
2228                 }
2229                 if (fVector2d) {
2230                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2231                   else {
2232                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2233                   }
2234                 }
2235               }
2236             }
2237           }
2238         }
2239         // Two zones voisines sinon rien!
2240         // Case 3
2241         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2242           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2243             // One of the two fast all the think
2244             if (k2 > (k1 + fDifference)) {
2245               //we choose to fill the fd1 with all the values
2246               fNumberUsedPh[1]++;
2247               for (Int_t i = 0; i < fTimeMax; i++) {
2248                 if (fHisto2d) {
2249                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2250                   else {
2251                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2252                   }
2253                 }
2254                 if (fVector2d) {
2255                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2256                   else {
2257                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2258                   }
2259                 }
2260               }
2261             }
2262             if ((k2+fDifference) < fTimeMax) {
2263               //we choose to fill the fd2 with all the values
2264               fNumberUsedPh[1]++;
2265               for (Int_t i = 0; i < fTimeMax; i++) {
2266                 if (fHisto2d) {
2267                   if(fFillWithZero) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2268                   else {
2269                     if(((Float_t) fPHValue[i] > 0.0)) fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2270                   }
2271                 }
2272                 if (fVector2d) {
2273                   if(fFillWithZero) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2274                   else {
2275                     if(((Float_t) fPHValue[i] > 0.0)) fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2276                   }
2277                 }
2278               }
2279             }
2280           }
2281         }
2282       }
2283       break;
2284     default: break;
2285     }
2286 }
2287 //////////////////////////////////////////////////////////////////////////////////////////
2288 // DAQ process functions
2289 /////////////////////////////////////////////////////////////////////////////////////////
2290 //_____________________________________________________________________
2291 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliTRDrawStreamBase *rawStream, Bool_t nocheck)
2292 {
2293   //
2294   // Event Processing loop - AliTRDrawStreamBase
2295   // TestBeam 2007 version
2296   // 0 timebin problem
2297   // 1 no input
2298   // 2 input
2299   //
2300   // Same algorithm as TestBeam but different reader
2301   //
2302
2303   Int_t withInput = 1;
2304
2305   Double_t phvalue[16][144][36];
2306   for(Int_t k = 0; k < 36; k++){
2307     for(Int_t j = 0; j < 16; j++){
2308       for(Int_t c = 0; c < 144; c++){
2309         phvalue[j][c][k] = 0.0;
2310       }
2311     }
2312   }
2313
2314   fDetectorPreviousTrack = -1;
2315   fMCMPrevious           = -1;
2316   fROBPrevious           = -1;
2317   Int_t nbtimebin = 0;
2318   Int_t baseline  = 10;
2319
2320
2321   if(!nocheck){
2322
2323     fTimeMax = 0;
2324
2325     while (rawStream->Next()) {
2326
2327       Int_t idetector = rawStream->GetDet();                            //  current detector
2328       Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2329       Int_t irob      = rawStream->GetROB();                            //  current ROB
2330
2331       //printf("Detector %d\n",idetector);
2332
2333       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2334
2335         // Fill
2336         withInput = TMath::Max(FillDAQ(phvalue),withInput);
2337
2338
2339         // reset
2340         for(Int_t k = 0; k < 36; k++){
2341           for(Int_t j = 0; j < 16; j++){
2342             for(Int_t c = 0; c < 144; c++){
2343               phvalue[j][c][k] = 0.0;
2344             }
2345           }
2346         }
2347       }
2348
2349       fDetectorPreviousTrack = idetector;
2350       fMCMPrevious           = imcm;
2351       fROBPrevious           = irob;
2352
2353       nbtimebin              = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2354       if(nbtimebin == 0) return 0;
2355       if((fTimeMax != 0) && (nbtimebin != fTimeMax)) return 0;
2356       fTimeMax          = nbtimebin;
2357
2358       //baseline          = rawStream->GetCommonAdditive();                // common additive baseline
2359       fNumberClustersf    = fTimeMax;
2360       fNumberClusters     = (Int_t)(0.6*fTimeMax);
2361
2362
2363       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
2364       Int_t  col        = rawStream->GetCol();
2365       Int_t row         = rawStream->GetRow();
2366
2367
2368       //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2369
2370
2371       for(Int_t itime = 0; itime < nbtimebin; itime++){
2372         phvalue[row][col][itime] = signal[itime]-baseline;
2373       }
2374     }
2375
2376     // fill the last one
2377     if(fDetectorPreviousTrack != -1){
2378
2379       // Fill
2380       withInput = TMath::Max(FillDAQ(phvalue),withInput);
2381
2382       // reset
2383       for(Int_t k = 0; k < 36; k++){
2384         for(Int_t j = 0; j < 16; j++){
2385           for(Int_t c = 0; c < 144; c++){
2386             phvalue[j][c][k] = 0.0;
2387           }
2388         }
2389       }
2390     }
2391
2392   }
2393   else{
2394
2395     while (rawStream->Next()) {
2396
2397       Int_t idetector = rawStream->GetDet();                            //  current detector
2398       Int_t imcm      = rawStream->GetMCM();                            //  current MCM
2399       Int_t irob      = rawStream->GetROB();                            //  current ROB
2400
2401       //printf("Detector %d\n",idetector);
2402
2403       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
2404
2405         // Fill
2406         withInput = TMath::Max(FillDAQ(phvalue),withInput);
2407
2408         // reset
2409         for(Int_t k = 0; k < 36; k++){
2410           for(Int_t j = 0; j < 16; j++){
2411             for(Int_t c = 0; c < 144; c++){
2412               phvalue[j][c][k] = 0.0;
2413             }
2414           }
2415         }
2416       }
2417
2418       fDetectorPreviousTrack = idetector;
2419       fMCMPrevious           = imcm;
2420       fROBPrevious           = irob;
2421
2422       //baseline          = rawStream->GetCommonAdditive();                //  common baseline
2423
2424       fTimeMax          = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
2425       fNumberClustersf    = fTimeMax;
2426       fNumberClusters     = (Int_t)(0.6*fTimeMax);
2427       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
2428       Int_t col         = rawStream->GetCol();
2429       Int_t row         = rawStream->GetRow();
2430
2431
2432       //printf("detector %d, signal[0] %d, signal[1] %d, signal[2] %d, baseline %d\n",idetector,signal[0],signal[1],signal[2], baseline);
2433
2434       for(Int_t itime = 0; itime < fTimeMax; itime++){
2435         phvalue[row][col][itime] = signal[itime]-baseline;
2436       }
2437     }
2438
2439     // fill the last one
2440     if(fDetectorPreviousTrack != -1){
2441
2442       // Fill
2443       withInput = TMath::Max(FillDAQ(phvalue),withInput);
2444
2445       // reset
2446       for(Int_t k = 0; k < 36; k++){
2447         for(Int_t j = 0; j < 16; j++){
2448           for(Int_t c = 0; c < 144; c++){
2449             phvalue[j][c][k] = 0.0;
2450           }
2451         }
2452       }
2453     }
2454   }
2455
2456   return withInput;
2457
2458 }
2459 //_____________________________________________________________________
2460 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
2461 {
2462   //
2463   //  Event processing loop - AliRawReader
2464   //  Testbeam 2007 version
2465   //
2466
2467   AliTRDrawStreamBase rawStream(rawReader);
2468
2469   rawReader->Select("TRD");
2470
2471   return ProcessEventDAQ(&rawStream, nocheck);
2472 }
2473
2474 //_________________________________________________________________________
2475 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
2476 #ifdef ALI_DATE
2477                                                eventHeaderStruct *event,
2478                                                Bool_t nocheck
2479 #else
2480                                                eventHeaderStruct* /*event*/,
2481                                                Bool_t /*nocheck*/
2482
2483 #endif
2484                                    )
2485 {
2486   //
2487   //  process date event
2488   //  Testbeam 2007 version
2489   //
2490 #ifdef ALI_DATE
2491     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
2492     Int_t result=ProcessEventDAQ(rawReader, nocheck);
2493     delete rawReader;
2494     return result;
2495 #else
2496     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
2497     return 0;
2498 #endif
2499
2500 }
2501 //////////////////////////////////////////////////////////////////////////////
2502 // Routine inside the DAQ process
2503 /////////////////////////////////////////////////////////////////////////////
2504 //_______________________________________________________________________
2505 Int_t AliTRDCalibraFillHisto::FillDAQ(Double_t phvalue[16][144][36]){
2506
2507   //
2508   // Look for the maximum by collapsing over the time
2509   // Sum over four pad col and two pad row
2510   //
2511
2512   Int_t used = 0;
2513
2514
2515   Int_t idect = fDetectorPreviousTrack;
2516   //printf("Enter Detector %d\n",fDetectorPreviousTrack);
2517   Double_t sum[36];
2518   for(Int_t tb = 0; tb < 36; tb++){
2519     sum[tb] = 0.0;
2520   }
2521
2522   //fGoodTracklet = kTRUE;
2523   //fDetectorPreviousTrack = 0;
2524
2525
2526   ///////////////////////////
2527   // look for maximum
2528   /////////////////////////
2529
2530   Int_t imaxRow = 0;
2531   Int_t imaxCol = 0;
2532   Double_t integralMax = -1;
2533
2534   for (Int_t ir = 1; ir <= 15; ir++)
2535     {
2536       for (Int_t ic = 2; ic <= 142; ic++)
2537         {
2538           Double_t integral = 0;
2539           for (Int_t ishiftR = 0; ishiftR < 2; ishiftR++)
2540             {
2541               for (Int_t ishiftC = -2; ishiftC < 2; ishiftC++)
2542                 {
2543                   if (ir + ishiftR >= 1 && ir + ishiftR <= 16 &&
2544                       ic + ishiftC >= 1 && ic + ishiftC <= 144)
2545                     {
2546
2547                       for(Int_t tb = 0; tb< fTimeMax; tb++){
2548                         integral += phvalue[ir + ishiftR-1][ic + ishiftC-1][tb];
2549                       }// addtb
2550                     } //addsignal
2551                 } //shiftC
2552             } // shiftR
2553
2554           if (integralMax < integral)
2555             {
2556               imaxRow = ir;
2557               imaxCol = ic;
2558               integralMax = integral;
2559             } // check max integral
2560         } //ic
2561     } // ir
2562
2563   //printf("imaxRow %d, imaxCol %d, fTimeMax %d, integralMax %f\n",imaxRow,imaxCol,fTimeMax, integralMax);
2564
2565   if((imaxRow == 0) || (imaxCol == 0)) {
2566     used=1;
2567     return used;
2568   }
2569   //CheckGoodTrackletV0(fDetectorPreviousTrack,imaxRow,imaxCol);
2570   //if(!fGoodTracklet) used = 1;;
2571
2572   //  /////////////////////////////////////////////////////
2573   // sum ober 2 row and 4 pad cols for each time bins
2574   //  ////////////////////////////////////////////////////
2575
2576
2577   for (Int_t ir = imaxRow - 1; ir < imaxRow + 1; ir++)
2578     {
2579       for (Int_t ic = imaxCol - 2; ic < imaxCol + 2; ic++)
2580         {
2581           for(Int_t it = 0; it < fTimeMax; it++){
2582             sum[it] += phvalue[ir][ic][it];
2583           }
2584         }//ic
2585     }//ir
2586
2587   Int_t nbcl = 0;
2588   Double_t sumcharge = 0.0;
2589   for(Int_t it = 0; it < fTimeMax; it++){
2590     sumcharge += sum[it];
2591     if(sum[it] > 20.0) nbcl++;
2592   }
2593
2594
2595   /////////////////////////////////////////////////////////
2596   // Debug
2597   ////////////////////////////////////////////////////////
2598   if(fDebugLevel > 0){
2599     if ( !fDebugStreamer ) {
2600       //debug stream
2601       TDirectory *backup = gDirectory;
2602       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2603       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2604     }
2605
2606     Double_t amph0 = sum[0];
2607     Double_t amphlast = sum[fTimeMax-1];
2608     Double_t rms      = TMath::RMS(fTimeMax,sum);
2609     Int_t    goodtracklet = (Int_t) fGoodTracklet;
2610     for(Int_t it = 0; it < fTimeMax; it++){
2611       Double_t clustera = sum[it];
2612
2613     (* fDebugStreamer) << "FillDAQa"<<
2614       "ampTotal="<<sumcharge<<
2615       "row="<<imaxRow<<
2616       "col="<<imaxCol<<
2617       "detector="<<idect<<
2618       "amph0="<<amph0<<
2619       "amphlast="<<amphlast<<
2620       "goodtracklet="<<goodtracklet<<
2621       "clustera="<<clustera<<
2622       "it="<<it<<
2623       "rms="<<rms<<
2624       "\n";
2625     }
2626   }
2627
2628   ////////////////////////////////////////////////////////
2629   // fill
2630   ///////////////////////////////////////////////////////
2631   if(sum[0] > 100.0) used = 1;
2632   if(nbcl < fNumberClusters) used = 1;
2633   if(nbcl > fNumberClustersf) used = 1;
2634
2635   //if(fDetectorPreviousTrack == 15){
2636   //  printf("rms %f and first time bin %f\n",TMath::RMS(fTimeMax,sum),sum[0]);
2637   //}
2638   //if((TMath::RMS(fTimeMax,sum) <= 10.0) && (sum[0] > 200.0)) return 1;
2639   if(used == 0){
2640     for(Int_t it = 0; it < fTimeMax; it++){
2641       if(fFillWithZero) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2642       else{
2643         if(sum[it] > 0.0) UpdateDAQ(fDetectorPreviousTrack,0,0,it,sum[it],fTimeMax);
2644       }
2645     }
2646
2647
2648     //((TH2I *)GetCH2d()->Fill(sumcharge/30.0,fDetectorPreviousTrack));
2649     used = 2;
2650     //printf("Pass Detector %d\n",fDetectorPreviousTrack);
2651
2652   }
2653
2654   return used;
2655
2656 }
2657 //____________Online trackling in AliTRDtrigger________________________________
2658 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Float_t signal, Int_t nbtimebins)
2659 {
2660   //
2661   // For the DAQ
2662   // Fill a simple average pulse height
2663   //
2664
2665
2666   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
2667
2668
2669   return kTRUE;
2670
2671 }
2672 //____________Write_____________________________________________________
2673 //_____________________________________________________________________
2674 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
2675 {
2676   //
2677   //  Write infos to file
2678   //
2679
2680   //For debugging
2681   if ( fDebugStreamer ) {
2682     delete fDebugStreamer;
2683     fDebugStreamer = 0x0;
2684   }
2685
2686   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
2687                ,fNumberTrack
2688                ,fNumberUsedCh[0]
2689                ,fNumberUsedCh[1]
2690                ,fNumberUsedPh[0]
2691                ,fNumberUsedPh[1]));
2692
2693   TDirectory *backup = gDirectory;
2694   TString option;
2695
2696   if ( append )
2697     option = "update";
2698   else
2699     option = "recreate";
2700
2701   TFile f(filename,option.Data());
2702
2703   TStopwatch stopwatch;
2704   stopwatch.Start();
2705   if(fVector2d) {
2706     f.WriteTObject(fCalibraVector);
2707   }
2708
2709   if (fCH2dOn ) {
2710     if (fHisto2d) {
2711       f.WriteTObject(fCH2d);
2712     }
2713   }
2714   if (fPH2dOn ) {
2715     if (fHisto2d) {
2716       f.WriteTObject(fPH2d);
2717     }
2718   }
2719   if (fPRF2dOn) {
2720     if (fHisto2d) {
2721         f.WriteTObject(fPRF2d);
2722     }
2723   }
2724   if(fLinearFitterOn){
2725     AnalyseLinearFitter();
2726     f.WriteTObject(fLinearVdriftFit);
2727   }
2728
2729   f.Close();
2730
2731   if ( backup ) backup->cd();
2732
2733   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
2734                ,stopwatch.RealTime(),stopwatch.CpuTime()));
2735 }
2736 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2737 // Stats stuff
2738 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
2739 //___________________________________________probe the histos__________________________________________________
2740 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
2741 {
2742   //
2743   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
2744   // debug mode with 2 for TH2I and 3 for TProfile2D
2745   // It gives a pointer to a Double_t[7] with the info following...
2746   // [0] : number of calibration groups with entries
2747   // [1] : minimal number of entries found
2748   // [2] : calibration group number of the min
2749   // [3] : maximal number of entries found
2750   // [4] : calibration group number of the max
2751   // [5] : mean number of entries found
2752   // [6] : mean relative error
2753   //
2754
2755   Double_t *info = new Double_t[7];
2756
2757   // Number of Xbins (detectors or groups of pads)
2758   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
2759   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
2760
2761   // Initialise
2762   Double_t nbwe = 0; //number of calibration groups with entries
2763   Double_t minentries = 0; //minimal number of entries found
2764   Double_t maxentries = 0; //maximal number of entries found
2765   Double_t placemin = 0; //calibration group number of the min
2766   Double_t placemax = -1; //calibration group number of the max
2767   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
2768   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
2769
2770   Double_t counter = 0;
2771
2772   //Debug
2773   TH1F *nbEntries = 0x0;//distribution of the number of entries
2774   TH1F *nbEntriesPerGroup = 0x0;//Number of entries per group
2775   TProfile *nbEntriesPerSp = 0x0;//Number of entries for one supermodule
2776
2777   // Beginning of the loop over the calibration groups
2778   for (Int_t idect = 0; idect < nbins; idect++) {
2779
2780     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
2781     projch->SetDirectory(0);
2782
2783     // Number of entries for this calibration group
2784     Double_t nentries = 0.0;
2785     if((i%2) == 0){
2786       for (Int_t k = 0; k < nxbins; k++) {
2787         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
2788       }
2789     }
2790     else{
2791       for (Int_t k = 0; k < nxbins; k++) {
2792         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
2793         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
2794           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
2795           counter++;
2796         }
2797       }
2798     }
2799
2800     //Debug
2801     if(i > 1){
2802       if((!((Bool_t)nbEntries)) && (nentries > 0)){
2803         nbEntries = new TH1F("Number of entries","Number of entries"
2804                                ,100,(Int_t)nentries/2,nentries*2);
2805         nbEntries->SetDirectory(0);
2806         nbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
2807                                ,nbins,0,nbins);
2808         nbEntriesPerGroup->SetDirectory(0);
2809         nbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
2810                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
2811         nbEntriesPerSp->SetDirectory(0);
2812       }
2813       if(nbEntries){
2814         if(nentries > 0) nbEntries->Fill(nentries);
2815         nbEntriesPerGroup->Fill(idect+0.5,nentries);
2816         nbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
2817       }
2818     }
2819
2820     //min amd max
2821     if(nentries > maxentries){
2822       maxentries = nentries;
2823       placemax = idect;
2824     }
2825     if(idect == 0) {
2826       minentries = nentries;
2827     }
2828     if(nentries < minentries){
2829       minentries = nentries;
2830       placemin = idect;
2831     }
2832     //nbwe
2833     if(nentries > 0) {
2834       nbwe++;
2835       meanstats += nentries;
2836     }
2837   }//calibration groups loop
2838
2839   if(nbwe > 0) meanstats /= nbwe;
2840   if(counter > 0) meanrelativerror /= counter;
2841
2842   AliInfo(Form("There are %f calibration groups with entries",nbwe));
2843   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
2844   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
2845   AliInfo(Form("The mean number of entries is %f",meanstats));
2846   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
2847
2848   info[0] = nbwe;
2849   info[1] = minentries;
2850   info[2] = placemin;
2851   info[3] = maxentries;
2852   info[4] = placemax;
2853   info[5] = meanstats;
2854   info[6] = meanrelativerror;
2855
2856   if(i > 1){
2857     gStyle->SetPalette(1);
2858     gStyle->SetOptStat(1111);
2859     gStyle->SetPadBorderMode(0);
2860     gStyle->SetCanvasColor(10);
2861     gStyle->SetPadLeftMargin(0.13);
2862     gStyle->SetPadRightMargin(0.01);
2863     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
2864     stat->Divide(2,1);
2865     stat->cd(1);
2866     nbEntries->Draw("");
2867     stat->cd(2);
2868     nbEntriesPerSp->SetStats(0);
2869     nbEntriesPerSp->Draw("");
2870     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
2871     stat1->cd();
2872     nbEntriesPerGroup->SetStats(0);
2873     nbEntriesPerGroup->Draw("");
2874   }
2875
2876   return info;
2877
2878 }
2879 //____________________________________________________________________________
2880 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
2881 {
2882   //
2883   // Return a Int_t[4] with:
2884   // 0 Mean number of entries
2885   // 1 median of number of entries
2886   // 2 rms of number of entries
2887   // 3 number of group with entries
2888   //
2889
2890   Double_t *stat      = new Double_t[4];
2891   stat[3]             = 0.0;
2892
2893   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
2894   Double_t *weight    = new Double_t[nbofgroups];
2895   Int_t    *nonul     = new Int_t[nbofgroups];
2896
2897   for(Int_t k = 0; k < nbofgroups; k++){
2898     if(fEntriesCH[k] > 0) {
2899       weight[k] = 1.0;
2900       nonul[(Int_t)stat[3]] = fEntriesCH[k];
2901       stat[3]++;
2902     }
2903     else weight[k] = 0.0;
2904   }
2905   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
2906   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
2907   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2908
2909   return stat;
2910
2911 }
2912 //____________________________________________________________________________
2913 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
2914 {
2915   //
2916   // Return a Int_t[4] with:
2917   // 0 Mean number of entries
2918   // 1 median of number of entries
2919   // 2 rms of number of entries
2920   // 3 number of group with entries
2921   //
2922
2923   Double_t *stat      = new Double_t[4];
2924   stat[3]             = 0.0;
2925
2926   Int_t    nbofgroups = 540;
2927   Double_t *weight    = new Double_t[nbofgroups];
2928   Int_t    *nonul     = new Int_t[nbofgroups];
2929
2930   for(Int_t k = 0; k < nbofgroups; k++){
2931     if(fEntriesLinearFitter[k] > 0) {
2932       weight[k] = 1.0;
2933       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
2934       stat[3]++;
2935     }
2936     else weight[k] = 0.0;
2937   }
2938   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
2939   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
2940   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
2941
2942   return stat;
2943
2944 }
2945 //////////////////////////////////////////////////////////////////////////////////////
2946 // Create Histos
2947 //////////////////////////////////////////////////////////////////////////////////////
2948 //_____________________________________________________________________________
2949 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
2950 {
2951   //
2952   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
2953   // If fNgroupprf is zero then no binning in tnp
2954   //
2955
2956   TString name("Nz");
2957   name += fCalibraMode->GetNz(2);
2958   name += "Nrphi";
2959   name += fCalibraMode->GetNrphi(2);
2960   name += "Ngp";
2961   name += fNgroupprf;
2962
2963   if(fNgroupprf != 0){
2964
2965     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2966                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
2967     fPRF2d->SetYTitle("Det/pad groups");
2968     fPRF2d->SetXTitle("Position x/W [pad width units]");
2969     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2970     fPRF2d->SetStats(0);
2971   }
2972   else{
2973     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
2974                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
2975     fPRF2d->SetYTitle("Det/pad groups");
2976     fPRF2d->SetXTitle("Position x/W [pad width units]");
2977     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
2978     fPRF2d->SetStats(0);
2979   }
2980
2981 }
2982
2983 //_____________________________________________________________________________
2984 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
2985 {
2986   //
2987   // Create the 2D histos
2988   //
2989
2990   TString name("Nz");
2991   name += fCalibraMode->GetNz(1);
2992   name += "Nrphi";
2993   name += fCalibraMode->GetNrphi(1);
2994
2995   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
2996                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
2997                          ,nn,0,nn);
2998   fPH2d->SetYTitle("Det/pad groups");
2999   fPH2d->SetXTitle("time [#mus]");
3000   fPH2d->SetZTitle("<PH> [a.u.]");
3001   fPH2d->SetStats(0);
3002
3003 }
3004 //_____________________________________________________________________________
3005 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
3006 {
3007   //
3008   // Create the 2D histos
3009   //
3010
3011   TString name("Nz");
3012   name += fCalibraMode->GetNz(0);
3013   name += "Nrphi";
3014   name += fCalibraMode->GetNrphi(0);
3015
3016   fCH2d = new TH2I("CH2d",(const Char_t *) name
3017                    ,fNumberBinCharge,0,300,nn,0,nn);
3018   fCH2d->SetYTitle("Det/pad groups");
3019   fCH2d->SetXTitle("charge deposit [a.u]");
3020   fCH2d->SetZTitle("counts");
3021   fCH2d->SetStats(0);
3022   fCH2d->Sumw2();
3023
3024 }
3025 //////////////////////////////////////////////////////////////////////////////////
3026 // Set relative scale
3027 /////////////////////////////////////////////////////////////////////////////////
3028 //_____________________________________________________________________________
3029 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
3030 {
3031   //
3032   // Set the factor that will divide the deposited charge
3033   // to fit in the histo range [0,300]
3034   //
3035
3036   if (RelativeScale > 0.0) {
3037     fRelativeScale = RelativeScale;
3038   }
3039   else {
3040     AliInfo("RelativeScale must be strict positif!");
3041   }
3042
3043 }
3044 //////////////////////////////////////////////////////////////////////////////////
3045 // Quick way to fill a histo
3046 //////////////////////////////////////////////////////////////////////////////////
3047 //_____________________________________________________________________
3048 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3049 {
3050   //
3051   // FillCH2d: Marian style
3052   //
3053
3054   //skip simply the value out of range
3055   if((y>=300.0) || (y<0.0)) return;
3056
3057   //Calcul the y place
3058   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3059   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3060
3061   //Fill
3062   fCH2d->GetArray()[place]++;
3063
3064 }
3065
3066 //////////////////////////////////////////////////////////////////////////////////
3067 // Geometrical functions
3068 ///////////////////////////////////////////////////////////////////////////////////
3069 //_____________________________________________________________________________
3070 Int_t AliTRDCalibraFillHisto::GetLayer(Int_t d) const
3071 {
3072   //
3073   // Reconstruct the layer number from the detector number
3074   //
3075
3076   return ((Int_t) (d % 6));
3077
3078 }
3079
3080 //_____________________________________________________________________________
3081 Int_t AliTRDCalibraFillHisto::GetStack(Int_t d) const
3082 {
3083   //
3084   // Reconstruct the stack number from the detector number
3085   //
3086   const Int_t kNlayer = 6;
3087
3088   return ((Int_t) (d % 30) / kNlayer);
3089
3090 }
3091
3092 //_____________________________________________________________________________
3093 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
3094 {
3095   //
3096   // Reconstruct the sector number from the detector number
3097   //
3098   Int_t fg = 30;
3099
3100   return ((Int_t) (d / fg));
3101
3102 }
3103 ///////////////////////////////////////////////////////////////////////////////////
3104 // Getter functions for DAQ of the CH2d and the PH2d
3105 //////////////////////////////////////////////////////////////////////////////////
3106 //_____________________________________________________________________
3107 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
3108 {
3109     //
3110     // return pointer to fPH2d TProfile2D
3111     // create a new TProfile2D if it doesn't exist allready
3112     //
3113     if ( fPH2d )
3114         return fPH2d;
3115
3116     // Some parameters
3117     fTimeMax = nbtimebin;
3118     fSf      = samplefrequency;
3119
3120     CreatePH2d(540);
3121
3122     return fPH2d;
3123 }
3124 //_____________________________________________________________________
3125 TH2I* AliTRDCalibraFillHisto::GetCH2d()
3126 {
3127     //
3128     // return pointer to fCH2d TH2I
3129     // create a new TH2I if it doesn't exist allready
3130     //
3131     if ( fCH2d )
3132         return fCH2d;
3133
3134     CreateCH2d(540);
3135
3136     return fCH2d;
3137 }
3138 ////////////////////////////////////////////////////////////////////////////////////////////
3139 // Drift velocity calibration
3140 ///////////////////////////////////////////////////////////////////////////////////////////
3141 //_____________________________________________________________________
3142 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3143 {
3144     //
3145     // return pointer to TLinearFitter Calibration
3146     // if force is true create a new TLinearFitter if it doesn't exist allready
3147     //
3148
3149   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3150     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3151   }
3152
3153   // if we are forced and TLinearFitter doesn't yet exist create it
3154
3155   // new TLinearFitter
3156   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3157   fLinearFitterArray.AddAt(linearfitter,detector);
3158   return linearfitter;
3159 }
3160
3161 //____________________________________________________________________________
3162 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3163 {
3164   //
3165   // Analyse array of linear fitter because can not be written
3166   // Store two arrays: one with the param the other one with the error param + number of entries
3167   //
3168
3169   for(Int_t k = 0; k < 540; k++){
3170     TLinearFitter *linearfitter = GetLinearFitter(k);
3171     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3172       TVectorD  *par  = new TVectorD(2);
3173       TVectorD   pare = TVectorD(2);
3174       TVectorD  *parE = new TVectorD(3);
3175       linearfitter->Eval();
3176       linearfitter->GetParameters(*par);
3177       linearfitter->GetErrors(pare);
3178       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3179       (*parE)[0] = pare[0]*ppointError;
3180       (*parE)[1] = pare[1]*ppointError;
3181       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3182       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3183       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3184     }
3185   }
3186 }
3187