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