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 "AliTRDRawStreamV2.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->GetPadCol();
 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(AliTRDRawStreamV2 *rawStream, Bool_t nocheck)
1208 {
1209   //
1210   // Event Processing loop - AliTRDRawStreamV2
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       //printf("detector %d, nbtimebin %d, iTimeBin %d\n",fDetectorPreviousTrack,nbtimebin,iTimeBin);
1257
1258       Int_t fin     = TMath::Min(fTimeMax,(iTimeBin+3));
1259       Int_t n       = 0;
1260       for(Int_t itime = iTimeBin; itime < fin; itime++){
1261         // should extract baseline here!
1262         if(signal[n]>13) phvalue[itime] = signal[n];
1263         //printf("signal is %d for %d\n",signal[n],n);
1264         n++;
1265       }
1266     }
1267
1268     // fill the last one
1269     if(fDetectorPreviousTrack != -1){
1270       if(TMath::Mean(fTimeMax,phvalue)>20.0){
1271         withInput = 2;
1272         for(Int_t k = 0; k < fTimeMax; k++){
1273           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],fTimeMax);
1274           phvalue[k] = 10;
1275           //row[k]     = -1;
1276           //col[k]     = -1;
1277         }
1278       }
1279     }
1280
1281   }
1282   else{
1283
1284     while (rawStream->Next()) {
1285
1286       Int_t idetector = rawStream->GetDet();                            //  current detector
1287       if((fDetectorPreviousTrack != idetector) && (fDetectorPreviousTrack != -1)){
1288         if(TMath::Mean(nbtimebin,phvalue)>20.0){
1289           withInput = 2;
1290           for(Int_t k = 0; k < nbtimebin; k++){
1291             UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1292             phvalue[k] = 10;
1293             //row[k]     = -1;
1294             //col[k]     = -1;
1295           }
1296         }
1297       }
1298       fDetectorPreviousTrack = idetector;
1299       nbtimebin         = rawStream->GetNumberOfTimeBins();              //  number of time bins read from data
1300       Int_t iTimeBin    = rawStream->GetTimeBin();                       //  current time bin
1301       //row[iTimeBin]   = rawStream->GetRow();                           //  current row
1302       //col[iTimeBin]   = rawStream->GetCol();                           //  current col
1303       Int_t *signal     = rawStream->GetSignals();                       //  current ADC signal
1304       //printf("detector %d, nbtimebin %d, iTimeBin %d\n",fDetectorPreviousTrack,nbtimebin,iTimeBin);
1305
1306       Int_t fin     = TMath::Min(nbtimebin,(iTimeBin+3));
1307       Int_t n       = 0;
1308       for(Int_t itime = iTimeBin; itime < fin; itime++){
1309         // should extract baseline here!
1310         if(signal[n]>13) phvalue[itime] = signal[n];
1311         //printf("signal is %d for %d\n",signal[n],n);
1312         n++;
1313       }
1314     }
1315
1316     // fill the last one
1317     if(fDetectorPreviousTrack != -1){
1318       if(TMath::Mean(nbtimebin,phvalue)>20.0){
1319         withInput = 2;
1320         for(Int_t k = 0; k < nbtimebin; k++){
1321           UpdateDAQ(fDetectorPreviousTrack,0,0,k,phvalue[k],nbtimebin);
1322           phvalue[k] = 10;
1323           //row[k]     = -1;
1324           //col[k]     = -1;
1325         }
1326       }
1327     }
1328   }
1329
1330   return withInput;
1331
1332 }
1333 //_____________________________________________________________________
1334 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(AliRawReader *rawReader, Bool_t nocheck)
1335 {
1336   //
1337   //  Event processing loop - AliRawReader
1338   //
1339
1340
1341   AliTRDRawStreamV2 rawStream(rawReader);
1342
1343   rawReader->Select("TRD");
1344
1345   return ProcessEventDAQ(&rawStream, nocheck);
1346 }
1347 //_________________________________________________________________________
1348 Int_t AliTRDCalibraFillHisto::ProcessEventDAQ(
1349 #ifdef ALI_DATE
1350                                                eventHeaderStruct *event,
1351                                                Bool_t nocheck
1352 #else
1353                                                eventHeaderStruct* /*event*/,
1354                                                Bool_t /*nocheck*/
1355
1356 #endif
1357                                    )
1358 {
1359   //
1360   //  process date event
1361   //
1362 #ifdef ALI_DATE
1363     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1364     Int_t result=ProcessEventDAQ(rawReader, nocheck);
1365     delete rawReader;
1366     return result;
1367 #else
1368     Fatal("AliTRDCalibraFillHisto", "this class was compiled without DATE");
1369     return 0;
1370 #endif
1371
1372 }
1373 //____________Online trackling in AliTRDtrigger________________________________
1374 Bool_t AliTRDCalibraFillHisto::UpdateDAQ(Int_t det, Int_t /*row*/, Int_t /*col*/, Int_t timebin, Int_t signal, Int_t nbtimebins)
1375 {
1376   //
1377   // For the DAQ
1378   // Fill a simple average pulse height
1379   //
1380
1381   // Localisation of the Xbins involved
1382   //LocalisationDetectorXbins(det);
1383
1384   // Row  and position in the pad groups
1385   //Int_t posr = 0;
1386   //if (fCalibraMode->GetNnZ(1) != 0) {
1387   //  posr = (Int_t) row / fCalibraMode->GetNnZ(1);
1388   //}
1389
1390   // Col of cluster and position in the pad groups
1391   //Int_t posc = 0;
1392   //if (fCalibraMode->GetNnRphi(1) != 0) {
1393   //  posc = (Int_t) col / fCalibraMode->GetNnRphi(1);
1394   //}
1395
1396   //fPH2d->Fill((Float_t) timebin/fSf,(fCalibraMode->GetXbins(1)+posc*fCalibraMode->GetNfragZ(1)+posr)+0.5,(Float_t) signal);
1397
1398   ((TProfile2D *)GetPH2d(nbtimebins,fSf))->Fill((Float_t) timebin/fSf,det+0.5,(Float_t) signal);
1399
1400   return kTRUE;
1401
1402 }
1403 //____________Write_____________________________________________________
1404 //_____________________________________________________________________
1405 void AliTRDCalibraFillHisto::Write2d(const Char_t *filename, Bool_t append)
1406 {
1407   //
1408   //  Write infos to file
1409   //
1410
1411   //For debugging
1412   if ( fDebugStreamer ) {
1413     delete fDebugStreamer;
1414     fDebugStreamer = 0x0;
1415   }
1416
1417   AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
1418                ,fNumberTrack
1419                ,fNumberUsedCh[0]
1420                ,fNumberUsedCh[1]
1421                ,fNumberUsedPh[0]
1422                ,fNumberUsedPh[1]));
1423
1424   TDirectory *backup = gDirectory;
1425   TString option;
1426
1427   if ( append )
1428     option = "update";
1429   else
1430     option = "recreate";
1431
1432   TFile f(filename,option.Data());
1433
1434   TStopwatch stopwatch;
1435   stopwatch.Start();
1436   if(fVector2d) {
1437     f.WriteTObject(fCalibraVector);
1438   }
1439
1440   if (fCH2dOn ) {
1441     if (fHisto2d) {
1442       f.WriteTObject(fCH2d);
1443     }
1444   }
1445   if (fPH2dOn ) {
1446     if (fHisto2d) {
1447       f.WriteTObject(fPH2d);
1448     }
1449   }
1450   if (fPRF2dOn) {
1451     if (fHisto2d) {
1452         f.WriteTObject(fPRF2d);
1453     }
1454   }
1455   if(fLinearFitterOn){
1456     AnalyseLinearFitter();
1457     f.WriteTObject(fLinearVdriftFit);
1458   }
1459
1460   f.Close();
1461
1462   if ( backup ) backup->cd();
1463
1464   AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
1465                ,stopwatch.RealTime(),stopwatch.CpuTime()));
1466 }
1467 //___________________________________________probe the histos__________________________________________________
1468 Double_t *AliTRDCalibraFillHisto::StatH(TH2 *h, Int_t i)
1469 {
1470   //
1471   // Check the number of stats in h, 0 is TH2I 1 is TProfile2D
1472   // debug mode with 2 for TH2I and 3 for TProfile2D
1473   // It gives a pointer to a Double_t[7] with the info following...
1474   // [0] : number of calibration groups with entries
1475   // [1] : minimal number of entries found
1476   // [2] : calibration group number of the min
1477   // [3] : maximal number of entries found
1478   // [4] : calibration group number of the max
1479   // [5] : mean number of entries found
1480   // [6] : mean relative error
1481   //
1482
1483   Double_t *info = new Double_t[7];
1484
1485   // Number of Xbins (detectors or groups of pads)
1486   Int_t    nbins   = h->GetNbinsY(); //number of calibration groups
1487   Int_t    nxbins  = h->GetNbinsX(); //number of bins per histo
1488
1489   // Initialise
1490   Double_t nbwe = 0; //number of calibration groups with entries
1491   Double_t minentries = 0; //minimal number of entries found
1492   Double_t maxentries = 0; //maximal number of entries found
1493   Double_t placemin = 0; //calibration group number of the min
1494   Double_t placemax = -1; //calibration group number of the max
1495   Double_t meanstats = 0.0; //mean number of entries over the calibration group with at least ome entry
1496   Double_t meanrelativerror = 0.0; //mean relativ error in the TProfile2D
1497
1498   Double_t counter = 0;
1499
1500   //Debug
1501   TH1F *NbEntries = 0x0;//distribution of the number of entries
1502   TH1F *NbEntriesPerGroup = 0x0;//Number of entries per group
1503   TProfile *NbEntriesPerSp = 0x0;//Number of entries for one supermodule
1504
1505   // Beginning of the loop over the calibration groups
1506   for (Int_t idect = 0; idect < nbins; idect++) {
1507
1508     TH1I *projch = (TH1I *) h->ProjectionX("projch",idect+1,idect+1,(Option_t *)"e");
1509     projch->SetDirectory(0);
1510
1511     // Number of entries for this calibration group
1512     Double_t nentries = 0.0;
1513     if((i%2) == 0){
1514       for (Int_t k = 0; k < nxbins; k++) {
1515         nentries += h->GetBinContent(h->GetBin(k+1,idect+1));
1516       }
1517     }
1518     else{
1519       for (Int_t k = 0; k < nxbins; k++) {
1520         nentries += ((TProfile2D *)h)->GetBinEntries(h->GetBin(k+1,idect+1));
1521         if(h->GetBinContent(h->GetBin(k+1,idect+1)) != 0) {
1522           meanrelativerror += (h->GetBinError(h->GetBin(k+1,idect+1))/(TMath::Abs(h->GetBinContent(h->GetBin(k+1,idect+1)))));
1523           counter++;
1524         }
1525       }
1526     }
1527
1528     //Debug
1529     if(i > 1){
1530       if((!((Bool_t)NbEntries)) && (nentries > 0)){
1531         NbEntries = new TH1F("Number of entries","Number of entries"
1532                                ,100,(Int_t)nentries/2,nentries*2);
1533         NbEntries->SetDirectory(0);
1534         NbEntriesPerGroup = new TH1F("Number of entries per group","Number of entries per group"
1535                                ,nbins,0,nbins);
1536         NbEntriesPerGroup->SetDirectory(0);
1537         NbEntriesPerSp = new TProfile("Number of entries per supermodule","Number of entries per supermodule"
1538                                ,(Int_t)(nbins/18),0,(Int_t)(nbins/18));
1539         NbEntriesPerSp->SetDirectory(0);
1540       }
1541       if(NbEntries){
1542         if(nentries > 0) NbEntries->Fill(nentries);
1543         NbEntriesPerGroup->Fill(idect+0.5,nentries);
1544         NbEntriesPerSp->Fill((idect%((Int_t)(nbins/18)))+0.5,nentries);
1545       }
1546     }
1547
1548     //min amd max
1549     if(nentries > maxentries){
1550       maxentries = nentries;
1551       placemax = idect;
1552     }
1553     if(idect == 0) {
1554       minentries = nentries;
1555     }
1556     if(nentries < minentries){
1557       minentries = nentries;
1558       placemin = idect;
1559     }
1560     //nbwe
1561     if(nentries > 0) {
1562       nbwe++;
1563       meanstats += nentries;
1564     }
1565   }//calibration groups loop
1566
1567   if(nbwe > 0) meanstats /= nbwe;
1568   if(counter > 0) meanrelativerror /= counter;
1569
1570   AliInfo(Form("There are %f calibration groups with entries",nbwe));
1571   AliInfo(Form("The minimum number of entries is %f for the group %f",minentries,placemin));
1572   AliInfo(Form("The maximum number of entries is %f for the group %f",maxentries,placemax));
1573   AliInfo(Form("The mean number of entries is %f",meanstats));
1574   if((i%2) == 1) AliInfo(Form("The mean relative error is %f",meanrelativerror));
1575
1576   info[0] = nbwe;
1577   info[1] = minentries;
1578   info[2] = placemin;
1579   info[3] = maxentries;
1580   info[4] = placemax;
1581   info[5] = meanstats;
1582   info[6] = meanrelativerror;
1583
1584   if(i > 1){
1585     gStyle->SetPalette(1);
1586     gStyle->SetOptStat(1111);
1587     gStyle->SetPadBorderMode(0);
1588     gStyle->SetCanvasColor(10);
1589     gStyle->SetPadLeftMargin(0.13);
1590     gStyle->SetPadRightMargin(0.01);
1591     TCanvas *stat = new TCanvas("stat","",50,50,600,800);
1592     stat->Divide(2,1);
1593     stat->cd(1);
1594     NbEntries->Draw("");
1595     stat->cd(2);
1596     NbEntriesPerSp->SetStats(0);
1597     NbEntriesPerSp->Draw("");
1598     TCanvas *stat1 = new TCanvas("stat1","",50,50,600,800);
1599     stat1->cd();
1600     NbEntriesPerGroup->SetStats(0);
1601     NbEntriesPerGroup->Draw("");
1602   }
1603
1604   return info;
1605
1606 }
1607 //____________________________________________________________________________
1608 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberCH()
1609 {
1610   //
1611   // Return a Int_t[4] with:
1612   // 0 Mean number of entries
1613   // 1 median of number of entries
1614   // 2 rms of number of entries
1615   // 3 number of group with entries
1616   //
1617
1618   Double_t *stat      = new Double_t[4];
1619   stat[3]             = 0.0;
1620
1621   Int_t    nbofgroups = CalculateTotalNumberOfBins(0);
1622   Double_t *weight    = new Double_t[nbofgroups];
1623   Int_t    *nonul     = new Int_t[nbofgroups];
1624
1625   for(Int_t k = 0; k < nbofgroups; k++){
1626     if(fEntriesCH[k] > 0) {
1627       weight[k] = 1.0;
1628       nonul[(Int_t)stat[3]] = fEntriesCH[k];
1629       stat[3]++;
1630     }
1631     else weight[k] = 0.0;
1632   }
1633   stat[0]          = TMath::Mean(nbofgroups,fEntriesCH,weight);
1634   stat[1]          = TMath::Median(nbofgroups,fEntriesCH,weight);
1635   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
1636
1637   return stat;
1638
1639 }
1640 //____________________________________________________________________________
1641 Double_t *AliTRDCalibraFillHisto::GetMeanMedianRMSNumberLinearFitter() const
1642 {
1643   //
1644   // Return a Int_t[4] with:
1645   // 0 Mean number of entries
1646   // 1 median of number of entries
1647   // 2 rms of number of entries
1648   // 3 number of group with entries
1649   //
1650
1651   Double_t *stat      = new Double_t[4];
1652   stat[3]             = 0.0;
1653
1654   Int_t    nbofgroups = 540;
1655   Double_t *weight    = new Double_t[nbofgroups];
1656   Int_t    *nonul     = new Int_t[nbofgroups];
1657
1658   for(Int_t k = 0; k < nbofgroups; k++){
1659     if(fEntriesLinearFitter[k] > 0) {
1660       weight[k] = 1.0;
1661       nonul[(Int_t) stat[3]] = fEntriesLinearFitter[k];
1662       stat[3]++;
1663     }
1664     else weight[k] = 0.0;
1665   }
1666   stat[0]          = TMath::Mean(nbofgroups,fEntriesLinearFitter,weight);
1667   stat[1]          = TMath::Median(nbofgroups,fEntriesLinearFitter,weight);
1668   stat[2]          = TMath::RMS((Int_t)stat[3],nonul);
1669
1670   return stat;
1671
1672 }
1673 //_____________________________________________________________________________
1674 void AliTRDCalibraFillHisto::SetNumberGroupsPRF(Short_t numberGroupsPRF)
1675 {
1676   //
1677   // Should be between 0 and 6
1678   //
1679
1680   if ((numberGroupsPRF < 0) || (numberGroupsPRF > 6)) {
1681     AliInfo("The number of groups must be between 0 and 6!");
1682   }
1683   else {
1684     fNgroupprf = numberGroupsPRF;
1685   }
1686
1687 }
1688 //_____________________________________________________________________________
1689 void AliTRDCalibraFillHisto::SetRelativeScale(Float_t RelativeScale)
1690 {
1691   //
1692   // Set the factor that will divide the deposited charge
1693   // to fit in the histo range [0,300]
1694   //
1695
1696   if (RelativeScale > 0.0) {
1697     fRelativeScale = RelativeScale;
1698   }
1699   else {
1700     AliInfo("RelativeScale must be strict positif!");
1701   }
1702
1703 }
1704
1705 //_____________________________________________________________________________
1706 void AliTRDCalibraFillHisto::SetNz(Int_t i, Short_t Nz)
1707 {
1708   //
1709   // Set the mode of calibration group in the z direction for the parameter i
1710   //
1711
1712   if ((Nz >= 0) &&
1713       (Nz <  5)) {
1714     fCalibraMode->SetNz(i, Nz);
1715   }
1716   else {
1717     AliInfo("You have to choose between 0 and 4");
1718   }
1719
1720 }
1721
1722 //_____________________________________________________________________________
1723 void AliTRDCalibraFillHisto::SetNrphi(Int_t i, Short_t Nrphi)
1724 {
1725   //
1726   // Set the mode of calibration group in the rphi direction for the parameter i
1727   //
1728
1729   if ((Nrphi >= 0) &&
1730       (Nrphi <  7)) {
1731     fCalibraMode->SetNrphi(i ,Nrphi);
1732   }
1733   else {
1734     AliInfo("You have to choose between 0 and 6");
1735   }
1736
1737 }
1738 //____________Protected Functions______________________________________________
1739 //____________Create the 2D histo to be filled online__________________________
1740 //
1741 //_____________________________________________________________________________
1742 void AliTRDCalibraFillHisto::CreatePRF2d(Int_t nn)
1743 {
1744   //
1745   // Create the 2D histos: here we have 2*fNgroupprf bins in tnp of 0.2 amplitude each
1746   // If fNgroupprf is zero then no binning in tnp
1747   //
1748
1749   TString name("Nz");
1750   name += fCalibraMode->GetNz(2);
1751   name += "Nrphi";
1752   name += fCalibraMode->GetNrphi(2);
1753   name += "Ngp";
1754   name += fNgroupprf;
1755
1756   if(fNgroupprf != 0){
1757
1758     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1759                             ,2*fNgroupprf*fNumberBinPRF,-3.0*fNgroupprf,3.0*fNgroupprf,nn,0,nn);
1760     fPRF2d->SetYTitle("Det/pad groups");
1761     fPRF2d->SetXTitle("Position x/W [pad width units]");
1762     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1763     fPRF2d->SetStats(0);
1764   }
1765   else{
1766     fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
1767                             ,fNumberBinPRF,-1.5,1.5,nn,0,nn);
1768     fPRF2d->SetYTitle("Det/pad groups");
1769     fPRF2d->SetXTitle("Position x/W [pad width units]");
1770     fPRF2d->SetZTitle("Q_{i}/Q_{total}");
1771     fPRF2d->SetStats(0);
1772   }
1773
1774 }
1775
1776 //_____________________________________________________________________________
1777 void AliTRDCalibraFillHisto::CreatePH2d(Int_t nn)
1778 {
1779   //
1780   // Create the 2D histos
1781   //
1782
1783   TString name("Nz");
1784   name += fCalibraMode->GetNz(1);
1785   name += "Nrphi";
1786   name += fCalibraMode->GetNrphi(1);
1787
1788   fPH2d = new TProfile2D("PH2d",(const Char_t *) name
1789                          ,fTimeMax,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf
1790                          ,nn,0,nn);
1791   fPH2d->SetYTitle("Det/pad groups");
1792   fPH2d->SetXTitle("time [#mus]");
1793   fPH2d->SetZTitle("<PH> [a.u.]");
1794   fPH2d->SetStats(0);
1795
1796 }
1797 //_____________________________________________________________________________
1798 void AliTRDCalibraFillHisto::CreateCH2d(Int_t nn)
1799 {
1800   //
1801   // Create the 2D histos
1802   //
1803
1804   TString name("Nz");
1805   name += fCalibraMode->GetNz(0);
1806   name += "Nrphi";
1807   name += fCalibraMode->GetNrphi(0);
1808
1809   fCH2d = new TH2I("CH2d",(const Char_t *) name
1810                    ,fNumberBinCharge,0,300,nn,0,nn);
1811   fCH2d->SetYTitle("Det/pad groups");
1812   fCH2d->SetXTitle("charge deposit [a.u]");
1813   fCH2d->SetZTitle("counts");
1814   fCH2d->SetStats(0);
1815   fCH2d->Sumw2();
1816
1817 }
1818
1819 //____________Offine tracking in the AliTRDtracker_____________________________
1820 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackCH()
1821 {
1822   //
1823   // For the offline tracking or mcm tracklets
1824   // This function will be called in the functions UpdateHistogram...
1825   // to fill the info of a track for the relativ gain calibration
1826   //
1827
1828   Int_t nb            =  0;   // Nombre de zones traversees
1829   Int_t fd            = -1;   // Premiere zone non nulle
1830   Float_t totalcharge = 0.0;  // Total charge for the supermodule histo
1831
1832
1833
1834
1835   // See if the track goes through different zones
1836   for (Int_t k = 0; k < fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0); k++) {
1837     if (fAmpTotal[k] > 0.0) {
1838       totalcharge += fAmpTotal[k];
1839       nb++;
1840       if (nb == 1) {
1841         fd = k;
1842       }
1843     }
1844   }
1845
1846
1847   switch (nb)
1848     {
1849     case 1:
1850       fNumberUsedCh[0]++;
1851       fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1852       if (fHisto2d) {
1853         FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1854         //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1855       }
1856       if (fVector2d) {
1857         fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1858       }
1859       break;
1860     case 2:
1861       if ((fAmpTotal[fd]   > 0.0) &&
1862           (fAmpTotal[fd+1] > 0.0)) {
1863         // One of the two very big
1864         if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
1865           if (fHisto2d) {
1866             FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1867             //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1868           }
1869           if (fVector2d) {
1870             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1871           }
1872           fNumberUsedCh[1]++;
1873           fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1874         }
1875         if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
1876           if (fHisto2d) {
1877             FillCH2d(fCalibraMode->GetXbins(0)+fd+1,fAmpTotal[fd+1]/fRelativeScale);
1878             //fCH2d->Fill(fAmpTotal[fd+1]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+1.5);
1879           }
1880           if (fVector2d) {
1881             fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+1,fAmpTotal[fd+1]/fRelativeScale);
1882           }
1883           fNumberUsedCh[1]++;
1884           fEntriesCH[fCalibraMode->GetXbins(0)+fd+1]++;
1885         }
1886       }
1887       if (fCalibraMode->GetNfragZ(0) > 1) {
1888         if (fAmpTotal[fd] > 0.0) {
1889           if ((fd+fCalibraMode->GetNfragZ(0)) < (fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))) {
1890             if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > 0.0) {
1891               // One of the two very big
1892               if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]) {
1893                 if (fHisto2d) {
1894                   FillCH2d(fCalibraMode->GetXbins(0)+fd,fAmpTotal[fd]/fRelativeScale);
1895                   //fCH2d->Fill(fAmpTotal[fd]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+0.5);
1896                 }
1897                 if (fVector2d) {
1898                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd,fAmpTotal[fd]/fRelativeScale);
1899                 }
1900                 fNumberUsedCh[1]++;
1901                 fEntriesCH[fCalibraMode->GetXbins(0)+fd]++;
1902               }
1903               if (fAmpTotal[fd+fCalibraMode->GetNfragZ(0)] > fProcent*fAmpTotal[fd]) {
1904                 if (fHisto2d) {
1905                   FillCH2d(fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1906                   //fCH2d->Fill(fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale,fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)+0.5);
1907                 }
1908                 fNumberUsedCh[1]++;
1909                 fEntriesCH[fCalibraMode->GetXbins(0)+fd+fCalibraMode->GetNfragZ(0)]++;
1910                 if (fVector2d) {
1911                   fCalibraVector->UpdateVectorCH(fDetectorPreviousTrack,fd+fCalibraMode->GetNfragZ(0),fAmpTotal[fd+fCalibraMode->GetNfragZ(0)]/fRelativeScale);
1912                 }
1913               }
1914             }
1915           }
1916         }
1917       }
1918       break;
1919     default: break;
1920     }
1921 }
1922 //____________Offine tracking in the AliTRDtracker_____________________________
1923 void AliTRDCalibraFillHisto::FillTheInfoOfTheTrackPH()
1924 {
1925   //
1926   // For the offline tracking or mcm tracklets
1927   // This function will be called in the functions UpdateHistogram...
1928   // to fill the info of a track for the drift velocity  calibration
1929   //
1930
1931   Int_t nb  =  1; // Nombre de zones traversees 1, 2 ou plus de 3
1932   Int_t fd1 = -1; // Premiere zone non nulle
1933   Int_t fd2 = -1; // Deuxieme zone non nulle
1934   Int_t k1  = -1; // Debut de la premiere zone
1935   Int_t k2  = -1; // Debut de la seconde zone
1936
1937   // See if the track goes through different zones
1938   for (Int_t k = 0; k < fTimeMax; k++) {
1939     if (fPHValue[k] > 0.0) {
1940       if (fd1 == -1) {
1941         fd1 = fPHPlace[k];
1942         k1  = k;
1943       }
1944       if (fPHPlace[k] != fd1) {
1945         if (fd2 == -1) {
1946           k2  = k;
1947           fd2 = fPHPlace[k];
1948           nb  = 2;
1949         }
1950         if (fPHPlace[k] != fd2) {
1951           nb = 3;
1952         }
1953       }
1954     }
1955   }
1956
1957
1958   switch(nb)
1959     {
1960     case 1:
1961       fNumberUsedPh[0]++;
1962       for (Int_t i = 0; i < fTimeMax; i++) {
1963         if (fHisto2d) {
1964           fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1965         }
1966         if (fVector2d) {
1967           fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1968         }
1969       }
1970       break;
1971     case 2:
1972       if ((fd1 == fd2+1) ||
1973           (fd2 == fd1+1)) {
1974         // One of the two fast all the think
1975         if (k2 > (k1+fDifference)) {
1976           //we choose to fill the fd1 with all the values
1977           fNumberUsedPh[1]++;
1978           for (Int_t i = 0; i < fTimeMax; i++) {
1979             if (fHisto2d) {
1980               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
1981             }
1982             if (fVector2d) {
1983               fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
1984             }
1985           }
1986         }
1987         if ((k2+fDifference) < fTimeMax) {
1988           //we choose to fill the fd2 with all the values
1989           fNumberUsedPh[1]++;
1990           for (Int_t i = 0; i < fTimeMax; i++) {
1991             if (fHisto2d) {
1992               fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
1993             }
1994           if (fVector2d) {
1995             fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
1996           }
1997           }
1998         }
1999       }
2000       // Two zones voisines sinon rien!
2001       if (fCalibraMode->GetNfragZ(1) > 1) {
2002         // Case 2
2003         if ((fd1+fCalibraMode->GetNfragZ(1)) < (fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))) {
2004           if (fd2 == (fd1+fCalibraMode->GetNfragZ(1))) {
2005             // One of the two fast all the think
2006             if (k2 > (k1+fDifference)) {
2007               //we choose to fill the fd1 with all the values
2008               fNumberUsedPh[1]++;
2009               for (Int_t i = 0; i < fTimeMax; i++) {
2010                 if (fHisto2d) {
2011                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2012                 }
2013                 if (fVector2d) {
2014                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2015                 }
2016               }
2017             }
2018             if ((k2+fDifference) < fTimeMax) {
2019               //we choose to fill the fd2 with all the values
2020               fNumberUsedPh[1]++;
2021               for (Int_t i = 0; i < fTimeMax; i++) {
2022                 if (fHisto2d) {
2023                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2024                 }
2025                 if (fVector2d) {
2026                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2027                 }
2028               }
2029             }
2030           }
2031         }
2032         // Two zones voisines sinon rien!
2033         // Case 3
2034         if ((fd1 - fCalibraMode->GetNfragZ(1)) >= 0) {
2035           if (fd2 == (fd1 - fCalibraMode->GetNfragZ(1))) {
2036             // One of the two fast all the think
2037             if (k2 > (k1 + fDifference)) {
2038               //we choose to fill the fd1 with all the values
2039               fNumberUsedPh[1]++;
2040               for (Int_t i = 0; i < fTimeMax; i++) {
2041                 if (fHisto2d) {
2042                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd1)+0.5,(Float_t) fPHValue[i]);
2043                 }
2044                 if (fVector2d) {
2045                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd1,i,fPHValue[i]);
2046                 }
2047               }
2048             }
2049             if ((k2+fDifference) < fTimeMax) {
2050               //we choose to fill the fd2 with all the values
2051               fNumberUsedPh[1]++;
2052               for (Int_t i = 0; i < fTimeMax; i++) {
2053                 if (fHisto2d) {
2054                   fPH2d->Fill((Float_t) i/fSf,(fCalibraMode->GetXbins(1)+fd2)+0.5,(Float_t) fPHValue[i]);
2055                 }
2056                 if (fVector2d) {
2057                   fCalibraVector->UpdateVectorPH(fDetectorPreviousTrack,fd2,i,fPHValue[i]);
2058                 }
2059               }
2060             }
2061           }
2062         }
2063       }
2064       break;
2065     default: break;
2066     }
2067 }
2068 //____________Offine tracking in the AliTRDtracker_____________________________
2069 Bool_t AliTRDCalibraFillHisto::FindP1TrackPH()
2070 {
2071   //
2072   // For the offline tracking
2073   // This function will be called in the functions UpdateHistogram...
2074   // to fill the find the parameter P1 of a track for the drift velocity  calibration
2075   //
2076
2077
2078   //Number of points: if less than 3 return kFALSE
2079   Int_t Npoints = fListClusters->GetEntriesFast();
2080   if(Npoints <= 2) return kFALSE;
2081
2082   //Variables
2083   TLinearFitter linearFitterTracklet      = TLinearFitter(2,"pol1");        // TLinearFitter per tracklet
2084   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
2085   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
2086   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
2087   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
2088   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
2089   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!
2090   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
2091   Double_t pointError                 = 0.0;                                // error after straight line fit
2092   Int_t    detector                   = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); //detector
2093   Int_t    snpright                   = 1;                                  // if we took in the middle snp
2094   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
2095   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
2096   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
2097   Int_t     nbli                      = 0;                                  // number linear fitter points
2098   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2099
2100   linearFitterTracklet.StoreData(kFALSE);
2101   linearFitterTracklet.ClearPoints();
2102
2103   //if more than one row
2104   Int_t    rowp                       = -1;                              // if it crosses a pad row
2105
2106   //tiltingangle
2107   tiltingangle                        = padplane->GetTiltingAngle();
2108   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2109
2110   //Fill with points
2111   for(Int_t k = 0; k < Npoints; k++){
2112
2113     AliTRDcluster *cl                 = (AliTRDcluster *) fListClusters->At(k);
2114     Double_t ycluster                 = cl->GetY();
2115     Int_t time                        = cl->GetLocalTimeBin();
2116     Double_t timeis                   = time/fSf;
2117     //See if cross two pad rows
2118     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
2119     if(k==0) rowp                     = row;
2120     if(row != rowp) crossrow          = 1;
2121     //Take in the middle of the chamber
2122     //FollowBack
2123     if(time > (Int_t) 10) {
2124     //Follow
2125     //if(time < (Int_t) 11) {
2126       z   = cl->GetZ();
2127       y   = cl->GetY();
2128       snp = fPar2[time];
2129       tgl = fPar3[time];
2130     }
2131     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2132     nbli++;
2133   }
2134   //FollowBack
2135   if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() < 10) snpright = 0;
2136   //Follow
2137   //if(((AliTRDcluster *) fListClusters->At(0))->GetLocalTimeBin() >= 11) snpright = 0;
2138   if(nbli <= 2) return kFALSE;
2139
2140   // Do the straight line fit now
2141   TVectorD pars;
2142   linearFitterTracklet.Eval();
2143   linearFitterTracklet.GetParameters(pars);
2144   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2145   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
2146   dydt  = pars[1];
2147
2148   if( TMath::Abs(snp) <  1.){
2149     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2150   }
2151   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2152
2153   if(fDebugLevel > 0){
2154     if ( !fDebugStreamer ) {
2155       //debug stream
2156       TDirectory *backup = gDirectory;
2157       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2158       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2159     }
2160
2161     (* fDebugStreamer) << "VDRIFT0"<<
2162       "Npoints="<<Npoints<<
2163       "\n";
2164
2165
2166     (* fDebugStreamer) << "VDRIFT"<<
2167       "snpright="<<snpright<<
2168       "Npoints="<<Npoints<<
2169       "nbli="<<nbli<<
2170       "detector="<<detector<<
2171       "snp="<<snp<<
2172       "tnp="<<tnp<<
2173       "tgl="<<tgl<<
2174       "tnt="<<tnt<<
2175       "y="<<y<<
2176       "z="<<z<<
2177       "dydt="<<dydt<<
2178       "dzdx="<<dzdx<<
2179       "crossrow="<<crossrow<<
2180       "errorpar="<<errorpar<<
2181       "pointError="<<pointError<<
2182       "\n";
2183
2184   }
2185
2186   if(Npoints < fNumberClusters) return kFALSE;
2187   if(snpright == 0) return kFALSE;
2188   if(pointError >= 0.1) return kFALSE;
2189   if(crossrow == 1) return kFALSE;
2190
2191   if(fLinearFitterOn){
2192     //Add to the linear fitter of the detector
2193     if( TMath::Abs(snp) <  1.){
2194       Double_t x = tnp-dzdx*tnt;
2195       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2196       if(fLinearFitterDebugOn) {
2197         fLinearVdriftFit->Update(detector,x,pars[1]);
2198       }
2199       fEntriesLinearFitter[detector]++;
2200     }
2201   }
2202   //AliInfo("End of FindP1TrackPH with success!")
2203   return kTRUE;
2204
2205 }
2206 //____________Offine tracking in the AliTRDtracker_____________________________
2207 Bool_t AliTRDCalibraFillHisto::HandlePRF()
2208 {
2209   //
2210   // For the offline tracking
2211   // Fit the tracklet with a line and take the position as reference for the PRF
2212   //
2213
2214   //Number of points
2215   Int_t Npoints  = fListClusters->GetEntriesFast();                         // number of total points
2216   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit
2217   Int_t detector = ((AliTRDcluster *) fListClusters->At(0))->GetDetector(); // detector
2218
2219
2220   // To see the difference due to the fit
2221   Double_t *padPositions;
2222   padPositions = new Double_t[Npoints];
2223   for(Int_t k = 0; k < Npoints; k++){
2224     padPositions[k] = 0.0;
2225   }
2226
2227
2228   //Find the position by a fit
2229   TLinearFitter fitter(2,"pol1");
2230   fitter.StoreData(kFALSE);
2231   fitter.ClearPoints();
2232   for(Int_t k = 0;  k < Npoints; k++){
2233     //Take the cluster
2234     AliTRDcluster *cl  = (AliTRDcluster *) fListClusters->At(k);
2235     Short_t  *signals  = cl->GetSignals();
2236     Double_t     time  = cl->GetLocalTimeBin();
2237     //Calculate x if possible
2238     Float_t xcenter    = 0.0;
2239     Bool_t  echec      = kTRUE;
2240     if((time<=7) || (time>=21)) continue;
2241     // Center 3 balanced: position with the center of the pad
2242     if ((((Float_t) signals[3]) > 0.0) &&
2243         (((Float_t) signals[2]) > 0.0) &&
2244         (((Float_t) signals[4]) > 0.0)) {
2245       echec = kFALSE;
2246       // Security if the denomiateur is 0
2247       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2248            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2249         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2250           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2251                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2252       }
2253       else {
2254         echec = kTRUE;
2255       }
2256     }
2257     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2258     if(echec) continue;
2259     //if no echec: calculate with the position of the pad
2260     // Position of the cluster
2261     Double_t       padPosition = xcenter +  cl->GetPadCol();
2262     padPositions[k]            = padPosition;
2263     Nb3pc++;
2264     fitter.AddPoint(&time, padPosition,1);
2265   }//clusters loop
2266
2267   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
2268   if(Nb3pc < 3) return kFALSE;
2269   fitter.Eval();
2270   TVectorD line(2);
2271   fitter.GetParameters(line);
2272   Float_t  pointError  = -1.0;
2273   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
2274
2275
2276   // Now fill the PRF
2277   for(Int_t k = 0;  k < Npoints; k++){
2278     //Take the cluster
2279     AliTRDcluster *cl      = (AliTRDcluster *) fListClusters->At(k);
2280     Short_t  *signals      = cl->GetSignals();              // signal
2281     Double_t     time      = cl->GetLocalTimeBin();         // time bin
2282     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
2283     Float_t padPos         = cl->GetPadCol();               // middle pad
2284     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
2285     Float_t ycenter        = 0.0;                           // relative center charge
2286     Float_t ymin           = 0.0;                           // relative left charge
2287     Float_t ymax           = 0.0;                           // relative right charge
2288     Double_t tgl           = fPar3[(Int_t)time];            // dz/dl and not dz/dx
2289     Double_t pt            = fPar4[(Int_t)time];            // pt
2290     Float_t  dzdx          = 0.0;                           // dzdx
2291
2292
2293     //Requiere simply two pads clusters at least
2294     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2295        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2296       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2297       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2298       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
2299       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
2300     }
2301
2302     //calibration group
2303     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2304     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2305     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
2306     Double_t  snp          = fPar2[(Int_t)time];                        // sin angle in xy plan
2307     Float_t   xcl          = cl->GetY();                                // y cluster
2308     Float_t   qcl          = cl->GetQ();                                // charge cluster
2309     Int_t     plane        = GetPlane(detector);                        // plane
2310     Int_t     chamber      = GetChamber(detector);                      // chamber
2311     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2312     Double_t  x            = dpad;                                      // reconstructed position moved
2313     Float_t   Ep           = pointError;                                // error of fit
2314     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2315     Float_t   signal3      = (Float_t)signals[3];                       // signal
2316     Float_t   signal2      = (Float_t)signals[2];                       // signal
2317     Float_t   signal4      = (Float_t)signals[4];                       // signal
2318     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2319     Float_t tnp            = 0.0;
2320     if(TMath::Abs(snp) < 1.0){
2321       tnp = snp / (TMath::Sqrt(1-snp*snp));
2322       dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2323     }
2324
2325
2326     if(fDebugLevel > 0){
2327       if ( !fDebugStreamer ) {
2328         //debug stream
2329         TDirectory *backup = gDirectory;
2330         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2331         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2332       }
2333
2334       (* fDebugStreamer) << "PRF0"<<
2335         "caligroup="<<caligroup<<
2336         "detector="<<detector<<
2337         "plane="<<plane<<
2338         "chamber="<<chamber<<
2339         "Npoints="<<Npoints<<
2340         "Np="<<Nb3pc<<
2341         "Ep="<<Ep<<
2342         "snp="<<snp<<
2343         "tnp="<<tnp<<
2344         "tgl="<<tgl<<
2345         "dzdx="<<dzdx<<
2346         "pt="<<pt<<
2347         "padPos="<<padPos<<
2348         "padPosition="<<padPositions[k]<<
2349         "padPosTracklet="<<padPosTracklet<<
2350         "x="<<x<<
2351         "ycenter="<<ycenter<<
2352         "ymin="<<ymin<<
2353         "ymax="<<ymax<<
2354         "xcl="<<xcl<<
2355         "qcl="<<qcl<<
2356         "signal1="<<signal1<<
2357         "signal2="<<signal2<<
2358         "signal3="<<signal3<<
2359         "signal4="<<signal4<<
2360         "signal5="<<signal5<<
2361         "time="<<time<<
2362         "\n";
2363       x = xdiff;
2364       Int_t type=0;
2365       Float_t y = ycenter;
2366       (* fDebugStreamer) << "PRFALL"<<
2367         "caligroup="<<caligroup<<
2368         "detector="<<detector<<
2369         "plane="<<plane<<
2370         "chamber="<<chamber<<
2371         "Npoints="<<Npoints<<
2372         "Np="<<Nb3pc<<
2373         "Ep="<<Ep<<
2374         "type="<<type<<
2375         "snp="<<snp<<
2376         "tnp="<<tnp<<
2377         "tgl="<<tgl<<
2378         "dzdx="<<dzdx<<
2379         "pt="<<pt<<
2380         "padPos="<<padPos<<
2381         "padPosition="<<padPositions[k]<<
2382         "padPosTracklet="<<padPosTracklet<<
2383         "x="<<x<<
2384         "y="<<y<<
2385         "xcl="<<xcl<<
2386         "qcl="<<qcl<<
2387         "signal1="<<signal1<<
2388         "signal2="<<signal2<<
2389         "signal3="<<signal3<<
2390         "signal4="<<signal4<<
2391         "signal5="<<signal5<<
2392         "time="<<time<<
2393         "\n";
2394       x=-(xdiff+1);
2395       y = ymin;
2396       type=-1;
2397       (* fDebugStreamer) << "PRFALL"<<
2398         "caligroup="<<caligroup<<
2399         "detector="<<detector<<
2400         "plane="<<plane<<
2401         "chamber="<<chamber<<
2402         "Npoints="<<Npoints<<
2403         "Np="<<Nb3pc<<
2404         "Ep="<<Ep<<
2405         "type="<<type<<
2406         "snp="<<snp<<
2407         "tnp="<<tnp<<
2408         "tgl="<<tgl<<
2409         "dzdx="<<dzdx<<
2410         "pt="<<pt<<
2411         "padPos="<<padPos<<
2412         "padPosition="<<padPositions[k]<<
2413         "padPosTracklet="<<padPosTracklet<<
2414         "x="<<x<<
2415         "y="<<y<<
2416         "xcl="<<xcl<<
2417         "qcl="<<qcl<<
2418         "signal1="<<signal1<<
2419         "signal2="<<signal2<<
2420         "signal3="<<signal3<<
2421         "signal4="<<signal4<<
2422         "signal5="<<signal5<<
2423         "time="<<time<<
2424         "\n";
2425       x=1-xdiff;
2426       y = ymax;
2427       type=1;
2428
2429       (* fDebugStreamer) << "PRFALL"<<
2430         "caligroup="<<caligroup<<
2431         "detector="<<detector<<
2432         "plane="<<plane<<
2433         "chamber="<<chamber<<
2434         "Npoints="<<Npoints<<
2435         "Np="<<Nb3pc<<
2436         "Ep="<<Ep<<
2437         "type="<<type<<
2438         "snp="<<snp<<
2439         "tnp="<<tnp<<
2440         "tgl="<<tgl<<
2441         "dzdx="<<dzdx<<
2442         "pt="<<pt<<
2443         "padPos="<<padPos<<
2444         "padPosition="<<padPositions[k]<<
2445         "padPosTracklet="<<padPosTracklet<<
2446         "x="<<x<<
2447         "y="<<y<<
2448         "xcl="<<xcl<<
2449         "qcl="<<qcl<<
2450         "signal1="<<signal1<<
2451         "signal2="<<signal2<<
2452         "signal3="<<signal3<<
2453         "signal4="<<signal4<<
2454         "signal5="<<signal5<<
2455         "time="<<time<<
2456         "\n";
2457
2458     }
2459
2460     // some cuts
2461     if(Npoints < fNumberClusters) continue;
2462     if(Nb3pc <= 5) continue;
2463     if((time >= 21) || (time < 7)) continue;
2464     if(TMath::Abs(snp) >= 1.0) continue;
2465     if(qcl < 80) continue;
2466
2467     Bool_t echec   = kFALSE;
2468     Double_t shift = 0.0;
2469     //Calculate the shift in x coresponding to this tnp
2470     if(fNgroupprf != 0.0){
2471       shift      = -3.0*(fNgroupprf-1)-1.5;
2472       Double_t limithigh  = -0.2*(fNgroupprf-1);
2473       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2474       else{
2475         while(tnp > limithigh){
2476           limithigh += 0.2;
2477           shift += 3.0;
2478         }
2479       }
2480     }
2481     if (fHisto2d && !echec) {
2482       if(TMath::Abs(dpad) < 1.5) {
2483         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2484         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2485       }
2486       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2487         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2488         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2489       }
2490       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2491         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2492         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2493       }
2494     }
2495     //Not equivalent anymore here!
2496     if (fVector2d && !echec) {
2497       if(TMath::Abs(dpad) < 1.5) {
2498         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2499         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2500       }
2501       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2502         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2503         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2504       }
2505       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
2506         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2507         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2508       }
2509     }
2510   }
2511   return kTRUE;
2512
2513 }
2514 //____________Offine tracking in the AliTRDtracker_____________________________
2515 Bool_t AliTRDCalibraFillHisto::FindP1TrackPHtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2516 {
2517   //
2518   // For the offline tracking
2519   // This function will be called in the functions UpdateHistogram...
2520   // to fill the find the parameter P1 of a track for the drift velocity  calibration
2521   //
2522
2523
2524   //Number of points: if less than 3 return kFALSE
2525   Int_t Npoints = index1-index0;
2526   if(Npoints <= 2) return kFALSE;
2527
2528   //Variables
2529   TLinearFitter linearFitterTracklet  = TLinearFitter(2,"pol1");            // TLinearFitter per tracklet
2530   Double_t snp                        = 0.0;                                // sin angle in the plan yx track
2531   Double_t y                          = 0.0;                                // y clusters in the middle of the chamber
2532   Double_t z                          = 0.0;                                // z cluster  in the middle of the chamber
2533   Double_t dydt                       = 0.0;                                // dydt tracklet after straight line fit
2534   Double_t tnp                        = 0.0;                                // tan angle in the plan xy track
2535   Double_t tgl                        = 0.0;                                // dz/dl and not dz/dx!
2536   Double_t errorpar                   = 0.0;                                // error after straight line fit on dy/dt
2537   Double_t pointError                 = 0.0;                                // error after straight line fit
2538   Int_t    detector                   = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); //detector
2539   //Int_t    snpright                   = 1;                                  // if we took in the middle snp
2540   Int_t    crossrow                   = 0;                                  // if it crosses a pad row
2541   Double_t  tiltingangle              = 0;                                  // tiltingangle of the pad
2542   Float_t   dzdx                      = 0;                                  // dz/dx now from dz/dl
2543   Int_t     nbli                      = 0;                                  // number linear fitter points
2544   AliTRDpadPlane *padplane            = fGeo->GetPadPlane(GetPlane(detector),GetChamber(detector));
2545
2546   linearFitterTracklet.StoreData(kFALSE);
2547   linearFitterTracklet.ClearPoints();
2548
2549   //if more than one row
2550   Int_t    rowp                       = -1;                              // if it crosses a pad row
2551
2552   //tiltingangle
2553   tiltingangle                        = padplane->GetTiltingAngle();
2554   Float_t  tnt                        = TMath::Tan(tiltingangle/180.*TMath::Pi()); // tan tiltingangle
2555
2556   //Fill with points
2557   for(Int_t k = 0; k < Npoints; k++){
2558
2559     AliTRDcluster *cl                 = (AliTRDcluster *) t->GetCluster(k+index0);
2560     Double_t ycluster                 = cl->GetY();
2561     Int_t time                        = cl->GetLocalTimeBin();
2562     Double_t timeis                   = time/fSf;
2563     //See if cross two pad rows
2564     Int_t    row                      = padplane->GetPadRowNumber(cl->GetZ());
2565     if(k==0) rowp                     = row;
2566     if(row != rowp) crossrow          = 1;
2567     //Take in the middle of the chamber
2568     //FollowBack
2569     //if(time > (Int_t) 10) {
2570     //Follow
2571     if(time < (Int_t) 11) {
2572       z   = cl->GetZ();
2573       y   = cl->GetY();
2574     }
2575     linearFitterTracklet.AddPoint(&timeis,ycluster,1);
2576     nbli++;
2577   }
2578
2579   // take now the snp, tnp and tgl from the track
2580   snp = t->GetSnpPlane(GetPlane(detector));
2581   tgl = t->GetTglPlane(GetPlane(detector));
2582
2583   //FollowBack
2584   //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() < 10) snpright = 0;
2585   //Follow
2586   //if(((AliTRDcluster *) t->GetCluster(index0))->GetLocalTimeBin() >= 11) snpright = 0;
2587   if(nbli <= 2) return kFALSE;
2588
2589   // Do the straight line fit now
2590   TVectorD pars;
2591   linearFitterTracklet.Eval();
2592   linearFitterTracklet.GetParameters(pars);
2593   pointError  =  TMath::Sqrt(linearFitterTracklet.GetChisquare()/nbli);
2594   errorpar    =  linearFitterTracklet.GetParError(1)*pointError;
2595   dydt  = pars[1];
2596
2597   if( TMath::Abs(snp) <  1.){
2598     tnp = snp / (TMath::Sqrt(1-(snp*snp)));
2599   }
2600   dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2601
2602   if(fDebugLevel > 0){
2603     if ( !fDebugStreamer ) {
2604       //debug stream
2605       TDirectory *backup = gDirectory;
2606       fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2607       if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2608     }
2609
2610
2611     (* fDebugStreamer) << "VDRIFT"<<
2612       //"snpright="<<snpright<<
2613       "Npoints="<<Npoints<<
2614       "nbli="<<nbli<<
2615       "detector="<<detector<<
2616       "snp="<<snp<<
2617       "tnp="<<tnp<<
2618       "tgl="<<tgl<<
2619       "tnt="<<tnt<<
2620       "y="<<y<<
2621       "z="<<z<<
2622       "dydt="<<dydt<<
2623       "dzdx="<<dzdx<<
2624       "crossrow="<<crossrow<<
2625       "errorpar="<<errorpar<<
2626       "pointError="<<pointError<<
2627       "\n";
2628
2629   }
2630
2631   if(Npoints < fNumberClusters) return kFALSE;
2632   //if(snpright == 0) return kFALSE;
2633   if(pointError >= 0.1) return kFALSE;
2634   if(crossrow == 1) return kFALSE;
2635
2636   if(fLinearFitterOn){
2637     //Add to the linear fitter of the detector
2638     if( TMath::Abs(snp) <  1.){
2639       Double_t x = tnp-dzdx*tnt;
2640       (GetLinearFitter(detector,kTRUE))->AddPoint(&x,dydt);
2641       if(fLinearFitterDebugOn) {
2642         fLinearVdriftFit->Update(detector,x,pars[1]);
2643       }
2644       fEntriesLinearFitter[detector]++;
2645     }
2646   }
2647   //AliInfo("End of FindP1TrackPH with success!")
2648   return kTRUE;
2649
2650 }
2651 //____________Offine tracking in the AliTRDtracker_____________________________
2652 Bool_t AliTRDCalibraFillHisto::HandlePRFtrack(AliTRDtrack *t, Int_t index0, Int_t index1)
2653 {
2654   //
2655   // For the offline tracking
2656   // Fit the tracklet with a line and take the position as reference for the PRF
2657   //
2658
2659   //Number of points
2660   Int_t Npoints  = index1-index0;                                           // number of total points
2661   Int_t Nb3pc    = 0;                                                       // number of three pads clusters used for fit
2662   Int_t detector = ((AliTRDcluster *) t->GetCluster(index0))->GetDetector(); // detector
2663
2664
2665   // To see the difference due to the fit
2666   Double_t *padPositions;
2667   padPositions = new Double_t[Npoints];
2668   for(Int_t k = 0; k < Npoints; k++){
2669     padPositions[k] = 0.0;
2670   }
2671
2672
2673   //Find the position by a fit
2674   TLinearFitter fitter(2,"pol1");
2675   fitter.StoreData(kFALSE);
2676   fitter.ClearPoints();
2677   for(Int_t k = 0;  k < Npoints; k++){
2678     //Take the cluster
2679     AliTRDcluster *cl  = (AliTRDcluster *) t->GetCluster(k+index0);
2680     Short_t  *signals  = cl->GetSignals();
2681     Double_t     time  = cl->GetLocalTimeBin();
2682     //Calculate x if possible
2683     Float_t xcenter    = 0.0;
2684     Bool_t  echec      = kTRUE;
2685     if((time<=7) || (time>=21)) continue;
2686     // Center 3 balanced: position with the center of the pad
2687     if ((((Float_t) signals[3]) > 0.0) &&
2688         (((Float_t) signals[2]) > 0.0) &&
2689         (((Float_t) signals[4]) > 0.0)) {
2690       echec = kFALSE;
2691       // Security if the denomiateur is 0
2692       if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
2693            ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
2694         xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
2695           / (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
2696                         / ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
2697       }
2698       else {
2699         echec = kTRUE;
2700       }
2701     }
2702     if(TMath::Abs(xcenter) > 0.5) echec = kTRUE;
2703     if(echec) continue;
2704     //if no echec: calculate with the position of the pad
2705     // Position of the cluster
2706     Double_t       padPosition = xcenter +  cl->GetPadCol();
2707     padPositions[k]            = padPosition;
2708     Nb3pc++;
2709     fitter.AddPoint(&time, padPosition,1);
2710   }//clusters loop
2711
2712   //printf("Nb3pc %d, Npoints %d\n",Nb3pc,Npoints);
2713   if(Nb3pc < 3) return kFALSE;
2714   fitter.Eval();
2715   TVectorD line(2);
2716   fitter.GetParameters(line);
2717   Float_t  pointError  = -1.0;
2718   pointError  =  TMath::Sqrt(fitter.GetChisquare()/Nb3pc);
2719
2720   // Take the tgl and snp with the track t now
2721   Double_t  tgl = t->GetTglPlane(GetPlane(detector)); //dz/dl and not dz/dx
2722   Double_t  snp = t->GetSnpPlane(GetPlane(detector)); // sin angle in xy plan
2723   Float_t  dzdx = 0.0;                                // dzdx
2724   Float_t  tnp  = 0.0;
2725   if(TMath::Abs(snp) < 1.0){
2726     tnp = snp / (TMath::Sqrt(1-snp*snp));
2727     dzdx = tgl*TMath::Sqrt(1+tnp*tnp);
2728   }
2729
2730
2731   // Now fill the PRF
2732   for(Int_t k = 0;  k < Npoints; k++){
2733     //Take the cluster
2734     AliTRDcluster *cl      = (AliTRDcluster *) t->GetCluster(k+index0);
2735     Short_t  *signals      = cl->GetSignals();              // signal
2736     Double_t     time      = cl->GetLocalTimeBin();         // time bin
2737     Float_t padPosTracklet = line[0]+line[1]*time;          // reconstruct position from fit
2738     Float_t padPos         = cl->GetPadCol();               // middle pad
2739     Double_t dpad          = padPosTracklet - padPos;       // reconstruct position relative to middle pad from fit
2740     Float_t ycenter        = 0.0;                           // relative center charge
2741     Float_t ymin           = 0.0;                           // relative left charge
2742     Float_t ymax           = 0.0;                           // relative right charge
2743
2744
2745
2746     //Requiere simply two pads clusters at least
2747     if(((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[2]) > 0.0)) ||
2748        ((((Float_t) signals[3]) > 0.0) && (((Float_t) signals[4]) > 0.0))){
2749       Float_t sum     = ((Float_t) signals[2]) + ((Float_t) signals[3]) + ((Float_t) signals[4]);
2750       if(sum > 0.0) ycenter = ((Float_t) signals[3])/ sum;
2751       if(sum > 0.0) ymin    = ((Float_t) signals[2])/ sum;
2752       if(sum > 0.0) ymax    = ((Float_t) signals[4])/ sum;
2753     }
2754
2755     //calibration group
2756     Int_t    *rowcol       = CalculateRowCol(cl);                       // calcul col and row pad of the cluster
2757     Int_t     grouplocal   = CalculateCalibrationGroup(2,rowcol);       // calcul the corresponding group
2758     Int_t     caligroup    = fCalibraMode->GetXbins(2)+ grouplocal;     // calcul the corresponding group
2759     Float_t   xcl          = cl->GetY();                                // y cluster
2760     Float_t   qcl          = cl->GetQ();                                // charge cluster
2761     Int_t     plane        = GetPlane(detector);                        // plane
2762     Int_t     chamber      = GetChamber(detector);                      // chamber
2763     Double_t  xdiff        = dpad;                                      // reconstructed position constant
2764     Double_t  x            = dpad;                                      // reconstructed position moved
2765     Float_t   Ep           = pointError;                                // error of fit
2766     Float_t   signal1      = (Float_t)signals[1];                       // signal at the border
2767     Float_t   signal3      = (Float_t)signals[3];                       // signal
2768     Float_t   signal2      = (Float_t)signals[2];                       // signal
2769     Float_t   signal4      = (Float_t)signals[4];                       // signal
2770     Float_t   signal5      = (Float_t)signals[5];                       // signal at the border
2771
2772
2773
2774     if(fDebugLevel > 0){
2775       if ( !fDebugStreamer ) {
2776         //debug stream
2777         TDirectory *backup = gDirectory;
2778         fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
2779         if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
2780       }
2781
2782       (* fDebugStreamer) << "PRF0"<<
2783         "caligroup="<<caligroup<<
2784         "detector="<<detector<<
2785         "plane="<<plane<<
2786         "chamber="<<chamber<<
2787         "Npoints="<<Npoints<<
2788         "Np="<<Nb3pc<<
2789         "Ep="<<Ep<<
2790         "snp="<<snp<<
2791         "tnp="<<tnp<<
2792         "tgl="<<tgl<<
2793         "dzdx="<<dzdx<<
2794         "padPos="<<padPos<<
2795         "padPosition="<<padPositions[k]<<
2796         "padPosTracklet="<<padPosTracklet<<
2797         "x="<<x<<
2798         "ycenter="<<ycenter<<
2799         "ymin="<<ymin<<
2800         "ymax="<<ymax<<
2801         "xcl="<<xcl<<
2802         "qcl="<<qcl<<
2803         "signal1="<<signal1<<
2804         "signal2="<<signal2<<
2805         "signal3="<<signal3<<
2806         "signal4="<<signal4<<
2807         "signal5="<<signal5<<
2808         "time="<<time<<
2809         "\n";
2810       x = xdiff;
2811       Int_t type=0;
2812       Float_t y = ycenter;
2813       (* fDebugStreamer) << "PRFALL"<<
2814         "caligroup="<<caligroup<<
2815         "detector="<<detector<<
2816         "plane="<<plane<<
2817         "chamber="<<chamber<<
2818         "Npoints="<<Npoints<<
2819         "Np="<<Nb3pc<<
2820         "Ep="<<Ep<<
2821         "type="<<type<<
2822         "snp="<<snp<<
2823         "tnp="<<tnp<<
2824         "tgl="<<tgl<<
2825         "dzdx="<<dzdx<<
2826         "padPos="<<padPos<<
2827         "padPosition="<<padPositions[k]<<
2828         "padPosTracklet="<<padPosTracklet<<
2829         "x="<<x<<
2830         "y="<<y<<
2831         "xcl="<<xcl<<
2832         "qcl="<<qcl<<
2833         "signal1="<<signal1<<
2834         "signal2="<<signal2<<
2835         "signal3="<<signal3<<
2836         "signal4="<<signal4<<
2837         "signal5="<<signal5<<
2838         "time="<<time<<
2839         "\n";
2840       x=-(xdiff+1);
2841       y = ymin;
2842       type=-1;
2843       (* fDebugStreamer) << "PRFALL"<<
2844         "caligroup="<<caligroup<<
2845         "detector="<<detector<<
2846         "plane="<<plane<<
2847         "chamber="<<chamber<<
2848         "Npoints="<<Npoints<<
2849         "Np="<<Nb3pc<<
2850         "Ep="<<Ep<<
2851         "type="<<type<<
2852         "snp="<<snp<<
2853         "tnp="<<tnp<<
2854         "tgl="<<tgl<<
2855         "dzdx="<<dzdx<<
2856         "padPos="<<padPos<<
2857         "padPosition="<<padPositions[k]<<
2858         "padPosTracklet="<<padPosTracklet<<
2859         "x="<<x<<
2860         "y="<<y<<
2861         "xcl="<<xcl<<
2862         "qcl="<<qcl<<
2863         "signal1="<<signal1<<
2864         "signal2="<<signal2<<
2865         "signal3="<<signal3<<
2866         "signal4="<<signal4<<
2867         "signal5="<<signal5<<
2868         "time="<<time<<
2869         "\n";
2870       x=1-xdiff;
2871       y = ymax;
2872       type=1;
2873
2874       (* fDebugStreamer) << "PRFALL"<<
2875         "caligroup="<<caligroup<<
2876         "detector="<<detector<<
2877         "plane="<<plane<<
2878         "chamber="<<chamber<<
2879         "Npoints="<<Npoints<<
2880         "Np="<<Nb3pc<<
2881         "Ep="<<Ep<<
2882         "type="<<type<<
2883         "snp="<<snp<<
2884         "tnp="<<tnp<<
2885         "tgl="<<tgl<<
2886         "dzdx="<<dzdx<<
2887         "padPos="<<padPos<<
2888         "padPosition="<<padPositions[k]<<
2889         "padPosTracklet="<<padPosTracklet<<
2890         "x="<<x<<
2891         "y="<<y<<
2892         "xcl="<<xcl<<
2893         "qcl="<<qcl<<
2894         "signal1="<<signal1<<
2895         "signal2="<<signal2<<
2896         "signal3="<<signal3<<
2897         "signal4="<<signal4<<
2898         "signal5="<<signal5<<
2899         "time="<<time<<
2900         "\n";
2901
2902     }
2903
2904     // some cuts
2905     if(Npoints < fNumberClusters) continue;
2906     if(Nb3pc <= 5) continue;
2907     if((time >= 21) || (time < 7)) continue;
2908     if(TMath::Abs(snp) >= 1.0) continue;
2909     if(qcl < 80) continue;
2910
2911     Bool_t echec   = kFALSE;
2912     Double_t shift = 0.0;
2913     //Calculate the shift in x coresponding to this tnp
2914     if(fNgroupprf != 0.0){
2915       shift      = -3.0*(fNgroupprf-1)-1.5;
2916       Double_t limithigh  = -0.2*(fNgroupprf-1);
2917       if((tnp < (-0.2*fNgroupprf)) || (tnp > (0.2*fNgroupprf))) echec = kTRUE;
2918       else{
2919         while(tnp > limithigh){
2920           limithigh += 0.2;
2921           shift += 3.0;
2922         }
2923       }
2924     }
2925     if (fHisto2d && !echec) {
2926       if(TMath::Abs(dpad) < 1.5) {
2927         fPRF2d->Fill(shift+dpad,(caligroup+0.5),ycenter);
2928         fPRF2d->Fill(shift-dpad,(caligroup+0.5),ycenter);
2929       }
2930       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2931         fPRF2d->Fill(shift-(dpad+1.0),(caligroup+0.5),ymin);
2932         fPRF2d->Fill(shift+(dpad+1.0),(caligroup+0.5),ymin);
2933       }
2934       if((ymax > 0.0) && (TMath::Abs(dpad-1.0) < 1.5)) {
2935         fPRF2d->Fill(shift+1.0-dpad,(caligroup+0.5),ymax);
2936         fPRF2d->Fill(shift-1.0+dpad,(caligroup+0.5),ymax);
2937       }
2938     }
2939     //Not equivalent anymore here!
2940     if (fVector2d && !echec) {
2941       if(TMath::Abs(dpad) < 1.5) {
2942         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+dpad,ycenter);
2943         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-dpad,ycenter);
2944       }
2945       if((ymin > 0.0) && (TMath::Abs(dpad+1.0) < 1.5)) {
2946         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-(dpad+1.0),ymin);
2947         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+(dpad+1.0),ymin);
2948       }
2949       if((ymax > 0.0)  && (TMath::Abs(dpad-1.0) < 1.5)) {
2950         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift+1.0-dpad,ymax);
2951         fCalibraVector->UpdateVectorPRF(fDetectorPreviousTrack,grouplocal,shift-1.0+dpad,ymax);
2952       }
2953     }
2954   }
2955   return kTRUE;
2956
2957 }
2958 //
2959 //____________Some basic geometry function_____________________________________
2960 //
2961 //_____________________________________________________________________________
2962 Int_t AliTRDCalibraFillHisto::GetPlane(Int_t d) const
2963 {
2964   //
2965   // Reconstruct the plane number from the detector number
2966   //
2967
2968   return ((Int_t) (d % 6));
2969
2970 }
2971
2972 //_____________________________________________________________________________
2973 Int_t AliTRDCalibraFillHisto::GetChamber(Int_t d) const
2974 {
2975   //
2976   // Reconstruct the chamber number from the detector number
2977   //
2978   Int_t fgkNplan = 6;
2979
2980   return ((Int_t) (d % 30) / fgkNplan);
2981
2982 }
2983
2984 //_____________________________________________________________________________
2985 Int_t AliTRDCalibraFillHisto::GetSector(Int_t d) const
2986 {
2987   //
2988   // Reconstruct the sector number from the detector number
2989   //
2990   Int_t fg = 30;
2991
2992   return ((Int_t) (d / fg));
2993
2994 }
2995 //_____________________________________________________________________
2996 TProfile2D* AliTRDCalibraFillHisto::GetPH2d(Int_t nbtimebin, Float_t samplefrequency)
2997 {
2998     //
2999     // return pointer to fPH2d TProfile2D
3000     // create a new TProfile2D if it doesn't exist allready
3001     //
3002     if ( fPH2d )
3003         return fPH2d;
3004
3005     // Some parameters
3006     fTimeMax = nbtimebin;
3007     fSf      = samplefrequency;
3008
3009     /*
3010       AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
3011       ,fCalibraMode->GetNz(1)
3012       ,fCalibraMode->GetNrphi(1)));
3013
3014       // Calcul the number of Xbins
3015       Int_t Ntotal1 = 0;
3016       fCalibraMode->ModePadCalibration(2,1);
3017       fCalibraMode->ModePadFragmentation(0,2,0,1);
3018       fCalibraMode->SetDetChamb2(1);
3019       Ntotal1 += 6 * 18 * fCalibraMode->GetDetChamb2(1);
3020       fCalibraMode->ModePadCalibration(0,1);
3021       fCalibraMode->ModePadFragmentation(0,0,0,1);
3022       fCalibraMode->SetDetChamb0(1);
3023       Ntotal1 += 6 * 4 * 18 * fCalibraMode->GetDetChamb0(1);
3024       AliInfo(Form("Total number of Xbins: %d",Ntotal1));
3025
3026       CreatePH2d(Ntotal1);
3027     */
3028
3029     CreatePH2d(540);
3030
3031     return fPH2d;
3032 }
3033 //_____________________________________________________________________
3034 TLinearFitter* AliTRDCalibraFillHisto::GetLinearFitter(Int_t detector, Bool_t force)
3035 {
3036     //
3037     // return pointer to TLinearFitter Calibration
3038     // if force is true create a new TLinearFitter if it doesn't exist allready
3039     //
3040
3041   if ((!force) || (fLinearFitterArray.UncheckedAt(detector))){
3042     return (TLinearFitter*)fLinearFitterArray.UncheckedAt(detector);
3043   }
3044
3045   // if we are forced and TLinearFitter doesn't yet exist create it
3046
3047   // new TLinearFitter
3048   TLinearFitter *linearfitter = new TLinearFitter(2,"pol1");
3049   fLinearFitterArray.AddAt(linearfitter,detector);
3050   return linearfitter;
3051 }
3052
3053 //_____________________________________________________________________
3054 void  AliTRDCalibraFillHisto::FillCH2d(Int_t x, Float_t y)
3055 {
3056   //
3057   // FillCH2d: Marian style
3058   //
3059
3060   //skip simply the value out of range
3061   if((y>=300.0) || (y<0.0)) return;
3062
3063   //Calcul the y place
3064   Int_t yplace = (Int_t) (fNumberBinCharge*y/300.0)+1;
3065   Int_t place = (fNumberBinCharge+2)*(x+1)+yplace;
3066
3067   //Fill
3068   fCH2d->GetArray()[place]++;
3069
3070 }
3071
3072 //____________________________________________________________________________
3073 void AliTRDCalibraFillHisto::AnalyseLinearFitter()
3074 {
3075   //
3076   // Analyse array of linear fitter because can not be written
3077   // Store two arrays: one with the param the other one with the error param + number of entries
3078   //
3079
3080   for(Int_t k = 0; k < 540; k++){
3081     TLinearFitter *linearfitter = GetLinearFitter(k);
3082     if((linearfitter!=0) && (fEntriesLinearFitter[k]>10)){
3083       TVectorD  *par  = new TVectorD(2);
3084       TVectorD   pare = TVectorD(2);
3085       TVectorD  *parE = new TVectorD(3);
3086       linearfitter->Eval();
3087       linearfitter->GetParameters(*par);
3088       linearfitter->GetErrors(pare);
3089       Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter->GetChisquare())/fEntriesLinearFitter[k]);
3090       (*parE)[0] = pare[0]*ppointError;
3091       (*parE)[1] = pare[1]*ppointError;
3092       (*parE)[2] = (Double_t) fEntriesLinearFitter[k];
3093       ((TObjArray *)fLinearVdriftFit->GetPArray())->AddAt(par,k);
3094       ((TObjArray *)fLinearVdriftFit->GetEArray())->AddAt(parE,k);
3095     }
3096   }
3097 }
3098 //________________________________________________________________________________
3099 Int_t AliTRDCalibraFillHisto::Arrondi(Double_t x) const
3100 {
3101    // Partie entiere of the (x+0.5)
3102
3103   int i;
3104   if (x >= (-0.5)) {
3105     i = int(x + 0.5);
3106     //if (x + 0.5 == Float_t(i) && i & 1) i--;
3107   } else {
3108     i = int(x - 0.5);
3109     //if (x - 0.5 == Float_t(i) && i & 1) i++;
3110     if((x-0.5)==i) i++;
3111
3112   }
3113   return i;
3114 }