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