TRD/AliTRDCalibraFillHisto.cxx

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