TRD/AliTRDCalibraFillHisto.cxx

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