TRD/AliTRDCalibraFillHisto.cxx

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