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