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