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