TRD/AliTRDCalibraFillHisto.cxx

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