/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
* Contributors are mentioned in the code where appropriate. *
* *
* Permission to use, copy, modify and distribute this software and its *
* documentation strictly for non-commercial purposes is hereby granted *
* without fee, provided that the above copyright notice appears in all *
* copies and that both the copyright notice and this permission notice *
* appear in the supporting documentation. The authors make no claims *
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
/* $Id$ */
/////////////////////////////////////////////////////////////////////////////////
//
// AliTRDCalibra
//
// This class is for the TRD calibration of the relative gain factor, the drift velocity,
// the time 0 and the pad response function.
// It can be used for the calibration per chamber but also per group of pads and eventually per pad.
// The user has to choose with the functions SetNz and SetNrphi the precision of the calibration.
//Begin_Html
/*
Nz | 0 | 1 | 2 | 3 | 4 |
group of row pads per detector | 1 | 2 | 4 | 6(chamb2) 8(others chambers) | 12 (chamb2) 16 (chamb0) |
row pads per group | 12 (chamb2) 16 (chamb0) | 6 (chamb2) 8 (chamb0) | 3 (chamb2) 4 (chamb0) | 2 | 1 |
~distance [cm] | 106 (chamb2) 130 (chamb0) | 53 (chamb2) 65 (chamb0) | 26.5 (chamb2) 32.5 (chamb0) | 17 (chamb2) 17 (chamb0) | 9 (chamb2) 9 (chamb0) |
In the z direction
Nrphi | 0 | 1 | 2 | 3 | 4 | 5 | 6 |
group of col pads per detector | 1 | 2 | 4 | 8 | 16 | 36 | 144 |
col pads per group | 144 | 72 | 36 | 18 | 9 | 4 | 1 |
~distance [cm] | 113.4 | 56.7 | 25.3 | 14.3 | 7.25 | 3.2 | 0.8 |
In the rphi direction
*/
//End_Html
//
// Fill histograms or vectors
//----------------------------
//
// 2D Histograms (Histo2d) or vectors (Vector2d), then converted in Trees, will be filled
// from RAW DATA in a run or from reconstructed TRD tracks during the offline tracking
// in the function "FollowBackProlongation" (AliTRDtracker)
// Per default the functions to fill are off.
//
//Begin_Html
/*
Example of 2D histos for the relative gain (left) and the drift velocity (right) calibration of the sector 13
*/
//End_Html
//
// Fit the histograms to find the coefficients
//---------------------------------------------
//
// These 2D histograms or vectors (first converted in 1D histos) will be fitted
// if they have enough entries, otherwise the (default) value of the choosen database
// will be put. For the relative gain calibration the resulted factors will be globally
// normalized to the gain factors of the choosen database. It unables to precise
// previous calibration procedure.
// The function SetDebug enables the user to see:
// _fDebug = 0: nothing, only the values are written in the tree if wanted
// _fDebug = 1: a comparaison of the coefficients found and the default values
// in the choosen database.
// fCoef , histogram of the coefs as function of the calibration group number
// fDelta , histogram of the relative difference of the coef with the default
// value in the database as function of the calibration group number
// fError , dirstribution of this relative difference
// _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir)
// _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the
// pad row and col number
// _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with
// also the comparaison histograms of the 1 for this detector
//
//Begin_Html
/*
Example of fCoef for the relative gain calibration of the sector 13
Example of fDelta (right) and fError (left) for the relative gain calibration of the sector 13
*/
//End_Html
//
// Author:
// R. Bailhache (R.Bailhache@gsi.de)
//
//////////////////////////////////////////////////////////////////////////////////////
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "AliLog.h"
#include "AliCDBManager.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliRawReaderFile.h"
#include "AliRawReader.h"
#include "AliTRDCalibra.h"
#include "AliTRDcalibDB.h"
#include "AliTRDCommonParam.h"
#include "AliTRDmcmTracklet.h"
#include "AliTRDpadPlane.h"
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
#include "AliTRDdigit.h"
#include "AliTRDdigitsManager.h"
#include "AliTRD.h"
#include "AliTRDgeometry.h"
#include "./Cal/AliTRDCalROC.h"
#include "./Cal/AliTRDCalPad.h"
#include "./Cal/AliTRDCalDet.h"
#include "AliTRDrawData.h"
ClassImp(AliTRDCalibra)
AliTRDCalibra* AliTRDCalibra::fgInstance = 0;
Bool_t AliTRDCalibra::fgTerminated = kFALSE;
//_____________singleton implementation_________________________________________________
AliTRDCalibra *AliTRDCalibra::Instance()
{
//
// Singleton implementation
//
if (fgTerminated != kFALSE) {
return 0;
}
if (fgInstance == 0) {
fgInstance = new AliTRDCalibra();
}
return fgInstance;
}
//______________________________________________________________________________________
void AliTRDCalibra::Terminate()
{
//
// Singleton implementation
// Deletes the instance of this class
//
fgTerminated = kTRUE;
if (fgInstance != 0) {
delete fgInstance;
fgInstance = 0;
}
}
//______________________________________________________________________________________
AliTRDCalibra::AliTRDCalibra()
:TObject()
,fMITracking(kFALSE)
,fMcmTracking(kFALSE)
,fMcmCorrectAngle(kFALSE)
,fCH2dOn(kFALSE)
,fPH2dOn(kFALSE)
,fPRF2dOn(kFALSE)
,fHisto2d(kFALSE)
,fVector2d(kFALSE)
,fRelativeScale(0)
,fCountRelativeScale(0)
,fRelativeScaleAuto(kFALSE)
,fThresholdDigit(0)
,fThresholdClusterPRF1(0.0)
,fThresholdClusterPRF2(0.0)
,fCenterOfflineCluster(kFALSE)
,fTraMaxPad(kFALSE)
,fWriteNameCoef(0)
,fWriteName(0)
,fFitPHOn(kFALSE)
,fFitPHPeriode(0)
,fBeginFitCharge(0.0)
,fRangeFitPRF(0.0)
,fMeanChargeOn(kFALSE)
,fFitChargeBisOn(kFALSE)
,fT0Shift(0.0)
,fAccCDB(kFALSE)
,fNumberFit(0)
,fNumberEnt(0)
,fStatisticMean(0.0)
,fDebug(0)
,fFitVoir(0)
,fPRF(0)
,fGain(0)
,fT0(0)
,fVdrift(0)
,fVdriftDetector(0)
,fVdriftPad(0x0)
,fT0Detector(0)
,fT0Pad(0x0)
,fPRFDetector(0)
,fPRFPad(0x0)
,fCoefCH(0x0)
,fDetectorAliTRDtrack(kFALSE)
,fChamberAliTRDtrack(-1)
,fDetectorPreviousTrack(-1)
,fGoodTrack(kTRUE)
,fAmpTotal(0x0)
,fPHPlace(0x0)
,fPHValue(0x0)
,fNumberClusters(0)
,fProcent(0.0)
,fDifference(0)
,fNumberTrack(0)
,fCoefPRFE(0)
,fCoefPRFDB(0)
,fTimeMax(0)
,fSf(0.0)
,fScaleFitFactor(0.0)
,fMinEntries(0)
,fEntriesCurrent(0)
,fL3P0(0.0)
,fL3P2(0.0)
,fG3P2(0.0)
,fVectorPH(0)
,fPlaPH(0)
,fNumberBinCharge(0)
,fVectorCH(0)
,fPlaCH(0)
,fVectorFitCH(0)
,fNumberBinPRF(0)
,fVectorPRF(0)
,fPlaPRF(0)
,fPH2d(0x0)
,fPRF2d(0x0)
,fCH2d(0x0)
,fRebin(0)
{
//
// Default constructor
//
for (Int_t i = 0; i < 3; i++) {
fNz[i] = 0;
fNrphi[i] = 0;
}
for (Int_t k = 0; k < 3; k++) {
fNtotal[k] = 0;
fDetChamb2[k] = 0;
fDetChamb0[k] = 0;
}
// Write
for (Int_t i = 0; i < 3; i++) {
fWriteCoef[i] = kFALSE;
fWrite[i] = kFALSE;
}
// Debug Mode
for (Int_t k = 0; k < 3; k++) {
fDet[k] = 0;
}
for (Int_t i = 0; i < 3; i++) {
fPhd[i] = 0.0;
}
// Init
Init();
}
//______________________________________________________________________________________
AliTRDCalibra::AliTRDCalibra(const AliTRDCalibra &c)
:TObject(c)
,fMITracking(kFALSE)
,fMcmTracking(kFALSE)
,fMcmCorrectAngle(kFALSE)
,fCH2dOn(kFALSE)
,fPH2dOn(kFALSE)
,fPRF2dOn(kFALSE)
,fHisto2d(kFALSE)
,fVector2d(kFALSE)
,fRelativeScale(0)
,fCountRelativeScale(0)
,fRelativeScaleAuto(kFALSE)
,fThresholdDigit(0)
,fThresholdClusterPRF1(0.0)
,fThresholdClusterPRF2(0.0)
,fCenterOfflineCluster(kFALSE)
,fTraMaxPad(kFALSE)
,fWriteNameCoef(0)
,fWriteName(0)
,fFitPHOn(kFALSE)
,fFitPHPeriode(0)
,fBeginFitCharge(0.0)
,fRangeFitPRF(0.0)
,fMeanChargeOn(kFALSE)
,fFitChargeBisOn(kFALSE)
,fT0Shift(0.0)
,fAccCDB(kFALSE)
,fNumberFit(0)
,fNumberEnt(0)
,fStatisticMean(0.0)
,fDebug(0)
,fFitVoir(0)
,fPRF(0)
,fGain(0)
,fT0(0)
,fVdrift(0)
,fVdriftDetector(0)
,fVdriftPad(0x0)
,fT0Detector(0)
,fT0Pad(0x0)
,fPRFDetector(0)
,fPRFPad(0x0)
,fCoefCH(0x0)
,fDetectorAliTRDtrack(kFALSE)
,fChamberAliTRDtrack(-1)
,fDetectorPreviousTrack(-1)
,fGoodTrack(kTRUE)
,fAmpTotal(0x0)
,fPHPlace(0x0)
,fPHValue(0x0)
,fNumberClusters(0)
,fProcent(0.0)
,fDifference(0)
,fNumberTrack(0)
,fCoefPRFE(0)
,fCoefPRFDB(0)
,fTimeMax(0)
,fSf(0.0)
,fScaleFitFactor(0.0)
,fMinEntries(0)
,fEntriesCurrent(0)
,fL3P0(0.0)
,fL3P2(0.0)
,fG3P2(0.0)
,fVectorPH(0)
,fPlaPH(0)
,fNumberBinCharge(0)
,fVectorCH(0)
,fPlaCH(0)
,fVectorFitCH(0)
,fNumberBinPRF(0)
,fVectorPRF(0)
,fPlaPRF(0)
,fPH2d(0x0)
,fPRF2d(0x0)
,fCH2d(0x0)
,fRebin(0)
{
//
// Copy constructor
//
}
//____________________________________________________________________________________
AliTRDCalibra::~AliTRDCalibra()
{
//
// AliTRDCalibra destructor
//
ClearHistos();
ClearTree();
}
//_____________________________________________________________________________
void AliTRDCalibra::Destroy()
{
//
// Delete instance
//
if (fgInstance) {
delete fgInstance;
fgInstance = 0x0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::ClearHistos()
{
//
// Delete the histos
//
if (fPH2d) {
delete fPH2d;
fPH2d = 0x0;
}
if (fCH2d) {
delete fCH2d;
fCH2d = 0x0;
}
if (fPRF2d) {
delete fPRF2d;
fPRF2d = 0x0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::ClearTree()
{
//
// Delete the trees
//
if (fPRF) {
delete fPRF;
fPRF = 0x0;
}
if (fGain) {
delete fGain;
fGain = 0x0;
}
if (fT0) {
delete fT0;
fT0 = 0x0;
}
if (fVdrift) {
delete fVdrift;
fVdrift = 0x0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::Init()
{
//
// Init some default values
//
// How to fill the 2D
fThresholdDigit = 5;
fThresholdClusterPRF1 = 2.0;
fThresholdClusterPRF2 = 3.0;
// Store the Info
fNumberBinCharge = 100;
fNumberBinPRF = 20;
// Write
fWriteName = "TRD.calibration.root";
fWriteNameCoef = "TRD.coefficient.root";
// Fit
fFitPHPeriode = 1;
fBeginFitCharge = 3.5;
fRangeFitPRF = 0.5;
fMinEntries = 800;
fT0Shift = 0.143397;
// Internal variables
// Fill the 2D histos in the offline tracking
fDetectorPreviousTrack = -1;
fChamberAliTRDtrack = -1;
fGoodTrack = kTRUE;
fProcent = 6.0;
fDifference = 17;
fNumberClusters = 18;
fNumberTrack = 0;
fNumberUsedCh[0] = 0;
fNumberUsedCh[1] = 0;
fNumberUsedPh[0] = 0;
fNumberUsedPh[1] = 0;
// Variables in the loop
for (Int_t k = 0; k < 4; k++) {
fChargeCoef[k] = 1.0;
fVdriftCoef[k] = 1.5;
fT0Coef[k] = -1.0;
}
for (Int_t i = 0; i < 2; i++) {
fPRFCoef[i] = -1.0;
}
// Pad calibration
for (Int_t i = 0; i < 3; i++) {
fRowMin[i] = -1;
fRowMax[i] = -1;
fColMax[i] = -1;
fColMin[i] = -1;
fNnZ[i] = -1;
fNnRphi[i] = -1;
fNfragZ[i] = -1;
fNfragRphi[i] = -1;
fXbins[i] = -1;
}
// Local database to be changed
fRebin = 1;
}
//____________Functions fit Online CH2d________________________________________
Bool_t AliTRDCalibra::FitCHOnline(TH2I *ch)
{
//
// Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
// calibration group normalized the resulted coefficients (to 1 normally)
// and write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
TAxis *xch = ch->GetXaxis();
Int_t nbins = xch->GetNbins();
TAxis *yph = ch->GetYaxis();
Int_t nybins = yph->GetNbins();
if (!InitFit(nbins,0)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorCH) {
fVectorCH->Clear();
}
if (fPlaCH) {
fPlaCH->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(0);
// Beginning of the loop betwwen dect1 and dect2
for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
TH1I *projch = (TH1I *) ch->ProjectionY("projch",idect+1,idect+1,(Option_t *)"e");
projch->SetDirectory(0);
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,0);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect, 0);
// Number of entries for this calibration group
Double_t nentries = 0.0;
for (Int_t k = 0; k < nybins; k++) {
nentries += ch->GetBinContent(ch->GetBin(idect+1,k+1));
}
if (nentries > 0) {
fNumberEnt++;
}
// Rebin and statistic stuff
// Rebin
if (fRebin > 1) {
projch = ReBin((TH1I *) projch);
}
// This detector has not enough statistics or was off
if (nentries < fMinEntries) {
// Fill with the default infos
NotEnoughStatistic(idect,0);
// Memory!!!
if (fDebug != 2) {
delete projch;
}
continue;
}
// Statistics of the group fitted
AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
fStatisticMean += nentries;
fNumberFit++;
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitChargeBisOn) {
FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefChargeDB();
}
// Fill Infos Fit
FillInfosFit(idect,0);
// Memory!!!
if (fDebug != 2) {
delete projch;
}
} // Boucle object
// Normierungcharge
if (fDebug != 2) {
NormierungCharge();
}
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)) {
PlotWriteCH();
}
if ((fDebug == 4) ||
(fDebug == 3)) {
PlotCHDB();
}
// Mean Statistic
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean/fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
ConvertVectorFitCHTree();
if (fWriteCoef[0]) {
WriteFitInfos(0);
}
return kTRUE;
}
//____________Functions fit Online CH2d________________________________________
Bool_t AliTRDCalibra::FitCHOnline()
{
//
// Reconstruct a 1D histo from the vectorCH for each calibration group,
// fit the histo, normalized the resulted coefficients (to 1 normally)
// and write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,0)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// Init fCountDet and fCount
InitfCountDetAndfCount(0);
// Beginning of the loop between dect1 and dect2
for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
// Search if the group is in the VectorCH
Int_t place = SearchInVector(idect,0);
// Is in
TH1F *projch = 0x0;
TString name("CH");
name += idect;
if (place != -1) {
// Variable
AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
// Retrieve
fCHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
projch = ConvertVectorCTHisto(fCHInfo,(const char *) name);
projch->SetDirectory(0);
delete fCHInfo;
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,0);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,0);
// Number of entries
Double_t nentries = 0.0;
if (projch) {
for (Int_t k = 0; k < fNumberBinCharge; k++) {
nentries += projch->GetBinContent(k+1);
}
}
if (nentries > 0) {
fNumberEnt++;
}
// Rebin and statistic stuff
// Rebin
if ((fRebin > 1) &&
(place != -1)) {
projch = ReBin((TH1F *) projch);
}
// This detector has not enough statistics or was not found in VectorCH
if ((place == -1) ||
((place != -1) &&
(nentries < fMinEntries))) {
// Fill with the default infos
NotEnoughStatistic(idect,0);
// Memory!!!
if (fDebug != 2) {
delete projch;
}
continue;
}
// Statistic of the histos fitted
AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
fNumberFit++;
fStatisticMean += nentries;
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitChargeBisOn) {
FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefChargeDB();
}
// Fill Infos Fit
FillInfosFit(idect,0);
// Memory!!!
if (fDebug != 2) {
delete projch;
}
} // Boucle object
// Normierungcharge
if (fDebug != 2) {
NormierungCharge();
}
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)){
PlotWriteCH();
}
if((fDebug == 4) ||
(fDebug == 3)){
PlotCHDB();
}
// Mean Statistics
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
ConvertVectorFitCHTree();
if (fWriteCoef[0]) {
WriteFitInfos(0);
}
return kTRUE;
}
//____________Functions fit Online CH2d________________________________________
Bool_t AliTRDCalibra::FitCHOnline(TTree *tree)
{
//
// Look if the calibration group can be found in the tree, if yes take the
// histo, fit it, normalized the resulted coefficients (to 1 normally) and
// write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,0)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorCH) {
fVectorCH->Clear();
}
if (fPlaCH) {
fPlaCH->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(0);
TH1F *projch = 0x0;
tree->SetBranchAddress("histo",&projch);
TObjArray *vectorplace = ConvertTreeVector(tree);
// Beginning of the loop between dect1 and dect2
for (Int_t idect = fDect1[0]; idect < fDect2[0]; idect++) {
//Search if the group is in the VectorCH
Int_t place = SearchInTreeVector(vectorplace,idect);
// Is in
if (place != -1) {
// Variable
tree->GetEntry(place);
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,0);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,0);
// Number of entries
Double_t nentries = 0.0;
if (projch) {
for (Int_t k = 0; k < projch->GetXaxis()->GetNbins(); k++) {
nentries += projch->GetBinContent(k+1);
}
}
if (nentries > 0) {
fNumberEnt++;
}
// Rebin and statistic stuff
// Rebin
if ((fRebin > 1) &&
(place != -1)) {
projch = ReBin((TH1F *) projch);
}
// This detector has not enough statistics or was not found in VectorCH
if((place == -1) ||
((place != -1) &&
(nentries < fMinEntries))) {
// Fill with the default infos
NotEnoughStatistic(idect,0);
continue;
}
// Statistics of the group fitted
AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
fNumberFit++;
fStatisticMean += nentries;
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitChargeBisOn) {
FitBisCH((TH1 *) projch,(Int_t) (idect-fDect1[0]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefChargeDB();
}
// Fill Infos Fit
FillInfosFit(idect,0);
} // Boucle object
// Normierungcharge
if (fDebug != 2) {
NormierungCharge();
}
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)){
PlotWriteCH();
}
if ((fDebug == 4) ||
(fDebug == 3)){
PlotCHDB();
}
// Mean Statistic
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
ConvertVectorFitCHTree();
if (fWriteCoef[0]) {
WriteFitInfos(0);
}
return kTRUE;
}
//________________functions fit Online PH2d____________________________________
Bool_t AliTRDCalibra::FitPHOnline(TProfile2D *ph)
{
//
// Take the 1D profiles (average pulse height), projections of the 2D PH
// on the Xaxis, for each calibration group
// Fit or use the slope of the average pulse height to reconstruct the
// drift velocity write the results in a tree
// A first calibration of T0 is also made using the same method (slope method)
//
// Number of Xbins (detectors or groups of pads)
TAxis *xph = ph->GetXaxis();
TAxis *yph = ph->GetYaxis();
Int_t nbins = xph->GetNbins();
Int_t nybins = yph->GetNbins();
if (!InitFit(nbins,1)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorPH) {
fVectorPH->Clear();
}
if (fPlaPH) {
fPlaPH->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(1);
// Beginning of the loop
for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
TH1D *projph = (TH1D *) ph->ProjectionY("projph",idect+1,idect+1,(Option_t *) "e");
projph->SetDirectory(0);
// Number of entries for this calibration group
Double_t nentries = 0;
for (Int_t k = 0; k < nybins; k++) {
nentries += ph->GetBinEntries(ph->GetBin(idect+1,k+1));
}
if (nentries > 0) {
fNumberEnt++;
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,1);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,1);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if (nentries < fMinEntries) {
// Fill with the default values
NotEnoughStatistic(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
continue;
}
// Statistics of the histos fitted
AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
fNumberFit++;
fStatisticMean += nentries;
// Calcul of "real" coef
CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitPHOn) {
FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefVdriftDB();
FillCoefT0DB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
} // Boucle object
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)) {
PlotWritePH();
PlotWriteT0();
}
if ((fDebug == 4) ||
(fDebug == 3)) {
PlotPHDB();
PlotT0DB();
}
// Mean Statistic
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if(fWriteCoef[1]) {
WriteFitInfos(1);
}
return kTRUE;
}
//____________Functions fit Online PH2d________________________________________
Bool_t AliTRDCalibra::FitPHOnline()
{
//
// Reconstruct the average pulse height from the vectorPH for each
// calibration group
// Fit or use the slope of the average pulse height to reconstruct the
// drift velocity write the results in a tree
// A first calibration of T0 is also made using the same method (slope method)
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,1)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// Init fCountDet and fCount
InitfCountDetAndfCount(1);
// Beginning of the loop
for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
// Search if the group is in the VectorCH
Int_t place = SearchInVector(idect,1);
// Is in
TH1F *projph = 0x0;
TString name("PH");
name += idect;
if (place != -1) {
//Entries
fNumberEnt++;
// Variable
AliTRDPInfo *fPHInfo = new AliTRDPInfo();
// Retrieve
fPHInfo = (AliTRDPInfo *) fVectorPH->At(place);
projph = CorrectTheError((TGraphErrors *) ConvertVectorPHisto(fPHInfo,(const char *) name));
projph->SetDirectory(0);
delete fPHInfo;
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,1);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,1);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if ((place == -1) ||
((place != -1) &&
(fEntriesCurrent < fMinEntries))) {
// Fill with the default values
NotEnoughStatistic(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
continue;
}
// Statistic of the histos fitted
AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
fNumberFit++;
fStatisticMean += fEntriesCurrent;
// Calcul of "real" coef
CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitPHOn) {
FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefVdriftDB();
FillCoefT0DB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
} // Boucle object
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)) {
PlotWritePH();
PlotWriteT0();
}
if ((fDebug == 4) ||
(fDebug == 3)) {
PlotPHDB();
PlotT0DB();
}
// Mean Statistic
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if (fWriteCoef[1]) {
WriteFitInfos(1);
}
return kTRUE;
}
//____________Functions fit Online PH2d________________________________________
Bool_t AliTRDCalibra::FitPHOnline(TTree *tree)
{
//
// Look if the calibration group can be found in the tree, if yes take the
// histo, fit it, and write the results in a tree
// A first calibration of T0 is also made using the same method (slope method)
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,1)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorPH) {
fVectorPH->Clear();
}
if (fPlaPH) {
fPlaPH->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(1);
TGraphErrors *projphtree = 0x0;
tree->SetBranchAddress("histo",&projphtree);
TObjArray *vectorplace = ConvertTreeVector(tree);
// Beginning of the loop
for (Int_t idect = fDect1[1]; idect < fDect2[1]; idect++) {
// Search if the group is in the VectorCH
Int_t place = SearchInTreeVector(vectorplace,idect);
TH1F *projph = 0x0;
// Is in
if (place != -1) {
//Entries
fNumberEnt++;
// Variable
tree->GetEntry(place);
projph = CorrectTheError(projphtree);
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,1);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,1);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if((place == -1) ||
((place != -1) &&
(fEntriesCurrent < fMinEntries))) {
// Fill with the default values
NotEnoughStatistic(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
continue;
}
// Statistics of the histos fitted
AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
fNumberFit++;
fStatisticMean += fEntriesCurrent;
// Calcul of "real" coef
CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
CalculT0CoefMean(fCountDet[1],(Int_t) (idect - fDect1[1]));
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPente((TH1 *) projph,(Int_t) (idect - fDect1[1]));
// Method fit bis
// idect is egal for fDebug = 0 and 2, only to fill the hist
if (fFitPHOn) {
FitPH((TH1 *) projph,(Int_t) (idect - fDect1[1]));
}
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefVdriftDB();
FillCoefT0DB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,1);
// Memory!!!
if (fDebug != 2) {
delete projph;
}
} // Boucle object
// Plot
// 0 no plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)){
PlotWritePH();
PlotWriteT0();
}
if ((fDebug == 4) ||
(fDebug == 3)){
PlotPHDB();
PlotT0DB();
}
// Mean Statistics
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if (fWriteCoef[1]) {
WriteFitInfos(1);
}
return kTRUE;
}
//____________Functions fit Online PRF2d_______________________________________
Bool_t AliTRDCalibra::FitPRFOnline(TProfile2D *prf)
{
//
// Take the 1D profiles (pad response function), projections of the 2D PRF
// on the Xaxis, for each calibration group
// Fit with a gaussian to reconstruct the sigma of the pad response function
// write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
TAxis *xprf = prf->GetXaxis();
TAxis *yprf = prf->GetYaxis();
Int_t nybins = yprf->GetNbins();
Int_t nbins = xprf->GetNbins();
if (!InitFit(nbins,2)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorPRF) {
fVectorPRF->Clear();
}
if (fPlaPRF) {
fPlaPRF->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(2);
// Beginning of the loop
for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
TH1D *projprf = (TH1D *) prf->ProjectionY("projprf",idect+1,idect+1,(Option_t *) "e");
projprf->SetDirectory(0);
// Number of entries for this calibration group
Double_t nentries = 0;
for (Int_t k = 0; k < nybins; k++) {
nentries += prf->GetBinEntries(prf->GetBin(idect+1,k+1));
}
if(nentries > 0) fNumberEnt++;
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,2);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,2);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if (nentries < fMinEntries) {
// Fill with the default values
NotEnoughStatistic(idect,2);
// Memory!
if (fDebug != 2) {
delete projprf;
}
continue;
}
// Statistics of the histos fitted
AliInfo(Form("For the group number %d there are %f stats",idect,nentries));
fNumberFit++;
fStatisticMean += nentries;
// Calcul of "real" coef
if ((fDebug == 1) ||
(fDebug == 4)) {
CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
}
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefPRFDB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,2);
// Memory!!!
if (fDebug != 2) {
delete projprf;
}
} // Boucle object
// Plot
// No plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)) {
PlotWritePRF();
}
if ((fDebug == 4) ||
(fDebug == 3)){
PlotPRFDB();
}
// Mean Statistic
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if (fWriteCoef[2]) {
WriteFitInfos(2);
}
return kTRUE;
}
//____________Functions fit Online PRF2d_______________________________________
Bool_t AliTRDCalibra::FitPRFOnline(TTree *tree)
{
//
// Look if the calibration group can be found in the tree, if yes take
// the histo, fit it, and write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,2)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// For memory
if (fVectorPRF) {
fVectorPRF->Clear();
}
if (fPlaPRF) {
fPlaPRF->Clear();
}
// Init fCountDet and fCount
InitfCountDetAndfCount(2);
TGraphErrors *projprftree = 0x0;
tree->SetBranchAddress("histo",&projprftree);
TObjArray *vectorplace = ConvertTreeVector(tree);
// Beginning of the loop
for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
// Search if the group is in the VectorCH
Int_t place = SearchInTreeVector(vectorplace,idect);
// Is in
TH1F *projprf = 0x0;
if (place != -1) {
//Entries
fNumberEnt++;
// Variable
tree->GetEntry(place);
projprf = CorrectTheError(projprftree);
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,2);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,2);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if ((place == -1) ||
((place != -1) &&
(fEntriesCurrent < fMinEntries))) {
// Fill with the default values
NotEnoughStatistic(idect,2);
// Memory!!!
if (fDebug != 2) {
delete projprf;
}
continue;
}
// Statistics of the histos fitted
AliInfo(Form("For the group number %d there are %d stats",idect,fEntriesCurrent));
fNumberFit++;
fStatisticMean += fEntriesCurrent;
// Calcul of "real" coef
if ((fDebug == 1) ||
(fDebug == 4)){
CalculPRFCoefMean(fCountDet[2],(Int_t) (idect - fDect1[2]));
}
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPRF((TH1 *) projprf,(Int_t) (idect - fDect1[2]));
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefPRFDB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,2);
// Memory!!!
if (fDebug != 2) {
delete projprf;
}
} // Boucle object
// Plot
// No plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)){
PlotWritePRF();
}
if ((fDebug == 4) ||
(fDebug == 3)){
PlotPRFDB();
}
// Mean Statistics
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if (fWriteCoef[2]) {
WriteFitInfos(2);
}
return kTRUE;
}
//____________Functions fit Online PRF2d_______________________________________
Bool_t AliTRDCalibra::FitPRFOnline()
{
//
// Reconstruct the 1D histo (pad response function) from the vectorPRD for
// each calibration group
// Fit with a gaussian to reconstruct the sigma of the pad response function
// write the results in a tree
//
// Number of Xbins (detectors or groups of pads)
if (!InitFit(0,2)) {
return kFALSE;
}
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberEnt = 0;
// Init fCountDet and fCount
InitfCountDetAndfCount(2);
// Beginning of the loop
for (Int_t idect = fDect1[2]; idect < fDect2[2]; idect++) {
// Search if the group is in the VectorCH
Int_t place = SearchInVector(idect,2);
// Is in
TH1F *projprf = 0x0;
TString name("PRF");
name += idect;
if (place != -1) {
//Entries
fNumberEnt++;
// Variable
AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
// Retrieve
fPRFInfo = (AliTRDPInfo *) fVectorPRF->At(place);
projprf = CorrectTheError((TGraphErrors *) ConvertVectorPHisto(fPRFInfo,(const char *)name));
projprf->SetDirectory(0);
delete fPRFInfo;
}
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
UpdatefCountDetAndfCount(idect,2);
// Reconstruction of the row and pad group: rowmin, row max ...
ReconstructFitRowMinRowMax(idect,2);
// Rebin and statistic stuff
// This detector has not enough statistics or was off
if ((place == -1) ||
((place != -1) &&
(fEntriesCurrent < fMinEntries))) {
// Fill with the default values
NotEnoughStatistic(idect,2);
// Memory
if (fDebug != 2) {
delete projprf;
}
continue;
}
// Statistic of the histos fitted
AliInfo(Form("For the group number %d there are %d stats", idect,fEntriesCurrent));
fNumberFit++;
fStatisticMean += fEntriesCurrent;
// Calcul of "real" coef
if ((fDebug == 1) ||
(fDebug == 4)) {
CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
}
// Method Mean and fit
// idect is egal for fDebug = 0 and 2, only to fill the hist
FitPRF((TH1 *) projprf,(Int_t) (idect-fDect1[2]));
// Visualise the detector for fDebug 3 or 4
// Here is the reconstruction of the pad and row group is used!
if (fDebug >= 3) {
FillCoefPRFDB();
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
FillInfosFit(idect,2);
// Memory!!!
if (fDebug != 2) {
delete projprf;
}
} // Boucle object
// Plot
// No plot, 1 and 4 error plot, 3 and 4 DB plot
if ((fDebug == 1) ||
(fDebug == 4)) {
PlotWritePRF();
}
if ((fDebug == 4) ||
(fDebug == 3)) {
PlotPRFDB();
}
// Mean Statistics
if (fNumberFit > 0) {
AliInfo(Form("There is a mean statistic of: %d",(Int_t) fStatisticMean / fNumberFit));
}
else {
AliInfo("There is no fit!");
}
// Write the things!
if (fWriteCoef[2]) {
WriteFitInfos(2);
}
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibra in the code_________
Bool_t AliTRDCalibra::Init2Dhistos()
{
//
// For the offline tracking
// This function will be called in the function AliReconstruction::Run()
// Init the calibration mode (Nz, Nrphi), the 2D histograms if fHisto2d = kTRUE,
//
// DB Setting
// Get cal
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam");
return kFALSE;
}
// Some parameters
fTimeMax = cal->GetNumberOfTimeBins();
fSf = parCom->GetSamplingFrequency();
if (fRelativeScaleAuto) {
fRelativeScale = 0;
}
else {
fRelativeScale = 20;
}
// Create the 2D histos corresponding to the pad groupCalibration mode
if (fCH2dOn) {
AliInfo(Form("The pad calibration mode for the relative gain calibration: Nz %d, and Nrphi %d"
,fNz[0]
,fNrphi[0]));
// Calcul the number of Xbins
fNtotal[0] = 0;
ModePadCalibration(2,0);
ModePadFragmentation(0,2,0,0);
fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
if (fDebug == 4) {
AliInfo(Form("For the chamber 2: %d",fDetChamb2[0]));
}
fNtotal[0] += 6 * 18 * fDetChamb2[0];
ModePadCalibration(0,0);
ModePadFragmentation(0,0,0,0);
fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
if (fDebug == 4) {
AliInfo(Form("For the other chamber 0: %d",fDetChamb0[0]));
}
fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
// Create the 2D histo
if (fHisto2d) {
CreateCH2d(fNtotal[0]);
}
if (fVector2d) {
fVectorCH = new TObjArray();
fPlaCH = new TObjArray();
}
// Variable
fAmpTotal = new Float_t[TMath::Max(fDetChamb2[0],fDetChamb0[0])];
for (Int_t k = 0; k < TMath::Max(fDetChamb2[0],fDetChamb0[0]); k++) {
fAmpTotal[k] = 0.0;
}
}
if (fPH2dOn) {
AliInfo(Form("The pad calibration mode for the drift velocity calibration: Nz %d, and Nrphi %d"
,fNz[1]
,fNrphi[1]));
// Calcul the number of Xbins
fNtotal[1] = 0;
ModePadCalibration(2,1);
ModePadFragmentation(0,2,0,1);
fDetChamb2[1] = fNfragZ[1]*fNfragRphi[1];
if (fDebug == 4) {
AliInfo(Form("For the chamber 2: %d",fDetChamb2[1]));
}
fNtotal[1] += 6 * 18 * fDetChamb2[1];
ModePadCalibration(0,1);
ModePadFragmentation(0,0,0,1);
fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
if (fDebug == 4) {
AliInfo(Form("For the chamber 0: %d",fDetChamb0[1]));
}
fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
// Create the 2D histo
if (fHisto2d) {
CreatePH2d(fNtotal[1]);
}
if (fVector2d) {
fVectorPH = new TObjArray();
fPlaPH = new TObjArray();
}
// Variable
fPHPlace = new Short_t[fTimeMax];
for (Int_t k = 0; k < fTimeMax; k++) {
fPHPlace[k] = -1;
}
fPHValue = new Float_t[fTimeMax];
for (Int_t k = 0; k < fTimeMax; k++) {
fPHValue[k] = -1.0;
}
}
if (fPRF2dOn) {
AliInfo(Form("The pad calibration mode for the PRF calibration: Nz %d, and Nrphi %d"
,fNz[2]
,fNrphi[2]));
// Calcul the number of Xbins
fNtotal[2] = 0;
ModePadCalibration(2,2);
ModePadFragmentation(0,2,0,2);
fDetChamb2[2] = fNfragZ[2] * fNfragRphi[2];
if (fDebug == 4) {
AliInfo(Form("For the chamber 2: %d",fDetChamb2[2]));
}
fNtotal[2] += 6 * 18 * fDetChamb2[2];
ModePadCalibration(0,2);
ModePadFragmentation(0,0,0,2);
fDetChamb0[2] = fNfragZ[2] * fNfragRphi[2];
if (fDebug == 4) {
AliInfo(Form("For the chamber 0: %d",fDetChamb0[2]));
}
fNtotal[2] += 6 * 4 * 18 * fDetChamb0[2];
AliInfo(Form("Total number of Xbins: %d",fNtotal[2]));
// Create the 2D histo
if (fHisto2d) {
CreatePRF2d(fNtotal[2]);
}
if (fVector2d) {
fVectorPRF = new TObjArray();
fPlaPRF = new TObjArray();
}
}
return kTRUE;
}
//____________Functions for filling the histos in the code_____________________
//____________Offine tracking in the AliTRDtracker_____________________________
Bool_t AliTRDCalibra::ResetTrack()
{
//
// For the offline tracking
// This function will be called in the function
// AliTRDtracker::FollowBackPropagation() at the beginning
// Reset the parameter to know we have a new TRD track
//
fDetectorAliTRDtrack = kFALSE;
return kTRUE;
}
//____________Offline tracking in the AliTRDtracker____________________________
Bool_t AliTRDCalibra::UpdateHistograms(AliTRDcluster *cl, AliTRDtrack *t)
{
//
// For the offline tracking
// This function will be called in the function
// AliTRDtracker::FollowBackPropagation() in the loop over the clusters
// of TRD tracks
// Fill the 2D histos or the vectors with the info of the clusters at
// the end of a detectors if the track is "good"
//
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam");
return kFALSE;
}
// Get the parameter object
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
// Localisation of the detector
Int_t detector = cl->GetDetector();
Int_t chamber = GetChamber(detector);
Int_t plane = GetPlane(detector);
// Fill the infos for the previous clusters if not the same
// detector anymore or if not the same track
if (((detector != fDetectorPreviousTrack) || (!fDetectorAliTRDtrack)) &&
(fDetectorPreviousTrack != -1)) {
fNumberTrack++;
// If the same track, then look if the previous detector is in
// the same plane, if yes: not a good track
if (fDetectorAliTRDtrack &&
(GetPlane(detector) <= GetPlane(fDetectorPreviousTrack))) {
fGoodTrack = kFALSE;
}
// Fill only if the track doesn't touch a masked pad or doesn't
// appear in the middle (fGoodTrack)
if (fGoodTrack) {
// Gain calibration
if (fCH2dOn) {
FillTheInfoOfTheTrackCH();
}
// PH calibration
if (fPH2dOn) {
FillTheInfoOfTheTrackPH();
}
} // if a good track
ResetfVariables();
} // Fill at the end the charge
// Calcul the position of the detector
if (detector != fDetectorPreviousTrack) {
LocalisationDetectorXbins(detector);
}
// Reset the good track for the PRF
Bool_t good = kTRUE;
// Localisation of the cluster
Double_t pos[3] = { 0.0, 0.0, 0.0 };
pos[0] = cl->GetX();
pos[1] = cl->GetY();
pos[2] = cl->GetZ();
Int_t time = cl->GetLocalTimeBin();
// Reset the detector
fDetectorPreviousTrack = detector;
fDetectorAliTRDtrack = kTRUE;
// Position of the cluster
AliTRDpadPlane *padplane = parCom->GetPadPlane(plane,chamber);
Int_t row = padplane->GetPadRowNumber(pos[2]);
Double_t offsetz = padplane->GetPadRowOffset(row,pos[2]);
Double_t offsettilt = padplane->GetTiltOffset(offsetz);
Int_t col = padplane->GetPadColNumber(pos[1] + offsettilt,offsetz);
// See if we are not near a masked pad
if (!IsPadOn(detector,col,row)) {
good = kFALSE;
fGoodTrack = kFALSE;
}
if (col > 0) {
if (!IsPadOn(detector,col-1,row)) {
fGoodTrack = kFALSE;
good = kFALSE;
}
}
if (col < 143) {
if (!IsPadOn(detector,col+1,row)) {
fGoodTrack = kFALSE;
good = kFALSE;
}
}
// Row of the cluster and position in the pad groups
Int_t posr[3] = { 0, 0, 0 };
if ((fCH2dOn) && (fNnZ[0] != 0)) {
posr[0] = (Int_t) row / fNnZ[0];
}
if ((fPH2dOn) && (fNnZ[1] != 0)) {
posr[1] = (Int_t) row / fNnZ[1];
}
if ((fPRF2dOn) && (fNnZ[2] != 0)) {
posr[2] = (Int_t) row / fNnZ[2];
}
// Col of the cluster and position in the pad groups
Int_t posc[3] = { 0, 0, 0 };
if ((fCH2dOn) && (fNnRphi[0] != 0)) {
posc[0] = (Int_t) col / fNnRphi[0];
}
if ((fPH2dOn) && (fNnRphi[1] != 0)) {
posc[1] = (Int_t) col / fNnRphi[1];
}
if ((fPRF2dOn) && (fNnRphi[2] != 0)) {
posc[2] = (Int_t) col / fNnRphi[2];
}
// Charge in the cluster
// For the moment take the abs
Float_t q = TMath::Abs(cl->GetQ());
Short_t *signals = cl->GetSignals();
// Correction due to the track angle
Float_t correction = 1.0;
Float_t normalisation = 6.67;
if ((q >0) && (t->GetNdedx() > 0)) {
correction = t->GetClusterdQdl((t->GetNdedx() - 1)) / (q * normalisation);
}
// Fill the fAmpTotal with the charge
if (fCH2dOn) {
if (!fTraMaxPad){
fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += q * correction;
}
else {
fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += ((Float_t) signals[3]) * correction;
}
}
// Fill the fPHPlace and value
if (fPH2dOn) {
fPHPlace[time] = posc[1]*fNfragZ[1]+posr[1];
if (!fTraMaxPad) {
fPHValue[time] = q * correction;
}
else {
fPHValue[time] = ((Float_t) signals[3]) * correction;
}
}
// Fill direct the PRF
if ((fPRF2dOn) && (good)) {
Float_t yminus = 0.0;
Float_t xcenter = 0.0;
Float_t ycenter = 0.0;
Float_t ymax = 0.0;
Bool_t echec = kFALSE;
if ((cl->From3pad()) && (!cl->IsUsed())) {
// Center 3 balanced
if ((((Float_t) signals[3]) > fThresholdClusterPRF2) &&
(((Float_t) signals[2]) > fThresholdClusterPRF2) &&
(((Float_t) signals[4]) > fThresholdClusterPRF2) &&
(((Float_t) signals[1]) < fThresholdClusterPRF1) &&
(((Float_t) signals[5]) < fThresholdClusterPRF1) &&
((((Float_t) signals[2])*((Float_t) signals[4])/(((Float_t) signals[3])*((Float_t) signals[3]))) < 0.06)) {
// Col correspond to signals[3]
if (fCenterOfflineCluster) {
xcenter = cl->GetCenter();
}
else {
// Security of the denomiateur is 0
if ((((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3]))) /
((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))) != 1.0) {
xcenter = 0.5 * (TMath::Log((Float_t) (((Float_t) signals[4]) / ((Float_t) signals[2]))))
/ (TMath::Log(((Float_t) (((Float_t) signals[3]) * ((Float_t) signals[3])))
/ ((Float_t) (((Float_t) signals[2]) * ((Float_t) signals[4])))));
}
else {
xcenter = -100.0;
}
}
if ((xcenter > -0.5) && (xcenter < 0.5)) {
ycenter = (Float_t) (((Float_t) signals[3])
/ (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
yminus = (Float_t) (((Float_t) signals[2])
/ (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
ymax = (Float_t) (((Float_t) signals[4])
/ (((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4]))));
if ((TMath::Abs(((Float_t) signals[2]) + ((Float_t) signals[3]) + (((Float_t) signals[4])) - q) < 10.0)) {
echec = kTRUE;
}
}
}
// Fill only if it is in the drift region!
if ((((Float_t) (((Float_t) time) / fSf)) > 0.3) && (echec)) {
if (fHisto2d) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),xcenter,ycenter);
if (xcenter < 0.0) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),-(xcenter+1.0),yminus);
}
if (xcenter > 0.0) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),1.0-xcenter,ymax);
}
}
if (fVector2d) {
UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],xcenter,ycenter);
if (xcenter < 0.0) {
UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],-(xcenter+1.0),yminus);
}
if (xcenter > 0.0) {
UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],1.0-xcenter,ymax);
}
}
} // If in the drift region
} // Cluster isole
} // PRF2dOn
return kTRUE;
}
//____________Online trackling in AliTRDtrigger________________________________
Bool_t AliTRDCalibra::UpdateHistogramcm(AliTRDmcmTracklet *trk)
{
//
// For the tracking
// This function will be called in the function AliTRDtrigger::TestTracklet
// before applying the pt cut on the tracklets
// Fill the infos for the tracklets fTrkTest if the tracklets is "good"
//
// Localisation of the Xbins involved
Int_t idect = trk->GetDetector();
LocalisationDetectorXbins(idect);
// Get the parameter object
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
// Reset
ResetfVariables();
// Row of the tracklet and position in the pad groups
Int_t row = trk->GetRow();
Int_t posr[3] = { 0, 0, 0 };
if ((fCH2dOn) && (fNnZ[0] != 0)) {
posr[0] = (Int_t) row / fNnZ[0];
}
if ((fPH2dOn) && (fNnZ[1] != 0)) {
posr[1] = (Int_t) row / fNnZ[1];
}
if ((fPRF2dOn) && (fNnZ[2] != 0)) {
posr[2] = (Int_t) row / fNnZ[2];
}
// Eventuelle correction due to track angle in z direction
Float_t correction = 1.0;
if (fMcmCorrectAngle) {
Float_t z = trk->GetRowz();
Float_t r = trk->GetTime0();
correction = r / TMath::Sqrt((r*r+z*z));
}
// Boucle sur les clusters
// Condition on number of cluster: don't come from the middle of the detector
if (trk->GetNclusters() >= fNumberClusters) {
for (Int_t icl = 0; icl < trk->GetNclusters(); icl++) {
Float_t amp[3] = { 0.0, 0.0, 0.0 };
Int_t time = trk->GetClusterTime(icl);
Int_t col = trk->GetClusterCol(icl);
amp[0] = trk->GetClusterADC(icl)[0] * correction;
amp[1] = trk->GetClusterADC(icl)[1] * correction;
amp[2] = trk->GetClusterADC(icl)[2] * correction;
if ((amp[0] < 0.0) ||
(amp[1] < 0.0) ||
(amp[2] < 0.0)) {
continue;
}
// Col of cluster and position in the pad groups
Int_t posc[3] = { 0, 0, 0 };
if ((fCH2dOn) && (fNnRphi[0] != 0)) {
posc[0] = (Int_t) col / fNnRphi[0];
}
if ((fPH2dOn) && (fNnRphi[1] != 0)) {
posc[1] = (Int_t) col / fNnRphi[1];
}
if ((fPRF2dOn) && (fNnRphi[2] != 0)) {
posc[2] = (Int_t) col / fNnRphi[2];
}
// See if we are not near a masked pad
Bool_t good = kTRUE;
if (!IsPadOn(idect,col,row)) {
fGoodTrack = kFALSE;
good = kFALSE;
}
if (col > 0) {
if (!IsPadOn(idect,col-1,row)) {
fGoodTrack = kFALSE;
good = kFALSE;
}
}
if (col < 143) {
if (!IsPadOn(idect,col+1,row)) {
fGoodTrack = kFALSE;
good = kFALSE;
}
}
// Total spectrum
if (fPH2dOn) {
fPHPlace[time] = posc[1] * fNfragZ[1] + posr[1];
}
if (!fTraMaxPad) {
if (fCH2dOn) {
fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += (Float_t) (amp[0]+amp[1]+amp[2]);
}
if (fPH2dOn) {
fPHValue[time] = (Float_t) (amp[0]+amp[1]+amp[2]);
}
}
else {
if (fCH2dOn) {
fAmpTotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += (Float_t) amp[1];
}
if (fPH2dOn) {
fPHValue[time] = amp[1];
}
}
// Fill PRF direct
if (fPRF2dOn && good) {
if ((amp[0] > fThresholdClusterPRF2) &&
(amp[1] > fThresholdClusterPRF2) &&
(amp[2] > fThresholdClusterPRF2) &&
((amp[0]*amp[2]/(amp[1]*amp[1])) < 0.06)) {
// Security of the denomiateur is 0
if ((((Float_t) (((Float_t) amp[1]) * ((Float_t) amp[1])))
/ ((Float_t) (((Float_t) amp[0]) * ((Float_t) amp[2])))) != 1.0) {
Float_t xcenter = 0.5 * (TMath::Log(amp[2] / amp[0]))
/ (TMath::Log((amp[1]*amp[1]) / (amp[0]*amp[2])));
Float_t ycenter = amp[1] / (amp[0] + amp[1] + amp[2]);
if ((xcenter > -0.5) &&
(xcenter < 0.5)) {
Float_t yminus = amp[0] / (amp[0]+amp[1]+amp[2]);
Float_t ymax = amp[2] / (amp[0]+amp[1]+amp[2]);
// Fill only if it is in the drift region!
if (((Float_t) time / fSf) > 0.3) {
if (fHisto2d) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),xcenter,ycenter);
if (xcenter < 0.0) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),-(xcenter+1.0),yminus);
}
if (xcenter > 0.0) {
fPRF2d->Fill((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]+0.5),(1.0-xcenter),ymax);
}
}
if (fVector2d) {
UpdateVectorPRF((fXbins[2]+posc[2]*fNfragZ[2]+posr[2]),xcenter,ycenter);
if (xcenter < 0.0) {
UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],-(xcenter+1.0),yminus);
}
if (xcenter > 0.0) {
UpdateVectorPRF(fXbins[2]+posc[2]*fNfragZ[2]+posr[2],(1.0-xcenter),ymax);
}
}
}
}
}
}
}
} // Boucle clusters
// Fill the charge
if (fCH2dOn && fGoodTrack) {
FillTheInfoOfTheTrackCH();
}
// PH calibration
if (fPH2dOn && fGoodTrack) {
FillTheInfoOfTheTrackPH();
}
} // Condition on number of clusters
return kTRUE;
}
//____________Functions for seeing if the pad is really okey___________________
//_____________________________________________________________________________
Bool_t AliTRDCalibra::SetModeCalibrationFromTObject(TObject *object, Int_t i)
{
//
// Set fNz[i] and fNrphi[i] of the AliTRDCalibra::Instance()
// corresponding to the given TObject
//
const char *nametitle = object->GetTitle();
// Some patterns
const Char_t *patternz0 = "Nz0";
const Char_t *patternz1 = "Nz1";
const Char_t *patternz2 = "Nz2";
const Char_t *patternz3 = "Nz3";
const Char_t *patternz4 = "Nz4";
const Char_t *patternrphi0 = "Nrphi0";
const Char_t *patternrphi1 = "Nrphi1";
const Char_t *patternrphi2 = "Nrphi2";
const Char_t *patternrphi3 = "Nrphi3";
const Char_t *patternrphi4 = "Nrphi4";
const Char_t *patternrphi5 = "Nrphi5";
const Char_t *patternrphi6 = "Nrphi6";
UShort_t testz = 0;
UShort_t testrphi = 0;
// Nz mode
if (strstr(nametitle,patternz0)) {
testz++;
fNz[i] = 0;
}
if (strstr(nametitle,patternz1)) {
testz++;
fNz[i] = 1;
}
if (strstr(nametitle,patternz2)) {
testz++;
fNz[i] = 2;
}
if (strstr(nametitle,patternz3)) {
testz++;
fNz[i] = 3;
}
if (strstr(nametitle,patternz4)) {
testz++;
fNz[i] = 4;
}
// Nrphi mode
if (strstr(nametitle,patternrphi0)) {
testrphi++;
fNrphi[i] = 0;
}
if (strstr(nametitle,patternrphi1)) {
testrphi++;
fNrphi[i] = 1;
}
if (strstr(nametitle,patternrphi2)) {
testrphi++;
fNrphi[i] = 2;
}
if (strstr(nametitle,patternrphi3)) {
testrphi++;
fNrphi[i] = 3;
}
if (strstr(nametitle,patternrphi4)) {
testrphi++;
fNrphi[i] = 4;
}
if (strstr(nametitle,patternrphi5)) {
testrphi++;
fNrphi[i] = 5;
}
if (strstr(nametitle,patternrphi6)) {
testrphi++;
fNrphi[i] = 6;
}
// Look if all is okey
if ((testz == 1) &&
(testrphi == 1)) {
return kTRUE;
}
else {
fNrphi[i] = 0;
fNz[i] = 0;
return kFALSE;
}
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::IsPadOn(Int_t detector, Int_t col, Int_t row) const
{
//
// Look in the choosen database if the pad is On.
// If no the track will be "not good"
//
// Get the parameter object
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
if (!cal->IsChamberInstalled(detector) ||
cal->IsChamberMasked(detector) ||
cal->IsPadMasked(detector,col,row)) {
return kFALSE;
}
else {
return kTRUE;
}
}
//____________Functions for plotting the 2D____________________________________
//_____________________________________________________________________________
void AliTRDCalibra::Plot2d()
{
//
// Plot the 2D histos
//
if (fCH2dOn) {
PlotCH2d();
}
if (fPH2dOn) {
PlotPH2d();
}
if (fPRF2dOn) {
PlotPRF2d();
}
}
//____________Writing the 2D___________________________________________________
//_____________________________________________________________________________
Bool_t AliTRDCalibra::Write2d()
{
//
// Write the 2D histograms or the vectors converted in trees in the file
// "TRD.calibration.root"
//
TFile *fout = TFile::Open(fWriteName,"RECREATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
return kFALSE;
}
AliInfo(Form("Numbertrack: %d Numberusedch[0]: %d, Numberusedch[1]: %d Numberusedph[0]: %d, Numberusedph[1]: %d"
,fNumberTrack
,fNumberUsedCh[0]
,fNumberUsedCh[1]
,fNumberUsedPh[0]
,fNumberUsedPh[1]));
TStopwatch stopwatch;
stopwatch.Start();
AliInfo("Write2d");
if ((fCH2dOn ) && (fWrite[0])) {
if (fHisto2d) {
fout->WriteTObject(fCH2d);
}
if (fVector2d) {
TString name("Nz");
name += fNz[0];
name += "Nrphi";
name += fNrphi[0];
TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const char *) name);
fout->WriteTObject(treeCH2d);
}
}
if ((fPH2dOn ) && (fWrite[1])) {
if (fHisto2d) {
fout->WriteTObject(fPH2d);
}
if (fVector2d) {
TString name("Nz");
name += fNz[1];
name += "Nrphi";
name += fNrphi[1];
TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const char *) name);
fout->WriteTObject(treePH2d);
}
}
if ((fPRF2dOn ) && (fWrite[2])) {
if (fHisto2d) {
fout->WriteTObject(fPRF2d);
}
if (fVector2d) {
TString name("Nz");
name += fNz[2];
name += "Nrphi";
name += fNrphi[2];
TTree *treePRF2d = ConvertVectorPTreeHisto(fVectorPRF,fPlaPRF,"treePRF2d",(const char *) name);
fout->WriteTObject(treePRF2d);
}
}
fout->Close();
AliInfo(Form("Execution time Write2d: R:%.2fs C:%.2fs"
,stopwatch.RealTime(),stopwatch.CpuTime()));
return kTRUE;
}
//_____________________________________________________________________________
AliTRDCalDet *AliTRDCalibra::CreateDetObjectTree(TTree *tree, Int_t i)
{
//
// It creates the AliTRDCalDet object from the tree of the coefficient
// for the calibration i (i != 2)
// It takes the mean value of the coefficients per detector
// This object has to be written in the database
//
// Create the DetObject
AliTRDCalDet *object = 0x0;
if (i == 0) {
object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
}
if (i == 1) {
object = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
}
else {
object = new AliTRDCalDet("ChamberT0","T0 (detector value)");
}
// Read the Tree
Int_t detector = -1;
Float_t values[2304];
tree->SetBranchAddress("detector",&detector);
if (i == 0) {
tree->SetBranchAddress("gainPad",values);
}
if (i == 1) {
tree->SetBranchAddress("vdrift" ,values);
}
if (i == 3) {
tree->SetBranchAddress("t0" ,values);
}
// For calculating the mean
Float_t mean = 0.0;
Int_t nto = 0;
Int_t numberofentries = tree->GetEntries();
if (numberofentries != 540) {
AliInfo("The tree is not complete");
}
for (Int_t det = 0; det < numberofentries; ++det) {
tree->GetEntry(det);
if (GetChamber(detector) == 2) {
nto = 1728;
}
else {
nto = 2304;
}
mean = 0.0;
if(i != 3){
for (Int_t k = 0; k < nto; k++) {
mean += TMath::Abs(values[k]) / nto;
}
}
else {
for (Int_t k = 0; k < nto; k++) {
if(k == 0) mean = values[k];
if(mean > values[k]) mean = values[k];
}
}
object->SetValue(detector,mean);
}
return object;
}
//_____________________________________________________________________________
TObject *AliTRDCalibra::CreatePadObjectTree(TTree *tree, Int_t i
, AliTRDCalDet *detobject)
{
//
// It Creates the AliTRDCalPad object from the tree of the
// coefficient for the calibration i (i != 2)
// You need first to create the object for the detectors,
// where the mean value is put.
// This object has to be written in the database
//
// Create the DetObject
AliTRDCalPad *object = 0x0;
if (i == 0) {
object = new AliTRDCalPad("GainFactor","GainFactor (local variations)");
}
if (i == 1) {
object = new AliTRDCalPad("LocalVdrift","TRD drift velocities (local variations)");
}
else {
object = new AliTRDCalPad("LocalT0","T0 (local variations)");
}
// Read the Tree
Int_t detector = -1;
Float_t values[2304];
tree->SetBranchAddress("detector",&detector);
if (i == 0) {
tree->SetBranchAddress("gainPad",values);
}
if (i == 1) {
tree->SetBranchAddress("vdrift" ,values);
}
if (i == 3) {
tree->SetBranchAddress("t0" ,values);
}
// Variables
Float_t mean = 0.0;
Int_t numberofentries = tree->GetEntries();
if (numberofentries != 540) {
AliInfo("The tree is not complete");
}
for (Int_t det = 0; det < numberofentries; ++det) {
tree->GetEntry(det);
AliTRDCalROC *calROC = object->GetCalROC(detector);
mean = detobject->GetValue(detector);
if ((mean == 0) && (i != 3)) {
continue;
}
Int_t rowMax = calROC->GetNrows();
Int_t colMax = calROC->GetNcols();
for (Int_t row = 0; row < rowMax; ++row) {
for (Int_t col = 0; col < colMax; ++col) {
if(i != 3) calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)])/mean);
else calROC->SetValue(col,row,values[(Int_t) (col*rowMax+row)]-mean);
} // Col
} // Row
}
return object;
}
//_____________________________________________________________________________
TObject *AliTRDCalibra::CreatePadObjectTree(TTree *tree)
{
//
// It Creates the AliTRDCalPad object from the tree of the
// coefficient for the calibration PRF (i = 2)
// This object has to be written in the database
//
// Create the DetObject
AliTRDCalPad *object = new AliTRDCalPad("PRFWidth","PRFWidth");
// Read the Tree
Int_t detector = -1;
Float_t values[2304];
tree->SetBranchAddress("detector",&detector);
tree->SetBranchAddress("width" ,values);
// Variables
Int_t numberofentries = tree->GetEntries();
if (numberofentries != 540) {
AliInfo("The tree is not complete");
}
for (Int_t det = 0; det < numberofentries; ++det) {
tree->GetEntry(det);
AliTRDCalROC *calROC = object->GetCalROC(detector);
Int_t rowMax = calROC->GetNrows();
Int_t colMax = calROC->GetNcols();
for (Int_t row = 0; row < rowMax; ++row) {
for (Int_t col = 0; col < colMax; ++col) {
calROC->SetValue(col,row,TMath::Abs(values[(Int_t) (col*rowMax+row)]));
} // Col
} // Row
}
return object;
}
//_____________________________________________________________________________
void AliTRDCalibra::SetRelativeScale(Float_t RelativeScale)
{
//
// Set the factor that will divide the deposited charge
// to fit in the histo range [0,300]
//
if (RelativeScale > 0.0) {
fRelativeScale = RelativeScale;
}
else {
AliInfo("RelativeScale must be strict positif!");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetNz(Int_t i, Short_t Nz)
{
//
// Set the mode of calibration group in the z direction for the parameter i
//
if ((Nz >= 0) &&
(Nz < 5)) {
fNz[i] = Nz;
}
else {
AliInfo("You have to choose between 0 and 4");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetNrphi(Int_t i, Short_t Nrphi)
{
//
// Set the mode of calibration group in the rphi direction for the parameter i
//
if ((Nrphi >= 0) &&
(Nrphi < 7)) {
fNrphi[i] = Nrphi;
}
else {
AliInfo("You have to choose between 0 and 6");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetPeriodeFitPH(Int_t periodeFitPH)
{
//
// Set FitPH if 1 then each detector will be fitted
//
if (periodeFitPH > 0) {
fFitPHPeriode = periodeFitPH;
}
else {
AliInfo("periodeFitPH must be higher than 0!");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetBeginFitCharge(Float_t beginFitCharge)
{
//
// The fit of the deposited charge distribution begins at
// histo->Mean()/beginFitCharge
// You can here set beginFitCharge
//
if (beginFitCharge > 0) {
fBeginFitCharge = beginFitCharge;
}
else {
AliInfo("beginFitCharge must be strict positif!");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetT0Shift(Float_t t0Shift)
{
//
// The t0 calculated with the maximum positif slope is shift from t0Shift
// You can here set t0Shift
//
if (t0Shift > 0) {
fT0Shift = t0Shift;
}
else {
AliInfo("t0Shift must be strict positif!");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetRangeFitPRF(Float_t rangeFitPRF)
{
//
// The fit of the PRF is from -rangeFitPRF to rangeFitPRF
// You can here set rangeFitPRF
//
if ((rangeFitPRF > 0) &&
(rangeFitPRF <= 1.0)) {
fRangeFitPRF = rangeFitPRF;
}
else {
AliInfo("rangeFitPRF must be between 0 and 1.0");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::SetRebin(Short_t rebin)
{
//
// Rebin with rebin time less bins the Ch histo
// You can set here rebin that should divide the number of bins of CH histo
//
if (rebin > 0) {
fRebin = rebin;
AliInfo("You have to be sure that fRebin divides fNumberBinCharge used!");
}
else {
AliInfo("You have to choose a positiv value!");
}
}
//_____________________________________________________________________________
TTree *AliTRDCalibra::Sum2Trees(const Char_t *filename1
, const Char_t *filename2
, const Char_t *variablecali)
{
//
// It returns the sum of two trees with the name variablecali
// in the files filenam1 and filename2 equivalent of merging two 2D histos
// The name of the resulting tree is the same as the two input trees
// variablecali can be treeCH2d, treePH2d or treePRF2d
//
// Variables
TChain *treeChain = new TChain(variablecali);
TObjArray *vectorplace = new TObjArray();
TObjArray *where = new TObjArray();
// First tree
// Take the tree
TFile *file1 = new TFile(filename1,"READ");
TTree *tree1 = (TTree *) file1->Get(variablecali);
gDirectory = gROOT;
// Take the places
vectorplace = ConvertTreeVector(tree1);
// Say where it is in tree 1
for (Int_t jui = 0; jui < (Int_t) vectorplace->GetEntriesFast(); jui++) {
AliTRDPlace *placejui = new AliTRDPlace();
placejui->SetPlace(jui);
TObjArray *chainplace = new TObjArray();
chainplace->Add((TObject *) placejui);
where->Add((TObject *) chainplace);
}
// Add to the chain
treeChain->Add(filename1);
delete file1;
// Second tree
// Take the tree
TFile *file2 = new TFile(filename2,"READ");
TTree *tree2 = (TTree *) file2->Get(variablecali);
gDirectory = gROOT;
// Take the places
TObjArray *vector2 = ConvertTreeVector(tree2);
Int_t j = treeChain->GetEntries();
for (Int_t jui = 0; jui < (Int_t) vector2->GetEntriesFast(); jui++) {
// Search if already found
Int_t place = SearchInTreeVector(vectorplace,((AliTRDPlace *) vector2->At(jui))->GetPlace());
// Create a new element in the two std vectors
if (place == -1) {
AliTRDPlace *placejjui = new AliTRDPlace();
placejjui->SetPlace((j+jui));
TObjArray *chainplace = new TObjArray();
chainplace->Add((TObject *) placejjui);
vectorplace->Add((TObject *) (vector2->At(jui)));
where->Add((TObject *) chainplace);
}
// Update the element at the place "place" in the std vector whereinthechain
else {
AliTRDPlace *placejjui = new AliTRDPlace();
placejjui->SetPlace((j+jui));
TObjArray *chainplace = ((TObjArray *) where->At(place));
chainplace->Add((TObject *) placejjui);
where->AddAt((TObject *) chainplace,place);
}
}
// Add to the Chain
treeChain->Add(filename2);
delete file2;
// Take care of the profile
const Char_t *pattern = "P";
TTree *tree = 0x0;
if (!strstr(variablecali,pattern)) {
// Ready to read the chain
TH1F *his = 0x0;
treeChain->SetBranchAddress("histo",&his);
// Initialise the final tree
Int_t group = -1;
TH1F *histsum = 0x0;
tree = new TTree(variablecali,variablecali);
tree->Branch("groupnumber",&group,"groupnumber/I");
tree->Branch("histo","TH1F",&histsum,32000,0);
// Init histsum
if (treeChain->GetEntries() < 1) {
return tree1;
}
for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
TObjArray *chainplace = ((TObjArray *) where->At(h));
treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
//Init for the first time
if (h == 0) {
histsum = new TH1F("","",his->GetXaxis()->GetNbins()
,his->GetXaxis()->GetBinLowEdge(1)
,his->GetXaxis()->GetBinUpEdge(his->GetXaxis()->GetNbins()));
histsum->Sumw2();
}
// Reset for each new group
histsum->SetEntries(0.0);
for (Int_t l = 0; l <= histsum->GetXaxis()->GetNbins(); l++) {
histsum->SetBinContent(l,0.0);
histsum->SetBinError(l,0.0);
}
histsum->Add(his,1);
if ((Int_t) chainplace->GetEntriesFast() > 1) {
for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
histsum->Add(his,1);
}
}
tree->Fill();
}
}
else {
// Ready to read the chain
TGraphErrors *his = 0x0;
treeChain->SetBranchAddress("histo",&his);
// Initialise the final tree
Int_t group = -1;
TGraphErrors *histsum = 0x0;
Double_t *xref = 0x0;
tree = new TTree(variablecali,variablecali);
tree->Branch("groupnumber",&group,"groupnumber/I");
tree->Branch("histo","TGraphErrors",&histsum,32000,0);
// Init histsum
if (treeChain->GetEntries() < 1) {
return tree1;
}
for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
TObjArray *chainplace = ((TObjArray *) where->At(h));
treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
//Init or reset for a new group
Int_t nbins = his->GetN();
Double_t *x;
x = new Double_t[nbins];
xref = his->GetX();
Double_t *ex;
ex = new Double_t[nbins];
Double_t *y;
y = new Double_t[nbins];
Double_t *ey;
ey = new Double_t[nbins];
for (Int_t lo = 0; lo < nbins; lo++) {
x[lo] = xref[lo];
ex[lo] = 0.0;
y[lo] = 0.0;
ey[lo] = 0.0;
}
delete histsum;
histsum = new TGraphErrors(nbins,x,y,ex,ey);
// Add the first
histsum = AddProfiles(his,histsum);
if ((Int_t) chainplace->GetEntriesFast() > 1) {
for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
histsum = AddProfiles(his,histsum);
}
}
tree->Fill();
}
}
return tree;
}
//____________Function fill 2D for the moment out of the code__________________
//____________Function fill 2D all objects from digits_________________________
Bool_t AliTRDCalibra::Create2DDiSimOnline(Int_t iev1, Int_t iev2)
{
//
// Only for simulations, after the simulation, create the 2D histos
// from the digits stored in the file "TRD.Digits.root"
// Only for CH and PH
//
const Int_t kNplan = 6;
const Int_t kNcham = 5;
// RunLoader and so on
if (gAlice) {
delete gAlice->GetRunLoader();
delete gAlice;
gAlice = 0;
}
AliRunLoader *rl = AliRunLoader::Open("galice.root");
if (!rl) {
return kFALSE;
}
rl->LoadgAlice();
gAlice = rl->GetAliRun();
if (!gAlice) {
return kFALSE;
}
// Import the Trees for the event nEvent in the file
rl->LoadKinematics();
rl->GetEvent(0);
rl->LoadHeader();
AliLoader *loader = rl->GetLoader("TRDLoader");
if (!loader) {
AliInfo("No TRDLLoader found!");
return kFALSE;
}
// Get the pointer to the TRD detector
AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
if (!trd) {
AliInfo("No TRD detector found");
return kFALSE;
}
// Get the pointer to the geometry object
AliTRDgeometry *geo;
if (trd) {
geo = trd->GetGeometry();
}
else {
AliInfo("No TRD geometry found");
return kFALSE;
}
// DB Setting
AliCDBManager *man = AliCDBManager::Instance();
if (!man) {
AliInfo("Could not get CDB Manager");
return kFALSE;
}
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam");
return kFALSE;
}
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
// Some parameters
fTimeMax = cal->GetNumberOfTimeBins();
fSf = (Float_t) parCom->GetSamplingFrequency();
if (fRelativeScaleAuto) {
fRelativeScale = 0;
}
else {
if (fRelativeScale <= 0.0) {
AliInfo("You have to set the relativescale factor per hand!");
return kFALSE;
}
}
// Create the 2D histos corresponding to the pad group calibration mode
if (fCH2dOn) {
AliInfo(Form("We will fill the CH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
,fNz[0]
,fNrphi[0]));
// Calcul the number of Xbins
fNtotal[0] = 0;
ModePadCalibration(2,0);
ModePadFragmentation(0,2,0,0);
fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
fNtotal[0] += 6 * 18 * fDetChamb2[0];
ModePadCalibration(0,0);
ModePadFragmentation(0,0,0,0);
fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
// Create the 2D histo
if (fHisto2d) {
CreateCH2d(fNtotal[0]);
}
if (fVector2d) {
fVectorCH = new TObjArray();
fPlaCH = new TObjArray();
}
}
if (fPH2dOn) {
AliInfo(Form("We will fill the PH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
,fNz[1]
,fNrphi[1]));
// Calcul the number of Xbins
fNtotal[1] = 0;
ModePadCalibration(2,1);
ModePadFragmentation(0,2,0,1);
fDetChamb2[1] = fNfragZ[1] * fNfragRphi[1];
fNtotal[1] += 6 * 18 * fDetChamb2[1];
ModePadCalibration(0,1);
ModePadFragmentation(0,0,0,1);
fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
// Create the 2D histo
if (fHisto2d) {
CreatePH2d(fNtotal[1]);
}
if (fVector2d) {
fVectorPH = new TObjArray();
fPlaPH = new TObjArray();
}
}
loader->LoadDigits();
AliInfo("LoadDigits ");
AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
//iev2 egal to the max if 0
if (iev2 == 0) {
iev2 = rl->GetNumberOfEvents();
AliInfo(Form("Total number of events: %d",iev2));
}
// Loop on event
for (Int_t ievent = iev1; ievent < iev2; ievent++) {
AliInfo(Form("Process event %d",ievent));
rl->GetEvent(ievent);
if (!loader->TreeD()) {
AliInfo("loader Loading Digits ... ");
loader->LoadDigits();
}
digitsManager->ReadDigits(loader->TreeD());
AliInfo("digitsManager Read Digits Done");
// Read the digits from the file
if (!(digitsManager->ReadDigits(loader->TreeD()))) {
return kFALSE;
}
// Loop on detector
for (Int_t iSect = 0; iSect < 18; iSect++) {
for (Int_t iChamb = 0; iChamb < kNcham; iChamb++) {
for (Int_t iPlane = 0; iPlane < kNplan; iPlane++) {
// A little geometry:
Int_t iDet = geo->GetDetector(iPlane,iChamb,iSect);
Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
Int_t colMax = parCom->GetColMax(iPlane);
// Variables for the group
LocalisationDetectorXbins(iDet);
// In the cas of charge
Float_t *amptotal;
amptotal = new Float_t[fNfragRphi[0]*fNfragZ[0]];
if (fCH2dOn) {
for (Int_t k = 0; k < fNfragRphi[0]*fNfragZ[0]; k++) {
amptotal[k] = 0.0;
}
}
// Loop through the detector pixel
for (Int_t time = 0; time < fTimeMax; time++) {
for (Int_t col = 0; col < colMax; col++) {
for (Int_t row = 0; row < rowMax; row++) {
// Amplitude and position in pad group
AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
Int_t amp = digit->GetAmp();
Int_t posr[2] = {0,0};
Int_t posc[2] = {0,0};
if ((fCH2dOn) &&
(fNnZ[0] != 0)) {
posr[0] = (Int_t) row / fNnZ[0];
}
if ((fCH2dOn) &&
(fNnRphi[0] != 0)) {
posc[0] = (Int_t) col / fNnRphi[0];
}
if ((fPH2dOn) &&
(fNnZ[1] != 0)) {
posr[1] = (Int_t) row / fNnZ[1];
}
if ((fPH2dOn) &&
(fNnRphi[1] != 0)) {
posc[1] = (Int_t) col / fNnRphi[1];
}
// Total spectrum
if (fCH2dOn) {
if (amp < fThresholdDigit) {
amp = 0;
}
amptotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += amp;
}
if (fPH2dOn) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+posc[1]*fNfragZ[1]+posr[1])+0.5,(Float_t) time/fSf,(Double_t) amp);
}
if (fVector2d) {
UpdateVectorPH((fXbins[1]+posc[1]*fNfragZ[1]+posr[1]),time,(Double_t) amp);
}
}
// Memory stuff
delete digit;
} // Boucle row
} // Boucle col
} // Boucle time
if (fCH2dOn) {
// If automatic scale
if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
// Take only the one zone track
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if ((fCountRelativeScale < 100) && (amptotal[k] > 2.0)) {
fRelativeScale += amptotal[k]*0.014*0.01;
fCountRelativeScale++;
}
}
}
// We fill the CH2d after having scale with the first 100
if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
// Case of
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if (fHisto2d &&
(amptotal[k] > 0.0)) {
fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
}
if (fVector2d &&
(amptotal[k] > 0.0)) {
UpdateVectorCH(fXbins[0]+k ,amptotal[k]/fRelativeScale);
}
}
}
// No relative salce
if (!fRelativeScaleAuto) {
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if (fHisto2d &&
(amptotal[k] > 0.0)) {
fCH2d->Fill(fXbins[0]+k+0.5, amptotal[k]/fRelativeScale);
}
if (fVector2d &&
(amptotal[k] > 0.0)) {
UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
}
}
}
}
delete amptotal;
} // Boucle chamber
} // Boucle plane
} // Boucle sect
loader->UnloadDigits();
} // Boucle event
if (fDebug == 1) {
if (fPH2dOn && fHisto2d) {
PlotPH2d();
}
if (fCH2dOn && fHisto2d) {
PlotCH2d();
}
}
if (fWrite[0] || fWrite[1]) {
TFile *fout = TFile::Open(fWriteName,"RECREATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("WriteTObject(fCH2d);
}
if (fPH2dOn && fHisto2d && fWrite[1]) {
fout->WriteTObject(fPH2d);
}
if (fVector2d && fCH2dOn && fWrite[0]) {
TString name("Nz");
name += fNz[0];
name += "Nrphi";
name += fNrphi[0];
TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const char *) name);
fout->WriteTObject(treeCH2d);
}
if (fVector2d && fPH2dOn && fWrite[1]) {
TString name("Nz");
name += fNz[1];
name += "Nrphi";
name += fNrphi[1];
TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const char *) name);
fout->WriteTObject(treePH2d);
}
fout->Close();
}
return kTRUE;
}
//____________Function fill 2D all objects from Raw Data_______________________
Bool_t AliTRDCalibra::Create2DRaDaOnline(Int_t iev1, Int_t iev2)
{
//
// After having written the RAW DATA in the current directory, create the
// 2D histos from these RAW DATA
// Only for CH and PH
//
const Int_t kNplan = 6;
const Int_t kNcham = 5;
TString dirname(".");
// DB Setting
AliCDBManager *man = AliCDBManager::Instance();
if (!man) {
AliInfo("Could not get CDB Manager");
return kFALSE;
}
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam");
return kFALSE;
}
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
// Some parameters
fTimeMax = cal->GetNumberOfTimeBins();
fSf = (Float_t) parCom->GetSamplingFrequency();
if (fRelativeScaleAuto) {
fRelativeScale = 0;
}
else {
if (fRelativeScale <= 0.0) {
AliInfo("You have to set the relativescale factor per hand!");
return kFALSE;
}
}
// Create the 2D histo corresponding to the pad group calibration mode
if (fCH2dOn) {
AliInfo(Form("We will fill the CH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
,fNz[0]
,fNrphi[0]));
// Calcul the number of Xbins
fNtotal[0] = 0;
ModePadCalibration(2,0);
ModePadFragmentation(0,2,0,0);
fDetChamb2[0] = fNfragZ[0] * fNfragRphi[0];
fNtotal[0] += 6 * 18 * fDetChamb2[0];
ModePadCalibration(0,0);
ModePadFragmentation(0,0,0,0);
fDetChamb0[0] = fNfragZ[0] * fNfragRphi[0];
fNtotal[0] += 6 * 4 * 18 * fDetChamb0[0];
AliInfo(Form("Total number of Xbins: %d",fNtotal[0]));
// Create the 2D histo
if (fHisto2d) {
CreateCH2d(fNtotal[0]);
}
if (fVector2d) {
fVectorCH = new TObjArray();
fPlaCH = new TObjArray();
}
}
if(fPH2dOn) {
AliInfo(Form("We will fill the PH2d histo with the pad calibration mode: Nz %d, and Nrphi %d"
,fNz[1]
,fNrphi[1]));
// Calcul the number of Xbins
fNtotal[1] = 0;
ModePadCalibration(2,1);
ModePadFragmentation(0,2,0,1);
fDetChamb2[1] = fNfragZ[1] * fNfragRphi[1];
fNtotal[1] += 6 * 18 * fDetChamb2[1];
ModePadCalibration(0,1);
ModePadFragmentation(0,0,0,1);
fDetChamb0[1] = fNfragZ[1] * fNfragRphi[1];
fNtotal[1] += 6 * 4 * 18 * fDetChamb0[1];
AliInfo(Form("Total number of Xbins: %d",fNtotal[1]));
// Create the 2D histo
if (fHisto2d) {
CreatePH2d(fNtotal[1]);
}
if (fVector2d){
fVectorPH = new TObjArray();
fPlaPH = new TObjArray();
}
}
AliTRDrawData *rawdata = new AliTRDrawData();
AliInfo("AliTRDrawData object created ");
// Loop on events
for (Int_t ievent = iev1; ievent < iev2; ievent++) {
// AliRawReaderFile
AliRawReaderFile *readerfile = new AliRawReaderFile(dirname,ievent);
if (!readerfile) {
AliInfo("No readerfile found!");
return kFALSE;
}
AliTRDdigitsManager *digitsManager = rawdata->Raw2Digits((AliRawReader *) readerfile);
if (!digitsManager) {
AliInfo("No DigitsManager done!");
return kFALSE;
}
// Loop on detectors
for (Int_t iSect = 0; iSect < 18; iSect++) {
for (Int_t iPlane = 0; iPlane < kNplan; iPlane++) {
for (Int_t iChamb = 0; iChamb < kNcham; iChamb++) {
// A little geometry:
Int_t iDet = AliTRDgeometry::GetDetector(iPlane,iChamb,iSect);
Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
Int_t colMax = parCom->GetColMax(iPlane);
// Variables for the group
LocalisationDetectorXbins(iDet);
// In the cas of charge
Float_t *amptotal;
amptotal = new Float_t[fNfragRphi[0]*fNfragZ[0]];
if(fCH2dOn) {
for (Int_t k = 0; k < fNfragRphi[0]*fNfragZ[0]; k++) {
amptotal[k] = 0.0;
}
}
// Loop through the detector pixel
for (Int_t time = 0; time < fTimeMax; time++) {
for (Int_t col = 0; col < colMax; col++) {
for (Int_t row = 0; row < rowMax; row++) {
// Amplitude and position of the digit
AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
Int_t amp = digit->GetAmp();
Int_t posr[2] = { 0, 0 };
Int_t posc[2] = { 0, 0 };
if ((fCH2dOn) &&
(fNnZ[0] != 0)) {
posr[0] = (Int_t) row / fNnZ[0];
}
if ((fCH2dOn) &&
(fNnRphi[0] != 0)) {
posc[0] = (Int_t) col / fNnRphi[0];
}
if ((fPH2dOn) &&
(fNnZ[1] != 0)) {
posr[1] = (Int_t) row / fNnZ[1];
}
if ((fPH2dOn) &&
(fNnRphi[1] != 0)) {
posc[1] = (Int_t) col / fNnRphi[1];
}
// Total spectrum
if (fCH2dOn) {
if (amp < fThresholdDigit) {
amp = 0;
}
amptotal[(Int_t) (posc[0]*fNfragZ[0]+posr[0])] += amp;
}
if (fPH2dOn ) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+posc[1]*fNfragZ[1]+posr[1])+0.5,(Float_t)time/fSf,amp);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+posc[1]*fNfragZ[1]+posr[1],time,amp);
}
}
delete digit;
} // Boucle row
} // Boucle col
} // Boucle time
if (fCH2dOn) {
// If automatic scale
if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
// Take only the one zone track
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if ((fCountRelativeScale < 100) && (amptotal[k] > 2.0)) {
fRelativeScale += amptotal[k] * 0.014 * 0.01;
fCountRelativeScale++;
}
}
}
// We fill the CH2d after having scale with the first 100
if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
// Case of
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if (fHisto2d && (amptotal[k] > 0.0)) {
fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
}
if (fVector2d && (amptotal[k] > 0.0)) {
UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
}
}
}
// No relative salce
if (!fRelativeScaleAuto) {
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if (fHisto2d &&
(amptotal[k] > 0.0)) {
fCH2d->Fill(fXbins[0]+k+0.5,amptotal[k]/fRelativeScale);
}
if (fVector2d &&
(amptotal[k] > 0.0)) {
UpdateVectorCH(fXbins[0]+k, amptotal[k]/fRelativeScale);
}
}
}
}
delete amptotal;
} // Boucle chamber
} // Boucle plane
} // Boucle sect
delete digitsManager;
delete readerfile;
} // Boucle event
if (fDebug == 1) {
if (fPH2dOn && fHisto2d) {
PlotPH2d();
}
if (fCH2dOn && fHisto2d) {
PlotCH2d();
}
}
if (fWrite[0] || fWrite[1]) {
TFile *fout = TFile::Open(fWriteName,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("WriteTObject(fCH2d);
}
if (fPH2dOn && fHisto2d && fWrite[1]) {
fout->WriteTObject(fPH2d);
}
if (fVector2d && fCH2dOn && fWrite[0]) {
TString name("Nz");
name += fNz[0];
name += "Nrphi";
name += fNrphi[0];
TTree *treeCH2d = ConvertVectorCTTreeHisto(fVectorCH,fPlaCH,"treeCH2d",(const Char_t *) name);
fout->WriteTObject(treeCH2d);
}
if (fVector2d && fPH2dOn && fWrite[1]) {
TString name("Nz");
name += fNz[1];
name += "Nrphi";
name += fNrphi[1];
TTree *treePH2d = ConvertVectorPTreeHisto(fVectorPH,fPlaPH,"treePH2d",(const Char_t *) name);
fout->WriteTObject(treePH2d);
}
}
return kTRUE;
}
//____________Pad Calibration Public___________________________________________
//____________Define the number of pads per group for one detector and one calibration
void AliTRDCalibra::ModePadCalibration(Int_t iChamb, Int_t i)
{
//
// Definition of the calibration mode
// from Nz and Nrphi, the number of row and col pads per calibration groups are setted
//
fNnZ[i] = 0;
fNnRphi[i] = 0;
if ((fNz[i] == 0) && (iChamb == 2)) {
fNnZ[i] = 12;
}
if ((fNz[i] == 0) && (iChamb != 2)) {
fNnZ[i] = 16;
}
if ((fNz[i] == 1) && (iChamb == 2)) {
fNnZ[i] = 6;
}
if ((fNz[i] == 1) && (iChamb != 2)) {
fNnZ[i] = 8;
}
if ((fNz[i] == 2) && (iChamb == 2)) {
fNnZ[i] = 3;
}
if ((fNz[i] == 2) && (iChamb != 2)) {
fNnZ[i] = 4;
}
if (fNz[i] == 3) {
fNnZ[i] = 2;
}
if (fNz[i] == 4) {
fNnZ[i] = 1;
}
if (fNrphi[i] == 0) {
fNnRphi[i] = 144;
}
if (fNrphi[i] == 1) {
fNnRphi[i] = 72;
}
if (fNrphi[i] == 2) {
fNnRphi[i] = 36;
}
if (fNrphi[i] == 3) {
fNnRphi[i] = 18;
}
if (fNrphi[i] == 4) {
fNnRphi[i] = 9;
}
if (fNrphi[i] == 5) {
fNnRphi[i] = 4;
}
if (fNrphi[i] == 6) {
fNnRphi[i] = 1;
}
}
//____________Define the number of pad groups in one detector for one calibration
Bool_t AliTRDCalibra::ModePadFragmentation(Int_t iPlane,Int_t iChamb, Int_t iSect, Int_t i)
{
//
// Definition of the calibration mode
// From the number of row and col pads per calibration groups the
// number of calibration groups are setted
//
fNfragZ[i] = 0;
fNfragRphi[i] = 0;
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam Manager");
return kFALSE;
}
// A little geometry:
Int_t rowMax = parCom->GetRowMax(iPlane,iChamb,iSect);
Int_t colMax = parCom->GetColMax(iPlane);
// The fragmentation
if (fNnZ[i] != 0) {
fNfragZ[i] = (Int_t) rowMax / fNnZ[i];
}
if (fNnRphi[i] != 0) {
fNfragRphi[i] = (Int_t) colMax / fNnRphi[i];
}
return kTRUE;
}
//____________Protected Functions______________________________________________
//____________Create the 2D histo to be filled online__________________________
//
//_____________________________________________________________________________
void AliTRDCalibra::CreatePRF2d(Int_t nn)
{
//
// Create the 2D histos
//
TString name("Nz");
name += fNz[2];
name += "Nrphi";
name += fNrphi[2];
fPRF2d = new TProfile2D("PRF2d",(const Char_t *) name
,nn,0,nn,fNumberBinPRF,-1.0,1.0);
fPRF2d->SetXTitle("Det/pad groups");
fPRF2d->SetYTitle("Position x/W [pad width units]");
fPRF2d->SetZTitle("Q_{i}/Q_{total}");
fPRF2d->SetStats(0);
}
//_____________________________________________________________________________
void AliTRDCalibra::CreatePH2d(Int_t nn)
{
//
// Create the 2D histos
//
TString name("Nz");
name += fNz[1];
name += "Nrphi";
name += fNrphi[1];
fPH2d = new TProfile2D("PH2d",(const Char_t *) name
,nn,0,nn,fTimeMax
,-0.5/fSf,(Float_t) (fTimeMax-0.5)/fSf);
fPH2d->SetXTitle("Det/pad groups");
fPH2d->SetYTitle("time [#mus]");
fPH2d->SetZTitle(" [a.u.]");
fPH2d->SetStats(0);
}
//_____________________________________________________________________________
void AliTRDCalibra::CreateCH2d(Int_t nn)
{
//
// Create the 2D histos
//
TString name("Nz");
name += fNz[0];
name += "Nrphi";
name += fNrphi[0];
fCH2d = new TH2I("CH2d",(const Char_t *) name
,nn,0,nn,fNumberBinCharge,0,300);
fCH2d->SetXTitle("Det/pad groups");
fCH2d->SetYTitle("charge deposit [a.u]");
fCH2d->SetZTitle("counts");
fCH2d->SetStats(0);
fCH2d->Sumw2();
}
//____________Offine tracking in the AliTRDtracker_____________________________
void AliTRDCalibra::FillTheInfoOfTheTrackCH()
{
//
// For the offline tracking or mcm tracklets
// This function will be called in the functions UpdateHistogram...
// to fill the info of a track for the relativ gain calibration
//
Int_t nb = 0; // Nombre de zones traversees
Int_t fd = -1; // Premiere zone non nulle
// See if the track goes through different zones
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
if (fAmpTotal[k] > 0.0) {
nb++;
if (nb == 1) {
fd = k;
}
}
}
// If automatic scale
if ((fCountRelativeScale < 100) && (fRelativeScaleAuto)) {
// Take only the one zone track
if (nb == 1) {
fRelativeScale += fAmpTotal[fd] * 0.014 * 0.01;
fCountRelativeScale++;
}
}
// We fill the CH2d after having scale with the first 100
if ((fCountRelativeScale >= 100) && (fRelativeScaleAuto)) {
// Case of track with only one zone
if (nb == 1) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
}
} // Case 1 zone
// Case of track with two zones
if (nb == 2) {
// Two zones voisines sinon rien!
if ((fAmpTotal[fd] > 0.0) &&
(fAmpTotal[fd+1] > 0.0)) {
// One of the two very big
if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
}
}
if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd+1]/fRelativeScale);
}
}
}
} // Case 2 zones
}
// Fill with no automatic scale
if (!fRelativeScaleAuto) {
// Case of track with only one zone
if (nb == 1) {
fNumberUsedCh[0]++;
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
}
} // Case 1 zone
// Case of track with two zones
if (nb == 2) {
// Two zones voisines sinon rien!
// Case 1
if ((fAmpTotal[fd] > 0.0) &&
(fAmpTotal[fd+1] > 0.0)) {
// One of the two very big
if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+1]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
}
fNumberUsedCh[1]++;
}
if (fAmpTotal[fd+1] > fProcent*fAmpTotal[fd]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+1.5,fAmpTotal[fd+1]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd+1,fAmpTotal[fd+1]/fRelativeScale);
}
fNumberUsedCh[1]++;
}
}
// Case 2
if (fNfragZ[0] > 1) {
if (fAmpTotal[fd] > 0.0) {
if ((fd+fNfragZ[0]) < (fNfragZ[0]*fNfragRphi[0])) {
if (fAmpTotal[fd+fNfragZ[0]] > 0.0) {
// One of the two very big
if (fAmpTotal[fd] > fProcent*fAmpTotal[fd+fNfragZ[0]]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+0.5,fAmpTotal[fd]/fRelativeScale);
}
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd,fAmpTotal[fd]/fRelativeScale);
}
fNumberUsedCh[1]++;
}
if (fAmpTotal[fd+fNfragZ[0]] > fProcent*fAmpTotal[fd]) {
if (fHisto2d) {
fCH2d->Fill(fXbins[0]+fd+fNfragZ[0]+0.5,fAmpTotal[fd+fNfragZ[0]]/fRelativeScale);
}
fNumberUsedCh[1]++;
if (fVector2d) {
UpdateVectorCH(fXbins[0]+fd+fNfragZ[0],fAmpTotal[fd+fNfragZ[0]]/fRelativeScale);
}
}
}
}
}
}
} // Case 2 zones
}
}
//____________Offine tracking in the AliTRDtracker_____________________________
void AliTRDCalibra::ResetfVariables()
{
//
// Reset values of fAmpTotal, fPHValue and fPHPlace for
// the updateHistogram... functions
//
// Reset the good track
fGoodTrack = kTRUE;
// Reset the fAmpTotal where we put value
if (fCH2dOn) {
for (Int_t k = 0; k < fNfragZ[0]*fNfragRphi[0]; k++) {
fAmpTotal[k] = 0.0;
}
}
// Reset the fPHValue
if (fPH2dOn) {
for (Int_t k = 0; k < fTimeMax; k++) {
fPHValue[k] = -1.0;
fPHPlace[k] = -1;
}
}
}
//____________Offine tracking in the AliTRDtracker_____________________________
void AliTRDCalibra::FillTheInfoOfTheTrackPH()
{
//
// For the offline tracking or mcm tracklets
// This function will be called in the functions UpdateHistogram...
// to fill the info of a track for the drift velocity calibration
//
Int_t nb = 1; // Nombre de zones traversees 1, 2 ou plus de 3
Int_t fd1 = -1; // Premiere zone non nulle
Int_t fd2 = -1; // Deuxieme zone non nulle
Int_t k1 = -1; // Debut de la premiere zone
Int_t k2 = -1; // Debut de la seconde zone
// See if the track goes through different zones
for (Int_t k = 0; k < fTimeMax; k++) {
if (fPHValue[k] > 0.0) {
if (fd1 == -1) {
fd1 = fPHPlace[k];
k1 = k;
}
if (fPHPlace[k] != fd1) {
if (fd2 == -1) {
k2 = k;
fd2 = fPHPlace[k];
nb = 2;
}
if (fPHPlace[k] != fd2) {
nb = 3;
}
}
}
}
// Fill
// Case of track with only one zone
if (nb == 1) {
fNumberUsedPh[0]++;
for (Int_t i = 0; i < fTimeMax; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fDebug == 13) {
AliInfo(Form("WRITE nb %d ,place final: %d, fPHPlace[i]: %d, i: %d, fPHValue[i]: %f"
,nb,fXbins[1]+fPHPlace[i],fPHPlace[i],i,fPHValue[i]));
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
} // Case 1 zone
// Case of track with two zones
if (nb == 2) {
// Two zones voisines sinon rien!
// Case 1
if ((fd1 == fd2+1) ||
(fd2 == fd1+1)) {
// One of the two fast all the think
if (k2 > (k1+fDifference)) {
fNumberUsedPh[1]++;
for (Int_t i = k1; i < k2; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
if ((k2+fDifference) < fTimeMax) {
fNumberUsedPh[1]++;
for (Int_t i = k2; i < fTimeMax; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
}
// Two zones voisines sinon rien!
if (fNfragZ[1] > 1) {
// Case 2
if ((fd1+fNfragZ[1]) < (fNfragZ[1]*fNfragRphi[1])) {
if (fd2 == (fd1+fNfragZ[1])) {
// One of the two fast all the think
if (k2 > (k1+fDifference)) {
fNumberUsedPh[1]++;
for (Int_t i = k1; i < k2; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
if ((k2+fDifference) < fTimeMax) {
fNumberUsedPh[1]++;
for (Int_t i = k2; i < fTimeMax; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
}
}
// Two zones voisines sinon rien!
// Case 3
if ((fd1 - fNfragZ[1]) >= 0) {
if (fd2 == (fd1 - fNfragZ[1])) {
// One of the two fast all the think
if (k2 > (k1 + fDifference)) {
fNumberUsedPh[1]++;
for (Int_t i = k1; i < k2; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
if ((k2+fDifference) < fTimeMax) {
fNumberUsedPh[1]++;
for (Int_t i = k2; i < fTimeMax; i++) {
if (fPHValue[i] > 0.0) {
if (fHisto2d) {
fPH2d->Fill((fXbins[1]+fPHPlace[i])+0.5,(Float_t) i/fSf,(Float_t) fPHValue[i]);
}
if (fVector2d) {
UpdateVectorPH(fXbins[1]+fPHPlace[i],i,fPHValue[i]);
}
}
}
}
}
}
}
} // case 2 zones
}
//____________Set the pad calibration variables for the detector_______________
Bool_t AliTRDCalibra::LocalisationDetectorXbins(Int_t detector)
{
//
// For the detector calcul the first Xbins and set the number of row
// and col pads per calibration groups, the number of calibration
// groups in the detector.
//
// first Xbins of the detector
if (fCH2dOn) {
CalculXBins(detector,0);
}
if (fPH2dOn) {
CalculXBins(detector,1);
}
if (fPRF2dOn) {
CalculXBins(detector,2);
}
// fragmentation of idect
for (Int_t i = 0; i < 3; i++) {
ModePadCalibration((Int_t) GetChamber(detector),i);
ModePadFragmentation((Int_t) GetPlane(detector)
, (Int_t) GetChamber(detector)
, (Int_t) GetSector(detector),i);
}
return kTRUE;
}
//____________Plot the 2D histos filled Online_________________________________
//_____________________________________________________________________________
void AliTRDCalibra::PlotPH2d()
{
//
// Plot the 2D histo
//
TCanvas *cph2d = new TCanvas("cph2d","",50,50,600,800);
cph2d->cd();
fPH2d->Draw("LEGO");
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotCH2d()
{
//
// Plot the 2D histos
//
TCanvas *cch2d = new TCanvas("cch2d","",50,50,600,800);
cch2d->cd();
fCH2d->Draw("LEGO");
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotPRF2d()
{
//
// Plot the 2D histos
//
TCanvas *cPRF2d = new TCanvas("cPRF2d","",50,50,600,800);
cPRF2d->cd();
fPRF2d->Draw("LEGO");
}
//____________Fit______________________________________________________________
//____________Create histos if fDebug == 1 or fDebug >= 3______________________
//_____________________________________________________________________________
void AliTRDCalibra::InitArrayFitPH()
{
//
// Initialise fCoefVdrift[3] and fCoefVdriftE[2] to the right dimension
//
Int_t nbins = fDect2[1]-fDect1[1];
//Init the pointer to nbins
fCoefVdrift[0] = new Double_t[nbins];
fCoefVdrift[1] = new Double_t[nbins];
fCoefVdrift[2] = new Double_t[nbins];
fCoefVdriftE[0] = new Double_t[nbins];
fCoefVdriftE[1] = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k++){
fCoefVdriftE[0][k] = 0.0;
fCoefVdriftE[1][k] = 0.0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::InitArrayFitT0()
{
//
// Initialise fCoefT0[3] and fCoefT0E[2] to the right dimension
//
Int_t nbins = fDect2[1]-fDect1[1];
//Init the pointer to nbins
fCoefT0[0] = new Double_t[nbins];
fCoefT0[1] = new Double_t[nbins];
fCoefT0[2] = new Double_t[nbins];
fCoefT0E[0] = new Double_t[nbins];
fCoefT0E[1] = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k++){
fCoefT0E[0][k] = 0.0;
fCoefT0E[1][k] = 0.0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::InitArrayFitCH()
{
//
// Initialise fCoefCharge[4] and fCoefChargeE[3] to the right dimension
//
Int_t nbins = fDect2[0]-fDect1[0];
//Init the pointer to nbins
fCoefCharge[0] = new Double_t[nbins];
fCoefCharge[1] = new Double_t[nbins];
fCoefCharge[2] = new Double_t[nbins];
fCoefCharge[3] = new Double_t[nbins];
fCoefChargeE[0] = new Double_t[nbins];
fCoefChargeE[1] = new Double_t[nbins];
fCoefChargeE[2] = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k++){
fCoefChargeE[0][k] = 0.0;
fCoefChargeE[1][k] = 0.0;
fCoefChargeE[2][k] = 0.0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::InitArrayFitPRF()
{
//
// Initialise fCoefPRF[2] and fCoefPRFE to the right dimension
//
Int_t nbins = fDect2[2]-fDect1[2];
//Init the pointer to nbins
fCoefPRF[0] = new Double_t[nbins];
fCoefPRF[1] = new Double_t[nbins];
fCoefPRFE = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k++){
fCoefPRFE[k] = 0.0;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::CreateFitHistoPRFDB(Int_t rowMax, Int_t colMax)
{
//
// Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
//
fCoefPRFDB = new TH2F("coefPRF","",rowMax,0,rowMax,colMax,0,colMax);
fCoefPRFDB->SetStats(0);
fCoefPRFDB->SetXTitle("row Number");
fCoefPRFDB->SetYTitle("col Number");
fCoefPRFDB->SetZTitle("PRF width [pad width units]");
fCoefPRFDB->SetFillColor(6);
fCoefPRFDB->SetLineColor(6);
}
//_____________________________________________________________________________
void AliTRDCalibra::CreateFitHistoCHDB(Int_t rowMax, Int_t colMax)
{
//
// Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
//
fCoefChargeDB[0] = new TH2F("coefchargedb0","",rowMax,0,rowMax,colMax,0,colMax);
fCoefChargeDB[1] = new TH2F("coefchargedb1","",rowMax,0,rowMax,colMax,0,colMax);
fCoefChargeDB[2] = new TH2F("coefchargedb2","",rowMax,0,rowMax,colMax,0,colMax);
fCoefChargeDB[0]->SetStats(0);
fCoefChargeDB[1]->SetStats(0);
fCoefChargeDB[2]->SetStats(0);
fCoefChargeDB[0]->SetXTitle("row Number");
fCoefChargeDB[0]->SetYTitle("col Number");
fCoefChargeDB[1]->SetXTitle("row Number");
fCoefChargeDB[1]->SetYTitle("col Number");
fCoefChargeDB[2]->SetXTitle("row Number");
fCoefChargeDB[2]->SetYTitle("col Number");
fCoefChargeDB[0]->SetZTitle("f_{g} Fit method");
fCoefChargeDB[1]->SetZTitle("f_{g} Mean method");
fCoefChargeDB[2]->SetZTitle("f_{g} Fitbis method");
fCoefChargeDB[0]->SetFillColor(6);
fCoefChargeDB[0]->SetLineColor(6);
fCoefChargeDB[0]->SetLineColor(6);
fCoefChargeDB[1]->SetFillColor(2);
fCoefChargeDB[1]->SetLineColor(2);
fCoefChargeDB[1]->SetLineColor(2);
fCoefChargeDB[2]->SetFillColor(8);
fCoefChargeDB[2]->SetLineColor(8);
fCoefChargeDB[2]->SetLineColor(8);
}
//_____________________________________________________________________________
void AliTRDCalibra::CreateFitHistoPHDB(Int_t rowMax, Int_t colMax)
{
//
// Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
//
fCoefVdriftDB[0] = new TH2F("coefvdriftdb0","",rowMax,0,rowMax,colMax,0,colMax);
fCoefVdriftDB[1] = new TH2F("coefvdriftdb1","",rowMax,0,rowMax,colMax,0,colMax);
fCoefVdriftDB[0]->SetStats(0);
fCoefVdriftDB[1]->SetStats(0);
fCoefVdriftDB[0]->SetXTitle("row Number");
fCoefVdriftDB[0]->SetYTitle("col Number");
fCoefVdriftDB[1]->SetXTitle("row Number");
fCoefVdriftDB[1]->SetYTitle("col Number");
fCoefVdriftDB[0]->SetZTitle("v_{drift} Fit method");
fCoefVdriftDB[1]->SetZTitle("v_{drift} slope method");
fCoefVdriftDB[0]->SetFillColor(6);
fCoefVdriftDB[0]->SetLineColor(6);
fCoefVdriftDB[0]->SetLineColor(6);
fCoefVdriftDB[1]->SetFillColor(2);
fCoefVdriftDB[1]->SetLineColor(2);
fCoefVdriftDB[1]->SetLineColor(2);
}
//_____________________________________________________________________________
void AliTRDCalibra::CreateFitHistoT0DB(Int_t rowMax, Int_t colMax)
{
//
// Create the histos for fDebug = 3 and fDebug = 4 (Fit functions)
//
fCoefT0DB[0] = new TH2F("coefT0db0","",rowMax,0,rowMax,colMax,0,colMax);
fCoefT0DB[1] = new TH2F("coefT0db1","",rowMax,0,rowMax,colMax,0,colMax);
fCoefT0DB[0]->SetStats(0);
fCoefT0DB[1]->SetStats(0);
fCoefT0DB[0]->SetXTitle("row Number");
fCoefT0DB[0]->SetYTitle("col Number");
fCoefT0DB[1]->SetXTitle("row Number");
fCoefT0DB[1]->SetYTitle("col Number");
fCoefT0DB[0]->SetZTitle("t0 Fit method");
fCoefT0DB[1]->SetZTitle("t0 slope method");
fCoefT0DB[0]->SetFillColor(6);
fCoefT0DB[0]->SetLineColor(6);
fCoefT0DB[0]->SetLineColor(6);
fCoefT0DB[1]->SetFillColor(2);
fCoefT0DB[1]->SetLineColor(2);
fCoefT0DB[1]->SetLineColor(2);
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::FillVectorFitCH(Int_t countdet)
{
//
// For the Fit functions fill the vector FitCH special for the gain calibration
//
AliTRDFitCHInfo *fitCHInfo = new AliTRDFitCHInfo();
Int_t ntotal = 1;
if (GetChamber(countdet) == 2) {
ntotal = 1728;
}
else {
ntotal = 2304;
}
Float_t *coef = new Float_t[ntotal];
for (Int_t i = 0; i < ntotal; i++) {
coef[i] = fCoefCH[i];
}
Int_t detector = countdet;
// Set
fitCHInfo->SetCoef(coef);
fitCHInfo->SetDetector(detector);
fVectorFitCH->Add((TObject *) fitCHInfo);
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibra in the code_________
Bool_t AliTRDCalibra::InitFit(Int_t nbins, Int_t i)
{
//
// Init the calibration mode (Nz, Nrphi), the histograms for
// debugging the fit methods if fDebug > 0,
//
gStyle->SetPalette(1);
gStyle->SetOptStat(1111);
gStyle->SetPadBorderMode(0);
gStyle->SetCanvasColor(10);
gStyle->SetPadLeftMargin(0.13);
gStyle->SetPadRightMargin(0.01);
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam");
return kFALSE;
}
// Mode groups of pads: the total number of bins!
Int_t numberofbinsexpected = 0;
ModePadCalibration(2,i);
ModePadFragmentation(0,2,0,i);
fDetChamb2[i] = fNfragZ[i] * fNfragRphi[i];
if (fDebug == 1) {
AliInfo(Form("For the chamber 2: %d",fDetChamb2[i]));
}
numberofbinsexpected += 6 * 18 * fDetChamb2[i];
ModePadCalibration(0,i);
ModePadFragmentation(0,0,0,i);
fDetChamb0[i] = fNfragZ[i] * fNfragRphi[i];
if (fDebug == 1) {
AliInfo(Form("For the other chamber 0: %d",fDetChamb0[i]));
}
numberofbinsexpected += 6 * 4 * 18 * fDetChamb0[i];
// Quick verification that we have the good pad calibration mode if 2D histos!
if (nbins != 0) {
if (numberofbinsexpected != nbins) {
AliInfo("It doesn't correspond to the mode of pad group calibration!");
return kFALSE;
}
}
// Security for fDebug 3 and 4
if ((fDebug >= 3) &&
((fDet[0] > 5) ||
(fDet[1] > 4) ||
(fDet[2] > 17))) {
AliInfo("This detector doesn't exit!");
return kFALSE;
}
// Determine fDet1 and fDet2
fDect1[i] = -1;
fDect2[i] = -1;
if (fDebug == 2) {
fDect1[i] = fFitVoir;
fDect2[i] = fDect1[i] +1;
}
if (fDebug <= 1) {
fDect1[i] = 0;
fDect2[i] = numberofbinsexpected;
}
if (fDebug >= 3) {
CalculXBins(AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]),i);
fDect1[i] = fXbins[i];
CalculXBins((AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2])+1),i);
fDect2[i] = fXbins[i];
}
// Create the histos for debugging
// CH
if (i == 0) {
gDirectory = gROOT;
// Init the VectorFitCH
fVectorFitCH = new TObjArray();
fCoefCH = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCoefCH[k] = 0.0;
}
fScaleFitFactor = 0.0;
// Number of Xbins(detectors or groups of pads) if Vector2d
// Quick verification that we are not out of range!
if (fVectorCH && fPlaCH) {
if ((nbins == 0) &&
(fVectorCH->GetEntriesFast() > 0) &&
((Int_t) fPlaCH->GetEntriesFast() > 0)) {
if ((Int_t) fVectorCH->GetEntriesFast() > numberofbinsexpected) {
AliInfo("ch doesn't correspond to the mode of pad group calibration!");
return kFALSE;
}
if ((Int_t) fVectorCH->GetEntriesFast() != (Int_t) fPlaCH->GetEntriesFast()) {
AliInfo("VectorCH doesn't correspond to PlaCH!");
return kFALSE;
}
}
}
//
// Debugging: Create the histos
//
// fDebug == 0 nothing
// fDebug == 1
if (fDebug == 1) {
InitArrayFitCH();
}
// fDebug == 2 and fFitVoir no histo
if (fDebug == 2) {
if (fFitVoir < numberofbinsexpected) {
AliInfo(Form("We will see the fit of the object %d",fFitVoir));
}
else {
AliInfo("fFitVoir is out of range of the histo!");
return kFALSE;
}
}
// fDebug == 3 or 4 and fDet
if (fDebug >= 3) {
if ((fNz[0] == 0) && (fNrphi[0] == 0)) {
AliInfo("Do you really want to see one detector without pad groups?");
return kFALSE;
}
else {
AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
,fDet[0],fDet[1],fDet[2]));
// A little geometry:
Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
Int_t colMax = parCom->GetColMax(fDet[0]);
// Create the histos to visualise
CreateFitHistoCHDB(rowMax,colMax);
if (fDebug == 4) {
InitArrayFitCH();
}
}
}
}
// PH and T0
if (i == 1) {
// Number of Xbins (detectors or groups of pads) if vector2d
// Quick verification that we are not out of range!
if (fVectorPH && fPlaPH) {
if ((nbins == 0) &&
(fVectorPH->GetEntriesFast() > 0) &&
((Int_t) fPlaPH->GetEntriesFast() > 0)) {
if ((Int_t) fVectorPH->GetEntriesFast() > numberofbinsexpected) {
AliInfo("ph doesn't correspond to the mode of pad group calibration!");
return kFALSE;
}
if ((Int_t) fVectorPH->GetEntriesFast() != (Int_t) fPlaPH->GetEntriesFast()) {
AliInfo("VectorPH doesn't correspond to PlaPH!");
return kFALSE;
}
}
}
// Init tree
InitTreePH();
InitTreeT0();
//
// Debugging: Create the histos
//
// fDebug == 0 nothing
// fDebug == 1
if (fDebug == 1) {
// Create the histos replique de ph
InitArrayFitPH();
InitArrayFitT0();
}
// fDebug == 2 and fFitVoir no histo
if (fDebug == 2) {
if (fFitVoir < numberofbinsexpected) {
AliInfo(Form("We will see the fit of the object %d",fFitVoir));
}
else {
AliInfo("fFitVoir is out of range of the histo!");
return kFALSE;
}
}
// fDebug == 3 or 4 and fDet
if (fDebug >= 3) {
if ((fNz[1] == 0) &&
(fNrphi[1] == 0)) {
AliInfo("Do you really want to see one detector without pad groups?");
return kFALSE;
}
else {
AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
,fDet[0],fDet[1],fDet[2]));
// A little geometry:
Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
Int_t colMax = parCom->GetColMax(fDet[0]);
// Create the histos to visualise
CreateFitHistoPHDB(rowMax,colMax);
CreateFitHistoT0DB(rowMax,colMax);
if (fDebug == 4) {
InitArrayFitPH();
InitArrayFitT0();
}
}
}
}
// PRF
if (i == 2) {
// Number of Xbins(detectors or groups of pads) if vector2d
if (fVectorPRF && fPlaPRF){
if ((nbins == 0) &&
(fVectorPRF->GetEntriesFast() > 0) &&
(fPlaPRF->GetEntriesFast() > 0)) {
// Quick verification that we are not out of range!
if ((Int_t) fVectorPRF->GetEntriesFast() > numberofbinsexpected) {
AliInfo("ch doesn't correspond to the mode of pad group calibration!");
return kFALSE;
}
if ((Int_t) fVectorPRF->GetEntriesFast() != (Int_t) fPlaPRF->GetEntriesFast()) {
AliInfo("VectorPRF doesn't correspond to PlaCH!");
return kFALSE;
}
}
}
// Init tree
InitTreePRF();
//
// Debugging: Create the histos
//
// fDebug == 0 nothing
// fDebug == 1
if (fDebug == 1) {
// Create the histos replique de ch
InitArrayFitPRF();
}
// fDebug == 2 and fFitVoir no histo
if (fDebug == 2) {
if (fFitVoir < numberofbinsexpected) {
AliInfo(Form("We will see the fit of the object %d",fFitVoir));
}
else {
AliInfo("fFitVoir is out of range of the histo!");
return kFALSE;
}
}
// fDebug == 3 or 4 and fDet
if (fDebug >= 3) {
if ((fNz[2] == 0) &&
(fNrphi[2] == 0)) {
AliInfo("Do you really want to see one detector without pad groups?");
return kFALSE;
}
else {
AliInfo(Form("You will see the detector: iPlane %d, iChamb %d, iSect %d"
,fDet[0],fDet[1],fDet[2]));
// A little geometry:
Int_t rowMax = parCom->GetRowMax(fDet[0],fDet[1],fDet[2]);
Int_t colMax = parCom->GetColMax(fDet[0]);
// Create the histos to visualise
CreateFitHistoPRFDB(rowMax,colMax);
if (fDebug == 4) {
InitArrayFitPRF();
}
}
}
}
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibra in the code_________
void AliTRDCalibra::InitfCountDetAndfCount(Int_t i)
{
//
// Init the current detector where we are fCountDet and the
// next fCount for the functions Fit...
//
// Loop on the Xbins of ch!!
fCountDet[i] = -1; // Current detector
fCount[i] = 0; // To find the next detector
// If fDebug >= 3
if (fDebug >= 3) {
// Set countdet to the detector
fCountDet[i] = AliTRDgeometry::GetDetector(fDet[0],fDet[1],fDet[2]);
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
ModePadCalibration(fDet[1],i);
ModePadFragmentation(fDet[0],fDet[1],fDet[2],i);
// Set counter to write at the end of the detector
fCount[i] = fDect1[i] + fNfragZ[i]*fNfragRphi[i];
}
}
//____________Functions for initialising the AliTRDCalibra in the code_________
void AliTRDCalibra::UpdatefCountDetAndfCount(Int_t idect, Int_t i)
{
//
// See if we are in a new detector and update the
// variables fNfragZ and fNfragRphi if yes
//
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
// If fDebug == 1 or 0
if ((fDebug == 0) ||
(fDebug == 1)) {
if (fCount[i] == idect) {
// On en est au detector
fCountDet[i] += 1;
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
ModePadCalibration((Int_t) GetChamber(fCountDet[i]),i);
ModePadFragmentation((Int_t) GetPlane(fCountDet[i])
,(Int_t) GetChamber(fCountDet[i])
,(Int_t) GetSector(fCountDet[i]),i);
// Set for the next detector
fCount[i] += fNfragZ[i]*fNfragRphi[i];
}
}
}
//____________Functions for initialising the AliTRDCalibra in the code_________
void AliTRDCalibra::ReconstructFitRowMinRowMax(Int_t idect, Int_t i)
{
//
// Reconstruct the min pad row, max pad row, min pad col and
// max pad col of the calibration group for the Fit functions
//
if (fDebug < 2) {
ReconstructionRowPadGroup((Int_t) (idect-(fCount[i]-(fNfragZ[i]*fNfragRphi[i]))),i);
}
if (fDebug >= 3) {
ReconstructionRowPadGroup((Int_t) (idect-fDect1[i]),i);
}
}
//____________Functions for initialising the AliTRDCalibra in the code_________
Bool_t AliTRDCalibra::NotEnoughStatistic(Int_t idect, Int_t i)
{
//
// For the case where there are not enough entries in the histograms
// of the calibration group, the value present in the choosen database
// will be put. A negativ sign enables to know that a fit was not possible.
//
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam Manager");
return kFALSE;
}
// Get cal
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
if (fDebug != 2) {
AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
,idect-(fCount[i]-(fNfragZ[i]*fNfragRphi[i])),fCountDet[i]));
}
if (fDebug == 2) {
AliInfo("The element has not enough statistic to be fitted");
}
if ((i == 0) && (fDebug != 2)) {
// Calcul the coef from the database choosen
CalculChargeCoefMean(fCountDet[0],(Int_t) (idect-fDect1[0]),kFALSE);
// Fill the coefCH[2304] with negative value to say: not fitted
for (Int_t k = fRowMin[0]; k < fRowMax[0]; k++) {
for (Int_t j = fColMin[0]; j < fColMax[0]; j++) {
if (GetChamber(fCountDet[0]) == 2) {
fCoefCH[(Int_t)(j*12+k)] = -TMath::Abs(fChargeCoef[3]);
}
if (GetChamber(fCountDet[0]) != 2) {
fCoefCH[(Int_t)(j*16+k)] = -TMath::Abs(fChargeCoef[3]);
}
}
}
// Put the default value negative
if ((fDebug == 1) ||
(fDebug == 4)) {
if (fFitChargeBisOn) {
fCoefCharge[2][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
}
if (fMeanChargeOn) {
fCoefCharge[1][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
}
fCoefCharge[0][idect-fDect1[0]]=-TMath::Abs(fChargeCoef[3]);
}
// End of one detector
if ((idect == (fCount[0]-1))) {
FillVectorFitCH((Int_t) fCountDet[0]);
// Reset
for (Int_t k = 0; k < 2304; k++) {
fCoefCH[k] = 0.0;
}
}
}
if ((i == 1) && (fDebug != 2)) {
CalculVdriftCoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
CalculT0CoefMean(fCountDet[1],(Int_t) (idect-fDect1[1]));
// Put the default value (time0 can be negativ, so we stay with + )
if ((fDebug == 1) ||
(fDebug == 4)) {
if (fFitPHOn) {
fCoefVdrift[0][(idect-fDect1[1])] = -fVdriftCoef[2];
fCoefT0[0][(idect-fDect1[1])] = fT0Coef[2];
}
fCoefVdrift[1][(idect-fDect1[1])] = -fVdriftCoef[2];
fCoefT0[1][(idect-fDect1[1])] = fT0Coef[2];
}
// Put the default value
if (fDebug >= 3) {
fVdriftCoef[0] = fVdriftCoef[2];
fVdriftCoef[1] = fVdriftCoef[2];
FillCoefVdriftDB();
fT0Coef[0] = fT0Coef[2];
fT0Coef[1] = fT0Coef[2];
FillCoefT0DB();
}
// Fill the tree if end of a detector.
// The pointer to the branch stays with the default value negative!!!
// PH
// Pointer to the branch
for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
if (GetChamber(fCountDet[1]) == 2) {
fVdriftPad[(Int_t)(j*12+k)] = -TMath::Abs(fVdriftCoef[2]);
}
if (GetChamber(fCountDet[1]) != 2) {
fVdriftPad[(Int_t)(j*16+k)] = -TMath::Abs(fVdriftCoef[2]);
}
}
}
// End of one detector
if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
FillTreeVdrift((Int_t) fCountDet[1]);
}
// T0
// Fill the tree if end of a detector.
// The pointer to the branch stays with the default value positive!!!
// Pointer to the branch
for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
if (GetChamber(fCountDet[1]) == 2) {
fT0Pad[(Int_t)(j*12+k)] = fT0Coef[2];
}
if (GetChamber(fCountDet[1]) != 2) {
fT0Pad[(Int_t)(j*16+k)] = fT0Coef[2];
}
}
}
// End of one detector
if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
FillTreeT0((Int_t) fCountDet[1]);
}
}
if ((i == 2) && (fDebug != 2)) {
CalculPRFCoefMean(fCountDet[2],(Int_t) (idect-fDect1[2]));
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefPRF[0][(idect-fDect1[2])] = -fPRFCoef[1];
}
if (fDebug >= 3){
fPRFCoef[0] = fPRFCoef[1];
FillCoefPRFDB();
}
// Fill the tree if end of a detector.
// The pointer to the branch stays with the default value 1.5!!!
// Pointer to the branch
for (Int_t k = fRowMin[2]; k < fRowMax[2]; k++) {
for (Int_t j = fColMin[2]; j < fColMax[2]; j++) {
if((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)){
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = -fPRFCoef[1];
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = -fPRFCoef[1];
}
}
else {
if (fAccCDB) {
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = -((Float_t) cal->GetPRFWidth(fCountDet[2],j,k));
}
}
if (!fAccCDB) {
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = -((Float_t) GetPRFDefault(GetPlane(fCountDet[2])));
}
}
}
}
}
// End of one detector
if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
FillTreePRF((Int_t) fCountDet[2]);
}
}
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibra in the code_________
Bool_t AliTRDCalibra::FillInfosFit(Int_t idect, Int_t i)
{
//
// Fill the coefficients found with the fits or other
// methods from the Fit functions
//
// Get the parameter object
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam Manager");
return kFALSE;
}
// Get cal
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB");
return kFALSE;
}
if ((i == 0) && (fDebug != 2)) {
// Fill the coefCH[2304] with fChargeCoef[0]
// that would be negativ only if the fit failed totally
for (Int_t k = fRowMin[0]; k < fRowMax[0]; k++) {
for (Int_t j = fColMin[0]; j < fColMax[0]; j++) {
if (GetChamber(fCountDet[0]) == 2) {
fCoefCH[(Int_t)(j*12+k)] = fChargeCoef[0];
}
if (GetChamber(fCountDet[0]) != 2) {
fCoefCH[(Int_t)(j*16+k)] = fChargeCoef[0];
}
}
}
// End of one detector
if ((idect == (fCount[0]-1))) {
FillVectorFitCH((Int_t) fCountDet[0]);
// Reset
for (Int_t k = 0; k < 2304; k++) {
fCoefCH[k] = 0.0;
}
}
}
if ((i == 1) && (fDebug != 2)) {
// PH
// Pointer to the branch: fVdriftCoef[1] will ne negativ only if the fit failed totally
for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
if (GetChamber(fCountDet[1]) == 2) {
fVdriftPad[(Int_t)(j*12+k)]=fVdriftCoef[1];
}
if (GetChamber(fCountDet[1]) != 2) {
fVdriftPad[(Int_t)(j*16+k)]=fVdriftCoef[1];
}
}
}
// End of one detector
if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
FillTreeVdrift((Int_t) fCountDet[1]);
}
// T0
// Pointer to the branch: fT0Coef[1] will ne negativ only if the fit failed totally
for (Int_t k = fRowMin[1]; k < fRowMax[1]; k++) {
for (Int_t j = fColMin[1]; j < fColMax[1]; j++) {
if (GetChamber(fCountDet[1]) == 2) {
fT0Pad[(Int_t)(j*12+k)]=fT0Coef[1];
}
if (GetChamber(fCountDet[1]) != 2) {
fT0Pad[(Int_t)(j*16+k)]=fT0Coef[1];
}
}
}
// End of one detector
if ((idect == (fCount[1]-1)) && (fDebug != 2)) {
FillTreeT0((Int_t) fCountDet[1]);
}
}
if ((i == 2) && (fDebug != 2)) {
// Pointer to the branch
for (Int_t k = fRowMin[2]; k < fRowMax[2]; k++) {
for (Int_t j = fColMin[2]; j < fColMax[2]; j++) {
if ((parCom->GetColMax(GetPlane(fCountDet[2])) != (j+1)) && (j != 0)) {
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = fPRFCoef[0];
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = fPRFCoef[0];
}
}
else {
if (fAccCDB) {
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = (Float_t) cal->GetPRFWidth(fCountDet[2],j,k);
}
}
if (!fAccCDB) {
if (GetChamber(fCountDet[2]) == 2) {
fPRFPad[(Int_t)(j*12+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
}
if (GetChamber(fCountDet[2]) != 2) {
fPRFPad[(Int_t)(j*16+k)] = (Float_t) GetPRFDefault(GetPlane(fCountDet[2]));
}
}
}
}
}
// End of one detector
if ((idect == (fCount[2]-1)) && (fDebug != 2)) {
FillTreePRF((Int_t) fCountDet[2]);
}
}
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibra in the code_________
Bool_t AliTRDCalibra::WriteFitInfos(Int_t i)
{
//
// In the case the user wants to write a file with a tree of the found
// coefficients for the calibration before putting them in the database
//
TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
return kFALSE;
}
if ((i == 0) && (fDebug != 2)) {
// The DB stuff
if ((fDebug == 4) ||
(fDebug == 3)) {
WriteCHDB(fout);
}
// The tree
fout->WriteTObject(fGain,fGain->GetName(),(Option_t *) "writedelete");
}
if ((i == 1) && (fDebug != 2)) {
// The DB stuff
if ((fDebug == 4) ||
(fDebug == 3)) {
WritePHDB(fout);
}
// The tree
fout->WriteTObject(fVdrift,fVdrift->GetName(),(Option_t *) "writedelete");
// The DB stuff
if ((fDebug == 4) ||
(fDebug == 3)) {
WriteT0DB(fout);
}
// The tree
fout->WriteTObject(fT0,fT0->GetName(),(Option_t *) "writedelete");
}
if ((i == 2) && (fDebug != 2)) {
// The DB stuff
if ((fDebug == 4) ||
(fDebug == 3)) {
WritePRFDB(fout);
}
// The tree
fout->WriteTObject(fPRF,fPRF->GetName(),(Option_t *) "writedelete");
}
fout->Close();
return kTRUE;
}
//
//____________Fill Coef DB in case of visualisation of one detector____________
//
//_____________________________________________________________________________
void AliTRDCalibra::FillCoefVdriftDB()
{
//
// Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
fCoefVdriftDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[1]));
if (fFitPHOn ) {
fCoefVdriftDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fVdriftCoef[0]));
}
}
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FillCoefT0DB()
{
//
// Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
fCoefT0DB[1]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[1]));
if (fFitPHOn) {
fCoefT0DB[0]->SetBinContent(row+1,col+1,TMath::Abs(fT0Coef[0]));
}
}
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FillCoefChargeDB()
{
//
// Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
for (Int_t row = fRowMin[0]; row < fRowMax[0]; row++) {
for (Int_t col = fColMin[0]; col < fColMax[0]; col++) {
if (fMeanChargeOn) {
fCoefChargeDB[1]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[1]));
}
if (fFitChargeBisOn) {
fCoefChargeDB[2]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[2]));
}
fCoefChargeDB[0]->SetBinContent(row+1,col+1,TMath::Abs(fChargeCoef[0]));
}
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FillCoefPRFDB()
{
//
// Fill the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
for (Int_t row = fRowMin[2]; row < fRowMax[2]; row++) {
for (Int_t col = fColMin[2]; col < fColMax[2]; col++) {
fCoefPRFDB->SetBinContent(row+1,col+1,fPRFCoef[0]);
}
}
}
//
//____________Plot histos CoefPRF....__________________________________________
//
//_____________________________________________________________________________
void AliTRDCalibra::PlotWriteCH()
{
//
// Scale the coefficients to one, create the graph errors and write them if wanted
//
//TObjArray of the grapherrors and so on
TObjArray *listofgraphs = new TObjArray();
Int_t nbins = fDect2[0]-fDect1[0];
//Scale the coefs
//counter
Int_t counter[3];
counter[0] = 0;
counter[1] = 0;
counter[2] = 0;
Double_t sum = 0.0;
Double_t scale = 1.0;
// Scale the histo
Double_t *xValuesFitted = new Double_t[nbins];
Double_t *xValuesFittedMean = new Double_t[nbins];
Double_t *xValuesFittedBis = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValuesFitted[k] = -1;
xValuesFittedMean[k] = -1;
xValuesFittedBis[k] = -1;
}
for(Int_t l = 0; l < nbins; l++){
if(fCoefCharge[0][l] > 0){
fCoefCharge[0][l]=fCoefCharge[0][l]*fScaleFitFactor;
fCoefChargeE[0][l]=fCoefChargeE[0][l]*fScaleFitFactor;
xValuesFitted[counter[0]]=l;
counter[0]++;
}
}
if(fMeanChargeOn){
sum = 0.0;
for(Int_t l = 0; l < nbins; l++){
if(fCoefCharge[1][l] > 0){
sum += fCoefCharge[1][l];
xValuesFittedMean[counter[1]]= l;
counter[1]++;
}
}
scale = 1.0;
if(sum > 0.0) scale = counter[1]/sum;
for(Int_t l = 0; l < nbins; l++){
if(fCoefCharge[1][l] > 0){
fCoefCharge[1][l]=fCoefCharge[1][l]*scale;
fCoefChargeE[1][l]=fCoefChargeE[1][l]*scale;
}
}
}
if(fFitChargeBisOn){
sum = 0.0;
for(Int_t l = 0; l < nbins; l++){
if(fCoefCharge[2][l] > 0){
fCoefCharge[2][l]=fCoefCharge[2][l]*fScaleFitFactor;
fCoefChargeE[2][l]=fCoefChargeE[2][l]*fScaleFitFactor;
sum += fCoefCharge[2][l];
xValuesFittedBis[counter[2]]= l;
counter[2]++;
}
}
scale = 1.0;
if(sum > 0.0) scale = counter[2]/sum;
for(Int_t l = 0; l < nbins; l++){
if(fCoefCharge[2][l] > 0){
fCoefCharge[2][l]=fCoefCharge[2][l]*scale;
fCoefChargeE[2][l]=fCoefChargeE[2][l]*scale;
if(fCoefCharge[2][l] > 1.0) printf("for the group %d, I have the coef %f with the error %f\n",l,fCoefCharge[2][l],fCoefChargeE[2][l]);
}
}
}
//Create the X and Xerror
Double_t *xValues = new Double_t[nbins];
Double_t *xValuesE = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValues[k] = k;
xValuesE[k] = 0.0;
}
//Create the graph erros and plot them
TGraph *graphCharge3 = new TGraph(nbins,xValues,fCoefCharge[3]);
graphCharge3->SetName("coefcharge3");
graphCharge3->SetTitle("");
graphCharge3->GetXaxis()->SetTitle("Det/Pad groups");
graphCharge3->GetYaxis()->SetTitle("gain factor");
graphCharge3->SetLineColor(4);
listofgraphs->Add((TObject *)graphCharge3);
TGraphErrors *graphCharge0 = new TGraphErrors(nbins,xValues,fCoefCharge[0],xValuesE,fCoefChargeE[0]);
graphCharge0->SetName("coefcharge0");
graphCharge0->SetTitle("");
graphCharge0->GetXaxis()->SetTitle("Det/Pad groups");
graphCharge0->GetYaxis()->SetTitle("gain factor");
graphCharge0->SetMarkerColor(6);
graphCharge0->SetLineColor(6);
graphCharge0->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphCharge0);
TCanvas *cch1 = new TCanvas("cch1","",50,50,600,800);
cch1->cd();
TLegend *legch1 = new TLegend(0.4,0.6,0.89,0.89);
legch1->AddEntry(graphCharge3,"f_{g} simulated","l");
legch1->AddEntry(graphCharge0,"f_{g} fit","p");
graphCharge0->Draw("AP");
//graphCharge3->Draw("AL");
if (fMeanChargeOn) {
TGraphErrors *graphCharge1 = new TGraphErrors(nbins,xValues,fCoefCharge[1],xValuesE,fCoefChargeE[1]);
graphCharge1->SetName("coefcharge1");
graphCharge1->GetXaxis()->SetTitle("Det/Pad groups");
graphCharge1->GetYaxis()->SetTitle("gain factor");
graphCharge1->SetTitle("");
graphCharge1->SetMarkerColor(2);
graphCharge1->SetLineColor(2);
graphCharge1->SetMarkerStyle(24);
legch1->AddEntry(graphCharge1,"f_{g} mean","p");
graphCharge1->Draw("P");
listofgraphs->Add((TObject *)graphCharge1);
}
if (fFitChargeBisOn ) {
TGraphErrors *graphCharge2 = new TGraphErrors(nbins,xValues,fCoefCharge[2],xValuesE,fCoefChargeE[2]);
graphCharge2->SetName("coefcharge2");
graphCharge2->SetTitle("");
graphCharge2->GetXaxis()->SetTitle("Det/Pad groups");
graphCharge2->GetYaxis()->SetTitle("gain factor");
graphCharge2->SetMarkerColor(8);
graphCharge2->SetLineColor(8);
graphCharge2->SetMarkerStyle(25);
legch1->AddEntry(graphCharge2,"f_{g} fitbis","p");
graphCharge2->Draw("P");
listofgraphs->Add((TObject *)graphCharge2);
}
legch1->Draw("same");
//Create the arrays and the graphs for the delta
TCanvas *cch2 = new TCanvas("cch2","",50,50,600,800);
cch2->Divide(2,1);
cch2->cd(2);
Double_t *yValuesDelta = new Double_t[counter[0]];
for(Int_t k = 0; k < counter[0]; k++){
if(fCoefCharge[3][(Int_t)(xValuesFitted[k])] > 0.0) {
yValuesDelta[k] = (fCoefCharge[0][(Int_t)xValuesFitted[k]]-fCoefCharge[3][(Int_t)xValuesFitted[k]])/fCoefCharge[3][(Int_t)xValuesFitted[k]];
}
else yValuesDelta[k] = 0.0;
}
TGraph *graphDeltaCharge0 = new TGraph(counter[0],&xValuesFitted[0],yValuesDelta);
graphDeltaCharge0->SetName("deltacharge0");
graphDeltaCharge0->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaCharge0->GetYaxis()->SetTitle("#Deltag/g_{sim}");
graphDeltaCharge0->SetMarkerColor(6);
graphDeltaCharge0->SetTitle("");
graphDeltaCharge0->SetLineColor(6);
graphDeltaCharge0->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphDeltaCharge0);
TLegend *legch3 = new TLegend(0.4,0.6,0.89,0.89);
legch3->AddEntry(graphDeltaCharge0,"fit","p");
graphDeltaCharge0->Draw("AP");
cch2->cd(1);
TH1I *histoErrorCharge0 = new TH1I("errorcharge0","",100 ,-0.3,0.3);
histoErrorCharge0->SetXTitle("#Deltag/g_{sim}");
histoErrorCharge0->SetYTitle("counts");
histoErrorCharge0->SetLineColor(6);
histoErrorCharge0->SetLineStyle(1);
histoErrorCharge0->SetStats(0);
for(Int_t k = 0; k < counter[0]; k++){
histoErrorCharge0->Fill(yValuesDelta[k]);
}
TLegend *legch2 = new TLegend(0.4,0.6,0.89,0.89);
legch2->AddEntry(histoErrorCharge0,"f_{g} fit","l");
histoErrorCharge0->Draw();
listofgraphs->Add((TObject *)histoErrorCharge0);
if (fMeanChargeOn) {
cch2->cd(2);
Double_t *yValuesDeltaMean = new Double_t[counter[1]];
for(Int_t k = 0; k < counter[1]; k++){
if(fCoefCharge[3][(Int_t)xValuesFittedMean[k]] > 0.0) {
yValuesDeltaMean[k] = (fCoefCharge[1][(Int_t)xValuesFittedMean[k]]-fCoefCharge[3][(Int_t)xValuesFittedMean[k]])/fCoefCharge[3][(Int_t)xValuesFittedMean[k]];
}
else yValuesDeltaMean[k] = 0.0;
}
TGraph *graphDeltaCharge1 = new TGraph(counter[1],&xValuesFittedMean[0],yValuesDeltaMean);
graphDeltaCharge1->SetName("deltacharge1");
graphDeltaCharge1->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaCharge1->GetYaxis()->SetTitle("#Deltag/g_{sim}");
graphDeltaCharge1->SetMarkerColor(2);
graphDeltaCharge1->SetMarkerStyle(24);
graphDeltaCharge1->SetLineColor(2);
graphDeltaCharge1->SetTitle("");
legch3->AddEntry(graphDeltaCharge1,"mean","p");
graphDeltaCharge1->Draw("P");
listofgraphs->Add((TObject *)graphDeltaCharge1);
cch2->cd(1);
TH1I *histoErrorCharge1 = new TH1I("errorcharge1","",100 ,-0.3,0.3);
histoErrorCharge1->SetXTitle("#Deltag/g_{sim}");
histoErrorCharge1->SetYTitle("counts");
histoErrorCharge1->SetLineColor(2);
histoErrorCharge1->SetLineStyle(2);
histoErrorCharge1->SetStats(0);
for(Int_t k = 0; k < counter[1]; k++){
histoErrorCharge1->Fill(yValuesDeltaMean[k]);
}
legch2->AddEntry(histoErrorCharge1,"f_{g} mean","l");
histoErrorCharge1->Draw("same");
listofgraphs->Add((TObject *)histoErrorCharge1);
}
if (fFitChargeBisOn) {
cch2->cd(2);
Double_t *yValuesDeltaBis = new Double_t[counter[2]];
for(Int_t k = 0; k < counter[2]; k++){
if(fCoefCharge[3][(Int_t)xValuesFittedBis[k]] > 0.0) {
yValuesDeltaBis[k] = (fCoefCharge[2][(Int_t)xValuesFittedBis[k]]-fCoefCharge[3][(Int_t)xValuesFittedBis[k]])/fCoefCharge[3][(Int_t)xValuesFittedBis[k]];
}
else yValuesDeltaBis[k] = 0.0;
}
TGraph *graphDeltaCharge2 = new TGraph(counter[2],&xValuesFittedBis[0],yValuesDeltaBis);
graphDeltaCharge2->SetName("deltacharge2");
graphDeltaCharge2->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaCharge2->GetYaxis()->SetTitle("#Deltag/g_{sim}");
graphDeltaCharge2->SetMarkerColor(8);
graphDeltaCharge2->SetLineColor(8);
graphDeltaCharge2->SetMarkerStyle(25);
legch3->AddEntry(graphDeltaCharge2,"fit","p");
graphDeltaCharge2->SetTitle("");
graphDeltaCharge2->Draw("P");
listofgraphs->Add((TObject *)graphDeltaCharge2);
cch2->cd(1);
TH1I *histoErrorCharge2 = new TH1I("errorcharge2","",100 ,-0.3,0.3);
histoErrorCharge2->SetXTitle("#Deltag/g_{sim}");
histoErrorCharge2->SetYTitle("counts");
histoErrorCharge2->SetLineColor(8);
histoErrorCharge2->SetLineStyle(5);
histoErrorCharge2->SetLineWidth(3);
histoErrorCharge2->SetStats(0);
for(Int_t k = 0; k < counter[2]; k++){
histoErrorCharge2->Fill(yValuesDeltaBis[k]);
}
legch2->AddEntry(histoErrorCharge2,"f_{g} fitbis","l");
histoErrorCharge2->Draw("same");
listofgraphs->Add((TObject *)histoErrorCharge2);
}
cch2->cd(2);
legch3->Draw("same");
cch2->cd(1);
legch2->Draw("same");
//Write if wanted
if (fWriteCoef[0]){
TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
}
else{
for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
}
}
fout->Close();
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotWritePH()
{
//
// create the graph errors and write them if wanted
//
//TObjArray of the grapherrors and so on
TObjArray *listofgraphs = new TObjArray();
Int_t nbins = fDect2[1]-fDect1[1];
//See the number of fitted for delta
//counter
Int_t counter[2];
counter[0] = 0;
counter[1] = 0;
Double_t *xValuesFitted = new Double_t[nbins];
Double_t *xValuesFittedPH = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValuesFitted[k] = -1;
xValuesFittedPH[k] = -1;
}
for(Int_t l = 0; l < nbins; l++){
if(fCoefVdrift[1][l] > 0){
xValuesFitted[counter[1]]=l;
counter[1]++;
}
}
if(fFitPHOn){
for(Int_t l = 0; l < nbins; l++){
if(fCoefVdrift[0][l] > 0){
xValuesFittedPH[counter[0]]= l;
counter[0]++;
}
}
}
//Create the X and Xerror
Double_t *xValues = new Double_t[nbins];
Double_t *xValuesE = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValues[k] = k;
xValuesE[k] = 0.0;
}
//Create the graph erros and plot them
TGraph *graphVdrift2 = new TGraph(nbins,xValues,fCoefVdrift[2]);
graphVdrift2->SetName("coefvdrift2");
graphVdrift2->SetTitle("");
graphVdrift2->GetXaxis()->SetTitle("Det/Pad groups");
graphVdrift2->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
graphVdrift2->SetLineColor(4);
listofgraphs->Add((TObject *)graphVdrift2);
TGraphErrors *graphVdrift1 = new TGraphErrors(nbins,xValues,fCoefVdrift[1],xValuesE,fCoefVdriftE[1]);
graphVdrift1->SetName("coefvdrift1");
graphVdrift1->SetTitle("");
graphVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
graphVdrift1->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
graphVdrift1->SetMarkerColor(6);
graphVdrift1->SetLineColor(6);
graphVdrift1->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphVdrift1);
TCanvas *cph1 = new TCanvas("cph1","",50,50,600,800);
cph1->cd();
TLegend *legph1 = new TLegend(0.4,0.6,0.89,0.89);
legph1->AddEntry(graphVdrift2,"Vdrift simulated","l");
legph1->AddEntry(graphVdrift1,"Vdrift fit","p");
graphVdrift2->Draw("AL");
graphVdrift1->Draw("P");
if (fFitPHOn) {
TGraphErrors *graphVdrift0 = new TGraphErrors(nbins,xValues,fCoefVdrift[0],xValuesE,fCoefVdriftE[0]);
graphVdrift0->SetName("coefVdrift0");
graphVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
graphVdrift0->GetYaxis()->SetTitle("Vdrift [cm/#mus]");
graphVdrift0->SetTitle("");
graphVdrift0->SetMarkerColor(2);
graphVdrift0->SetLineColor(2);
graphVdrift0->SetMarkerStyle(24);
legph1->AddEntry(graphVdrift0,"v_{fit PH}","p");
graphVdrift0->Draw("P");
listofgraphs->Add((TObject *)graphVdrift0);
}
legph1->Draw("same");
//Create the arrays and the graphs for the delta
TCanvas *cph2 = new TCanvas("cph2","",50,50,600,800);
cph2->Divide(2,1);
cph2->cd(2);
Double_t *yValuesDelta = new Double_t[counter[1]];
for(Int_t k = 0; k < counter[1]; k++){
if(fCoefVdrift[2][(Int_t)(xValuesFitted[k])] > 0.0) {
yValuesDelta[k] = (fCoefVdrift[1][(Int_t)xValuesFitted[k]]-fCoefVdrift[2][(Int_t)xValuesFitted[k]])/fCoefVdrift[2][(Int_t)xValuesFitted[k]];
}
else yValuesDelta[k] = 0.0;
}
TGraph *graphDeltaVdrift1 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
graphDeltaVdrift1->SetName("deltavdrift1");
graphDeltaVdrift1->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaVdrift1->GetYaxis()->SetTitle("#Deltav/v_{sim}");
graphDeltaVdrift1->SetMarkerColor(6);
graphDeltaVdrift1->SetTitle("");
graphDeltaVdrift1->SetLineColor(6);
graphDeltaVdrift1->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphDeltaVdrift1);
TLegend *legph3 = new TLegend(0.4,0.6,0.89,0.89);
legph3->AddEntry(graphDeltaVdrift1,"v_{slope method}","p");
graphDeltaVdrift1->Draw("AP");
cph2->cd(1);
TH1I *histoErrorVdrift1 = new TH1I("errorvdrift1","",100 ,-0.05,0.05);
histoErrorVdrift1->SetXTitle("#Deltav/v_{sim}");
histoErrorVdrift1->SetYTitle("counts");
histoErrorVdrift1->SetLineColor(6);
histoErrorVdrift1->SetLineStyle(1);
histoErrorVdrift1->SetStats(0);
for(Int_t k = 0; k < counter[1]; k++){
histoErrorVdrift1->Fill(yValuesDelta[k]);
}
TLegend *legph2 = new TLegend(0.4,0.6,0.89,0.89);
legph2->AddEntry(histoErrorVdrift1,"v_{slope method}","l");
histoErrorVdrift1->Draw();
listofgraphs->Add((TObject *)histoErrorVdrift1);
if (fFitPHOn) {
cph2->cd(2);
Double_t *yValuesDeltaPH = new Double_t[counter[0]];
for(Int_t k = 0; k < counter[0]; k++){
if(fCoefVdrift[2][(Int_t)xValuesFittedPH[k]] > 0.0) {
yValuesDeltaPH[k] = (fCoefVdrift[0][(Int_t)xValuesFittedPH[k]]-fCoefVdrift[2][(Int_t)xValuesFittedPH[k]])/fCoefVdrift[2][(Int_t)xValuesFittedPH[k]];
}
else yValuesDeltaPH[k] = 0.0;
}
TGraph *graphDeltaVdrift0 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
graphDeltaVdrift0->SetName("deltavdrift0");
graphDeltaVdrift0->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaVdrift0->GetYaxis()->SetTitle("#Deltav/v_{sim}");
graphDeltaVdrift0->SetMarkerColor(2);
graphDeltaVdrift0->SetMarkerStyle(24);
graphDeltaVdrift0->SetLineColor(2);
graphDeltaVdrift0->SetTitle("");
legph3->AddEntry(graphDeltaVdrift0,"v_{fit PH}","p");
graphDeltaVdrift0->Draw("P");
listofgraphs->Add((TObject *)graphDeltaVdrift0);
cph2->cd(1);
TH1I *histoErrorVdrift0 = new TH1I("errorvdrift0","",100 ,-0.05,0.05);
histoErrorVdrift0->SetXTitle("#Deltav/v_{sim}");
histoErrorVdrift0->SetYTitle("counts");
histoErrorVdrift0->SetLineColor(2);
histoErrorVdrift0->SetLineStyle(2);
histoErrorVdrift0->SetStats(0);
for(Int_t k = 0; k < counter[0]; k++){
histoErrorVdrift0->Fill(yValuesDeltaPH[k]);
}
legph2->AddEntry(histoErrorVdrift0,"v_{fit PH}","l");
histoErrorVdrift0->Draw("same");
listofgraphs->Add((TObject *)histoErrorVdrift0);
}
cph2->cd(2);
legph3->Draw("same");
cph2->cd(1);
legph2->Draw("same");
//Write if wanted
if (fWriteCoef[1]){
TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
}
else{
for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
}
}
fout->Close();
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotWriteT0()
{
//
// create the graph errors and write them if wanted
//
//TObjArray of the grapherrors and so on
TObjArray *listofgraphs = new TObjArray();
Int_t nbins = fDect2[1]-fDect1[1];
//See the number of fitted for delta: here T0 can be negative, we don't use the sign but the error
//and the grapherrors of the coefficients contained the no fitted with error 0.0
//counter
Int_t counter[2];
counter[0] = 0;
counter[1] = 0;
Double_t *xValuesFitted = new Double_t[nbins];
Double_t *xValuesFittedPH = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValuesFitted[k] = -1;
xValuesFittedPH[k] = -1;
}
for(Int_t l = 0; l < nbins; l++){
if(fCoefT0E[1][l] != 0.0){
xValuesFitted[counter[1]]=l;
counter[1]++;
}
}
if(fFitPHOn){
for(Int_t l = 0; l < nbins; l++){
if(fCoefT0E[0][l] != 0.0){
xValuesFittedPH[counter[0]]= l;
counter[0]++;
}
}
}
//Create the X and Xerror
Double_t *xValues = new Double_t[nbins];
Double_t *xValuesE = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValues[k] = k;
xValuesE[k] = 0.0;
}
//Create the graph erros and plot them
TGraph *graphT02 = new TGraph(nbins,xValues,fCoefT0[2]);
graphT02->SetName("coeft02");
graphT02->SetTitle("");
graphT02->GetXaxis()->SetTitle("Det/Pad groups");
graphT02->GetYaxis()->SetTitle("T0 [time bins]");
graphT02->SetLineColor(4);
listofgraphs->Add((TObject *)graphT02);
TGraphErrors *graphT01 = new TGraphErrors(nbins,xValues,fCoefT0[1],xValuesE,fCoefT0E[1]);
graphT01->SetName("coeft01");
graphT01->SetTitle("");
graphT01->GetXaxis()->SetTitle("Det/Pad groups");
graphT01->GetYaxis()->SetTitle("T0 [time bins]");
graphT01->SetMarkerColor(6);
graphT01->SetLineColor(6);
graphT01->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphT01);
TCanvas *ct01 = new TCanvas("ct01","",50,50,600,800);
ct01->cd();
TLegend *legt01 = new TLegend(0.4,0.6,0.89,0.89);
legt01->AddEntry(graphT02,"T0 simulated","l");
legt01->AddEntry(graphT01,"T0 slope method","p");
graphT02->Draw("AL");
graphT01->Draw("P");
if (fFitPHOn) {
TGraphErrors *graphT00 = new TGraphErrors(nbins,xValues,fCoefT0[0],xValuesE,fCoefT0E[0]);
graphT00->SetName("coeft00");
graphT00->GetXaxis()->SetTitle("Det/Pad groups");
graphT00->GetYaxis()->SetTitle("T0 [time bins]");
graphT00->SetTitle("");
graphT00->SetMarkerColor(2);
graphT00->SetLineColor(2);
graphT00->SetMarkerStyle(24);
legt01->AddEntry(graphT00,"T0 fit","p");
graphT00->Draw("P");
listofgraphs->Add((TObject *)graphT00);
}
legt01->Draw("same");
//Create the arrays and the graphs for the delta
TCanvas *ct02 = new TCanvas("ct02","",50,50,600,800);
ct02->Divide(2,1);
ct02->cd(2);
Double_t *yValuesDelta = new Double_t[counter[1]];
for(Int_t k = 0; k < counter[1]; k++){
yValuesDelta[k] = (fCoefT0[1][(Int_t)xValuesFitted[k]]-fCoefT0[2][(Int_t)xValuesFitted[k]]);
}
TGraph *graphDeltaT01 = new TGraph(counter[1],&xValuesFitted[0],yValuesDelta);
graphDeltaT01->SetName("deltat01");
graphDeltaT01->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaT01->GetYaxis()->SetTitle("#Deltat0 [time bins]");
graphDeltaT01->SetMarkerColor(6);
graphDeltaT01->SetTitle("");
graphDeltaT01->SetLineColor(6);
graphDeltaT01->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphDeltaT01);
TLegend *legt03 = new TLegend(0.4,0.6,0.89,0.89);
legt03->AddEntry(graphDeltaT01,"T0_{slope method}","p");
graphDeltaT01->Draw("AP");
ct02->cd(1);
TH1I *histoErrorT01 = new TH1I("errort01","",100 ,-0.05,0.05);
histoErrorT01->SetXTitle("#Deltat0 [time bins]");
histoErrorT01->SetYTitle("counts");
histoErrorT01->SetLineColor(6);
histoErrorT01->SetLineStyle(1);
histoErrorT01->SetStats(0);
for(Int_t k = 0; k < counter[1]; k++){
histoErrorT01->Fill(yValuesDelta[k]);
}
TLegend *legt02 = new TLegend(0.4,0.6,0.89,0.89);
legt02->AddEntry(histoErrorT01,"T0_{slope method}","l");
histoErrorT01->Draw();
listofgraphs->Add((TObject *)histoErrorT01);
if (fFitPHOn) {
ct02->cd(2);
Double_t *yValuesDeltaPH = new Double_t[counter[0]];
for(Int_t k = 0; k < counter[0]; k++){
yValuesDeltaPH[k] = (fCoefT0[0][(Int_t)xValuesFittedPH[k]]-fCoefT0[2][(Int_t)xValuesFittedPH[k]]);
}
TGraph *graphDeltaT00 = new TGraph(counter[0],&xValuesFittedPH[0],yValuesDeltaPH);
graphDeltaT00->SetName("deltat00");
graphDeltaT00->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaT00->GetYaxis()->SetTitle("#Deltat0 [time bins]");
graphDeltaT00->SetMarkerColor(2);
graphDeltaT00->SetMarkerStyle(24);
graphDeltaT00->SetLineColor(2);
graphDeltaT00->SetTitle("");
legt03->AddEntry(graphDeltaT00,"T0_{fit PH}","p");
graphDeltaT00->Draw("P");
listofgraphs->Add((TObject *)graphDeltaT00);
ct02->cd(1);
TH1I *histoErrorT00 = new TH1I("errort00","",100 ,-0.05,0.05);
histoErrorT00->SetXTitle("#Deltat0 [time bins]");
histoErrorT00->SetYTitle("counts");
histoErrorT00->SetLineColor(2);
histoErrorT00->SetLineStyle(2);
histoErrorT00->SetStats(0);
for(Int_t k = 0; k < counter[0]; k++){
histoErrorT00->Fill(yValuesDeltaPH[k]);
}
legt02->AddEntry(histoErrorT00,"T0_{fit PH}","l");
histoErrorT00->Draw("same");
listofgraphs->Add((TObject *)histoErrorT00);
}
ct02->cd(2);
legt03->Draw("same");
ct02->cd(1);
legt02->Draw("same");
//Write if wanted
if (fWriteCoef[1]){
TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
}
else{
for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
}
}
fout->Close();
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotWritePRF()
{
//
// create the graph errors and write them if wanted
//
//TObjArray of the grapherrors and so on
TObjArray *listofgraphs = new TObjArray();
Int_t nbins = fDect2[2]-fDect1[2];
//See the number of fitted for delta
//counter
Int_t counter = 0;
Double_t *xValuesFitted = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValuesFitted[k] = -1;
}
for(Int_t l = 0; l < nbins; l++){
if(fCoefPRF[0][l] > 0){
xValuesFitted[counter]=l;
counter++;
}
}
//Create the X and Xerror
Double_t *xValues = new Double_t[nbins];
Double_t *xValuesE = new Double_t[nbins];
for(Int_t k = 0; k < nbins; k ++){
xValues[k] = k;
xValuesE[k] = 0.0;
}
//Create the graph erros and plot them
TGraph *graphPRF1 = new TGraph(nbins,xValues,fCoefPRF[1]);
graphPRF1->SetName("coefprf1");
graphPRF1->SetTitle("");
graphPRF1->GetXaxis()->SetTitle("Det/Pad groups");
graphPRF1->GetYaxis()->SetTitle("PRF width [p.u]");
graphPRF1->SetLineColor(4);
graphPRF1->SetMarkerColor(4);
graphPRF1->SetMarkerStyle(25);
graphPRF1->SetMarkerSize(0.7);
listofgraphs->Add((TObject *)graphPRF1);
TGraphErrors *graphPRF0 = new TGraphErrors(nbins,xValues,fCoefPRF[0],xValuesE,fCoefPRFE);
graphPRF0->SetName("coefprf0");
graphPRF0->SetTitle("");
graphPRF0->GetXaxis()->SetTitle("Det/Pad groups");
graphPRF0->GetYaxis()->SetTitle("PRF Width [p.u]");
graphPRF0->SetMarkerColor(6);
graphPRF0->SetLineColor(6);
graphPRF0->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphPRF0);
TCanvas *cprf1 = new TCanvas("cprf1","",50,50,600,800);
cprf1->cd();
TLegend *legprf1 = new TLegend(0.4,0.6,0.89,0.89);
legprf1->AddEntry(graphPRF1,"PRF width simulated","p");
legprf1->AddEntry(graphPRF0,"PRF fit","p");
graphPRF1->Draw("AP");
graphPRF0->Draw("P");
legprf1->Draw("same");
//Create the arrays and the graphs for the delta
TCanvas *cprf2 = new TCanvas("cprf2","",50,50,600,800);
cprf2->Divide(2,1);
cprf2->cd(2);
Double_t *yValuesDelta = new Double_t[counter];
for(Int_t k = 0; k < counter; k++){
if(fCoefPRF[1][(Int_t)xValuesFitted[k]] > 0.0){
yValuesDelta[k] = (fCoefPRF[0][(Int_t)xValuesFitted[k]]-fCoefPRF[1][(Int_t)xValuesFitted[k]])/(fCoefPRF[1][(Int_t)xValuesFitted[k]]);
}
}
TGraph *graphDeltaPRF = new TGraph(counter,&xValuesFitted[0],yValuesDelta);
graphDeltaPRF->SetName("deltaprf");
graphDeltaPRF->GetXaxis()->SetTitle("Det/Pad groups");
graphDeltaPRF->GetYaxis()->SetTitle("#Delta#sigma/#sigma_{sim}");
graphDeltaPRF->SetMarkerColor(6);
graphDeltaPRF->SetTitle("");
graphDeltaPRF->SetLineColor(6);
graphDeltaPRF->SetMarkerStyle(26);
listofgraphs->Add((TObject *)graphDeltaPRF);
TLegend *legprf3 = new TLegend(0.4,0.6,0.89,0.89);
legprf3->AddEntry(graphDeltaPRF,"#sigma_{fit}","p");
graphDeltaPRF->Draw("AP");
cprf2->cd(1);
TH1I *histoErrorPRF = new TH1I("errorprf1","",100 ,-0.5,0.5);
histoErrorPRF->SetXTitle("#Delta#sigma/#sigma_{sim}");
histoErrorPRF->SetYTitle("counts");
histoErrorPRF->SetLineColor(6);
histoErrorPRF->SetLineStyle(1);
histoErrorPRF->SetStats(0);
for(Int_t k = 0; k < counter; k++){
histoErrorPRF->Fill(yValuesDelta[k]);
}
TLegend *legprf2 = new TLegend(0.4,0.6,0.89,0.89);
legprf2->AddEntry(histoErrorPRF,"#sigma_{fit}","l");
histoErrorPRF->Draw();
listofgraphs->Add((TObject *)histoErrorPRF);
cprf2->cd(2);
legprf3->Draw("same");
cprf2->cd(1);
legprf2->Draw("same");
//Write if wanted
if (fWriteCoef[2]){
TFile *fout = TFile::Open(fWriteNameCoef,"UPDATE");
// Check if the file could be opened
if (!fout || !fout->IsOpen()) {
AliInfo("No File found!");
}
else{
for(Int_t k = 0; k < listofgraphs->GetEntriesFast(); k++){
fout->WriteTObject((TObject *) listofgraphs->At(k),((TObject *)listofgraphs->At(k))->GetName(),(Option_t *) "OverWrite");
}
}
fout->Close();
}
}
//
//____________Plot histos DB___________________________________________________
//
//_____________________________________________________________________________
void AliTRDCalibra::PlotCHDB()
{
//
// Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
TCanvas *cchdb = new TCanvas("cchdb","",50,50,600,800);
if ((fFitChargeBisOn) && (fMeanChargeOn)) {
cchdb->Divide(3,1);
cchdb->cd(1);
fCoefChargeDB[0]->Draw("LEGO");
cchdb->cd(2);
fCoefChargeDB[1]->Draw("LEGO");
cchdb->cd(3);
fCoefChargeDB[2]->Draw("LEGO");
}
if ((!fFitChargeBisOn) && (fMeanChargeOn)) {
cchdb->Divide(2,1);
cchdb->cd(1);
fCoefChargeDB[0]->Draw("LEGO");
cchdb->cd(2);
fCoefChargeDB[1]->Draw("LEGO");
}
else {
cchdb->cd();
fCoefChargeDB[0]->Draw("LEGO");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotPHDB()
{
//
// Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
TCanvas *cphdb = new TCanvas("cphdb","",50,50,600,800);
if (fFitPHOn) {
cphdb->Divide(2,1);
cphdb->cd(1);
fCoefVdriftDB[0]->Draw("LEGO");
cphdb->cd(2);
fCoefVdriftDB[1]->Draw("LEGO");
}
else {
cphdb->cd();
fCoefVdriftDB[1]->Draw("LEGO");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotT0DB()
{
//
// Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
TCanvas *ct0db = new TCanvas("ct0db","",50,50,600,800);
if (fFitPHOn ) {
ct0db->Divide(2,1);
ct0db->cd(1);
fCoefT0DB[0]->Draw("LEGO");
ct0db->cd(2);
fCoefT0DB[1]->Draw("LEGO");
}
else {
ct0db->cd();
fCoefT0DB[1]->Draw("LEGO");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::PlotPRFDB()
{
//
// Plot the histos for fDebug = 3 and fDebug = 4 to visualise the detector
//
TCanvas *cprfdb = new TCanvas("cprfdb","",50,50,600,800);
cprfdb->cd();
fCoefPRFDB->Draw("LEGO");
}
//
//____________Write DB Histos__________________________________________________
//
//_____________________________________________________________________________
void AliTRDCalibra::WriteCHDB(TFile *fout)
{
//
// If wanted, write the debug histos for fDebug = 3 and fDebug = 4
//
fout->WriteTObject(fCoefChargeDB[0],fCoefChargeDB[0]->GetName(),(Option_t *) "OverWrite");
if (fMeanChargeOn) {
fout->WriteTObject(fCoefChargeDB[1],fCoefChargeDB[1]->GetName(),(Option_t *) "OverWrite");
}
if (fFitChargeBisOn ) {
fout->WriteTObject(fCoefChargeDB[2],fCoefChargeDB[2]->GetName(),(Option_t *) "OverWrite");
}
}
//_____________________________________________________________________________
void AliTRDCalibra::WritePHDB(TFile *fout)
{
//
// If wanted, write the debug histos for fDebug = 3 and fDebug = 4
//
if (fFitPHOn) {
fout->WriteTObject(fCoefVdriftDB[0],fCoefVdriftDB[0]->GetName(),(Option_t *) "OverWrite");
}
fout->WriteTObject(fCoefVdriftDB[1],fCoefVdriftDB[1]->GetName(),(Option_t *) "OverWrite");
}
//_____________________________________________________________________________
void AliTRDCalibra::WriteT0DB(TFile *fout)
{
//
// If wanted, write the debug histos for fDebug = 3 and fDebug = 4
//
if (fFitPHOn) {
fout->WriteTObject(fCoefT0DB[0],fCoefT0DB[0]->GetName(),(Option_t *) "OverWrite");
}
fout->WriteTObject(fCoefT0DB[1],fCoefT0DB[1]->GetName(),(Option_t *) "OverWrite");
}
//_____________________________________________________________________________
void AliTRDCalibra::WritePRFDB(TFile *fout)
{
//
// If wanted, write the debug histos for fDebug = 3 and fDebug = 4
//
fout->WriteTObject(fCoefPRFDB,fCoefPRFDB->GetName(),(Option_t *) "OverWrite");
}
//
//____________Calcul Coef Mean_________________________________________________
//
//_____________________________________________________________________________
Bool_t AliTRDCalibra::CalculT0CoefMean(Int_t dect, Int_t idect)
{
//
// For the detector Dect calcul the mean time 0
// for the calibration group idect from the choosen database
//
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB Manager");
return kFALSE;
}
fT0Coef[2] = 0.0;
if ((fDebug != 2) && fAccCDB) {
for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
// Groups of pads
if ((fNz[1] > 0) &&
(fNrphi[1] > 0)) {
fT0Coef[2] += (Float_t) cal->GetT0(dect,col,row);
}
// Per detectors
else {
fT0Coef[2] += (Float_t) cal->GetT0Average(dect);
}
}
}
fT0Coef[2] = fT0Coef[2] / ((fColMax[1]-fColMin[1])*(fRowMax[1]-fRowMin[1]));
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefT0[2][idect] = fT0Coef[2];
}
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::CalculChargeCoefMean(Int_t dect, Int_t idect, Bool_t vrai)
{
//
// For the detector Dect calcul the mean gain factor
// for the calibration group idect from the choosen database
//
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB Manager");
return kFALSE;
}
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam Manager");
return kFALSE;
}
fChargeCoef[3] = 0.0;
if (fDebug != 2) {
for (Int_t row = fRowMin[0]; row < fRowMax[0]; row++) {
for (Int_t col = fColMin[0]; col < fColMax[0]; col++) {
// Groups of pads
if ((fNz[0] > 0) ||
(fNrphi[0] > 0)) {
if (fAccCDB) {
fChargeCoef[3] += (Float_t) cal->GetGainFactor(dect,col,row);
}
if (vrai && fAccCDB) {
fScaleFitFactor += (Float_t) cal->GetGainFactor(dect,col,row);
}
if (!fAccCDB) {
fChargeCoef[3] += 1.0;
}
if (vrai && (!fAccCDB)) {
fScaleFitFactor += 1.0;
}
}
// Per detectors
else {
if (fAccCDB) {
fChargeCoef[3] += (Float_t) cal->GetGainFactorAverage(dect);
}
if (vrai && fAccCDB) {
fScaleFitFactor += ((Float_t) cal->GetGainFactorAverage(dect));
}
if (!fAccCDB) {
fChargeCoef[3] += 1.0;
}
if (vrai && (!fAccCDB)) {
fScaleFitFactor += 1.0;
}
}
}
}
fChargeCoef[3] = fChargeCoef[3] / ((fColMax[0]-fColMin[0])*(fRowMax[0]-fRowMin[0]));
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefCharge[3][idect]=fChargeCoef[3];
}
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::CalculPRFCoefMean(Int_t dect, Int_t idect)
{
//
// For the detector Dect calcul the mean sigma of pad response
// function for the calibration group idect from the choosen database
//
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB Manager");
return kFALSE;
}
AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
if (!parCom) {
AliInfo("Could not get CommonParam Manager");
return kFALSE;
}
fPRFCoef[1] = 0.0;
Int_t cot = 0;
if (fDebug != 2) {
for (Int_t row = fRowMin[2]; row < fRowMax[2]; row++) {
for (Int_t col = fColMin[2]; col < fColMax[2]; col++) {
if ((parCom->GetColMax(GetPlane(dect)) != (col+1)) && (col != 0)) {
cot++;
if (fAccCDB) {
fPRFCoef[1] += (Float_t) cal->GetPRFWidth(dect,col,row);
}
if (!fAccCDB) {
fPRFCoef[1] += GetPRFDefault(GetPlane(dect));
}
}
}
}
if (cot > 0) {
fPRFCoef[1] = fPRFCoef[1]/cot;
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefPRF[1][idect] = fPRFCoef[1];
}
}
if (cot <= 0) {
if ((fDebug == 1) ||
(fDebug == 4)) {
if (fAccCDB) {
fCoefPRF[1][idect] = cal->GetPRFWidth(dect,fColMin[2],fRowMin[2]);
}
if (!fAccCDB) {
fCoefPRF[1][idect] = GetPRFDefault(GetPlane(dect));
}
}
}
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::CalculVdriftCoefMean(Int_t dect, Int_t idect)
{
//
// For the detector dect calcul the mean drift velocity for the
// calibration group idect from the choosen database
//
AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
if (!cal) {
AliInfo("Could not get calibDB Manager");
return kFALSE;
}
fVdriftCoef[2] = 0.0;
if (fDebug != 2) {
for (Int_t row = fRowMin[1]; row < fRowMax[1]; row++) {
for (Int_t col = fColMin[1]; col < fColMax[1]; col++) {
// Groups of pads
if ((fNz[1] > 0) ||
(fNrphi[1] > 0)) {
if (fAccCDB) {
fVdriftCoef[2] += (Float_t) cal->GetVdrift(dect,col,row);
}
if (!fAccCDB) {
fVdriftCoef[2] += 1.5;
}
}
// Per detectors
else {
if (fAccCDB) {
fVdriftCoef[2] += (Float_t) cal->GetVdriftAverage(dect);
}
if (!fAccCDB) {
fVdriftCoef[2] += 1.5;
}
}
}
}
fVdriftCoef[2] = fVdriftCoef[2] / ((fColMax[1]-fColMin[1])*(fRowMax[1]-fRowMin[1]));
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefVdrift[2][idect] = fVdriftCoef[2];
}
}
return kTRUE;
}
//_____________________________________________________________________________
Float_t AliTRDCalibra::GetPRFDefault(Int_t plane) const
{
//
// Default width of the PRF if there is no database as reference
//
if (plane == 0) {
return 0.515;
}
if (plane == 1) {
return 0.502;
}
if (plane == 2) {
return 0.491;
}
if (plane == 3) {
return 0.481;
}
if (plane == 4) {
return 0.471;
}
if (plane == 5) {
return 0.463;
}
else {
return 0.0;
}
}
//
//____________Pad group calibration mode_______________________________________
//
//_____________________________________________________________________________
void AliTRDCalibra::ReconstructionRowPadGroup(Int_t idect, Int_t i)
{
//
// For the calibration group idect in a detector calculate the
// first and last row pad and col pad.
// The pads in the interval will have the same calibrated coefficients
//
Int_t posc = -1;
Int_t posr = -1;
fRowMin[i] = -1;
fRowMax[i] = -1;
fColMin[i] = -1;
fColMax[i] = -1;
if (fNfragZ[i] != 0) {
posc = (Int_t) idect / fNfragZ[i];
}
if (fNfragRphi[i] != 0) {
posr = (Int_t) idect % fNfragZ[i];
}
fRowMin[i] = posr * fNnZ[i];
fRowMax[i] = (posr+1) * fNnZ[i];
fColMin[i] = posc * fNnRphi[i];
fColMax[i] = (posc+1) * fNnRphi[i];
}
//_____________________________________________________________________________
void AliTRDCalibra::CalculXBins(Int_t idect, Int_t i)
{
//
// For the detector idect calcul the first Xbins
//
fXbins[i] = 0;
if (fDebug == 4) {
AliInfo(Form("detector: %d", idect));
}
// In which sector?
Int_t sector = GetSector(idect);
fXbins[i] += sector*(6*fDetChamb2[i]+6*4*fDetChamb0[i]);
// In which chamber?
Int_t chamber = GetChamber(idect);
Int_t kc = 0;
while (kc < chamber) {
if (kc == 2) {
fXbins[i] += 6 * fDetChamb2[i];
}
else {
fXbins[i] += 6 * fDetChamb0[i];
}
kc ++;
}
// In which plane?
Int_t plane = GetPlane(idect);
if (chamber == 2) {
fXbins[i] += plane*fDetChamb2[i];
}
else {
fXbins[i] += plane*fDetChamb0[i];
}
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::SearchInVector(Int_t group, Int_t i) const
{
//
// Search if the calibration group "group" has already been
// initialised by a previous track in the vector
//
if (i == 0) {
for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
if (((AliTRDPlace *) fPlaCH->At(k))->GetPlace() == group) {
return k;
}
}
return -1;
}
if (i == 1) {
for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
if (((AliTRDPlace *) fPlaPH->At(k))->GetPlace() == group) {
return k;
}
}
return -1;
}
if (i == 2) {
for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
if (((AliTRDPlace *) fPlaPRF->At(k))->GetPlace() == group) {
return k;
}
}
return -1;
}
return -1;
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::SearchInTreeVector(TObjArray *vectorplace, Int_t group) const
{
//
// Search if the calibration group "group" is present in the tree
//
for (Int_t k = 0; k < (Int_t) vectorplace->GetEntriesFast(); k++) {
if (((AliTRDPlace *) vectorplace->At(k))->GetPlace() == group) {
return k;
}
}
return -1;
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::SearchBin(Float_t value, Int_t i) const
{
//
// Search the bin
//
Int_t reponse = 0;
Int_t fbinmin = 0;
Int_t fbinmax = (Int_t) value;
Int_t fNumberOfBin = -1;
// Charge
if (i == 0) {
fbinmax = 300;
fbinmin = 0;
fNumberOfBin = fNumberBinCharge;
}
// PRF
if (i == 2) {
fbinmax = 1;
fbinmin = -1;
fNumberOfBin = fNumberBinPRF;
}
// Return -1 if out
if ((value >= fbinmax) ||
(value < fbinmin)) {
return -1;
}
// Sinon
else {
reponse = (Int_t) ((fNumberOfBin*(value-fbinmin)) / (fbinmax-fbinmin));
}
return reponse;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::UpdateVectorCH(Int_t group, Float_t value)
{
//
// Fill the vector if a new calibration group "group" or update the
// values of the calibration group "group" if already here
//
// Search bin
Int_t bin = SearchBin(value,0);
// Out
if ((bin < 0) || (bin >= fNumberBinCharge)) {
return kFALSE;
}
// Search place
Int_t place = SearchInVector(group,0);
// New group
if (place == -1) {
AliTRDPlace *placegroup = new AliTRDPlace();
placegroup->SetPlace(group);
fPlaCH->Add((TObject *) placegroup);
// Variable
AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
UShort_t *entries = new UShort_t[fNumberBinCharge];
// Initialise first
for(Int_t k = 0; k < fNumberBinCharge; k++) {
entries[k] = 0;
}
// Add the value
entries[bin]= 1;
// Set
fCHInfo->SetEntries(entries);
// Set in the vector
fVectorCH->Add((TObject *) fCHInfo);
}
// Group already exits
else {
// Variable
AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
// Retrieve
fCHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
UShort_t *entries = fCHInfo->GetEntries();
// Add
entries[bin]++;
// Set
fCHInfo->SetEntries(entries);
// Update the vector
fVectorCH->AddAt((TObject *) fCHInfo,place);
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::UpdateVectorPRF(Int_t group, Float_t x, Float_t y)
{
//
// Fill the vector if a new calibration group "group" or update the
// values of the calibration group "group" if already here
//
// Search bin
Int_t bin = SearchBin(x,2);
// Out
if ((bin < 0) || (bin >= fNumberBinPRF)) {
return kFALSE;
}
// Search place
Int_t place = SearchInVector(group,2);
// New group
if (place == -1) {
AliTRDPlace *placegroup = new AliTRDPlace();
placegroup->SetPlace(group);
fPlaPRF->Add((TObject *) placegroup);
AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
Float_t *sum = new Float_t[fNumberBinPRF];
Float_t *sumsquare = new Float_t[fNumberBinPRF];
UShort_t *entries = new UShort_t[fNumberBinPRF];
// Initialise first
for (Int_t k = 0; k < fNumberBinPRF; k++) {
sum[k] = 0.0;
sumsquare[k] = 0.0;
entries[k] = 0;
}
// Add the value
sum[bin] += y;
sumsquare[bin] += y*y;
entries[bin]++;
// Set
fPRFInfo->SetSum(sum);
fPRFInfo->SetSumSquare(sumsquare);
fPRFInfo->SetEntries(entries);
// Set in the vector
fVectorPRF->Add((TObject *) fPRFInfo);
}
// Group already exits
else {
AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
// Retrieve
fPRFInfo = (AliTRDPInfo *) fVectorPRF->At(place);
Float_t *sum = fPRFInfo->GetSum();
Float_t *sumsquare = fPRFInfo->GetSumSquare();
UShort_t *entries = fPRFInfo->GetEntries();
// Add
Double_t calcul = (((Double_t) fPRFInfo->GetEntries()[bin])
* ((Double_t) fPRFInfo->GetSum()[bin]) + (Double_t) y)
/ (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
sum[bin] = (Float_t) calcul;
Double_t calculsquare = (((Double_t) fPRFInfo->GetSumSquare()[bin])
* ((Double_t) fPRFInfo->GetEntries()[bin]) + ((Double_t) y)*((Double_t) y))
/ (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
sumsquare[bin] = (Float_t) calculsquare;
entries[bin]++;
// Set
fPRFInfo->SetSum(sum);
fPRFInfo->SetSumSquare(sumsquare);
fPRFInfo->SetEntries(entries);
// Update the vector
fVectorPRF->AddAt((TObject *) fPRFInfo,place);
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::UpdateVectorPH(Int_t group, Int_t time, Float_t value)
{
//
// Fill the vector if a new calibration group "group" or update
// the values of the calibration group "group" if already here
//
// Search bin
Int_t bin = time;
// Out
if ((bin < 0) ||
(bin >= fTimeMax)) {
return kFALSE;
}
// Search place
Int_t place = SearchInVector(group,1);
// New group
if(place == -1){
AliTRDPlace *placegroup = new AliTRDPlace();
placegroup->SetPlace(group);
fPlaPH->Add((TObject *) placegroup);
AliTRDPInfo *fPHInfo = new AliTRDPInfo();
Float_t *sum = new Float_t[fTimeMax];
Float_t *sumsquare = new Float_t[fTimeMax];
UShort_t *entries = new UShort_t[fTimeMax];
// Initialise first
for (Int_t k = 0; k < fTimeMax; k++) {
sum[k] = 0.0;
sumsquare[k] = 0.0;
entries[k] = 0;
}
// Add the value
sum[bin] += value;
sumsquare[bin] += value*value;
entries[bin]++;
// Set
fPHInfo->SetSum(sum);
fPHInfo->SetSumSquare(sumsquare);
fPHInfo->SetEntries(entries);
// Set in the vector
fVectorPH->Add((TObject *) fPHInfo);
}
// Group already exits
else {
AliTRDPInfo *fPHInfo = new AliTRDPInfo();
// Retrieve
fPHInfo = (AliTRDPInfo *) fVectorPH->At(place);
Float_t *sum = fPHInfo->GetSum();
Float_t *sumsquare = fPHInfo->GetSumSquare();
UShort_t *entries = fPHInfo->GetEntries();
// Add
Double_t calcul = (((Double_t) fPHInfo->GetEntries()[bin])
* ((Double_t) fPHInfo->GetSum()[bin]) + (Double_t) value)
/ (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
sum[bin] = (Float_t) calcul;
Double_t calculsquare = ((((Double_t) fPHInfo->GetSumSquare()[bin])
* ((Double_t) fPHInfo->GetEntries()[bin]))
+ (((Double_t) value) * ((Double_t)value)))
/ (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
sumsquare[bin] = (Float_t) calculsquare;
entries[bin]++;
// Set
fPHInfo->SetSum(sum);
fPHInfo->SetSumSquare(sumsquare);
fPHInfo->SetEntries(entries);
// Update the vector
fVectorPH->AddAt((TObject *) fPHInfo,place);
}
return kTRUE;
}
//_____________________________________________________________________________
TGraphErrors *AliTRDCalibra::ConvertVectorPHisto(AliTRDPInfo *pInfo
, const Char_t *name) const
{
//
// Convert the PInfo in a 1D grapherror, name must contains "PRF"
// if PRF calibration and not "PRF" for Vdrift calibration
//
TGraphErrors *histo;
const Char_t *pattern1 = "PRF";
// Axis
Double_t *x;
Double_t *y;
Double_t *ex;
Double_t *ey;
Double_t step = 0.0;
Double_t min = 0.0;
// Ntimes
Int_t ntimes = 0;
if (strstr(name,pattern1)) {
ntimes = fNumberBinPRF;
}
else {
ntimes = fTimeMax;
}
x = new Double_t[ntimes]; // Xaxis
y = new Double_t[ntimes]; // Mean
ex = new Double_t[ntimes]; // Nentries
ey = new Double_t[ntimes]; // Sum of square/nentries
// Init histo
if (!strstr(name,pattern1)) {
step = 1.0 / fSf;
min = 0.0;
}
else {
step = (1.0 - (-1.0)) / fNumberBinPRF;
min = -1.0 + step / 2.0;
}
// Fill histo
for (Int_t k = 0; k < ntimes; k++) {
x[k] = min + k*step;
y[k] = 0.0;
ex[k] = 0.0;
ey[k] = 0.0;
// Fill only if there is more than 0 something
if (pInfo->GetEntries()[k] > 0) {
ex[k] = pInfo->GetEntries()[k];
y[k] = pInfo->GetSum()[k];
ey[k] = pInfo->GetSumSquare()[k];
}
}
// Define the TGraphErrors
histo = new TGraphErrors(ntimes,x,y,ex,ey);
histo->SetTitle(name);
return histo;
}
//_____________________________________________________________________________
TH1F *AliTRDCalibra::ConvertVectorCTHisto(AliTRDCTInfo *cTInfo
, const Char_t * name) const
{
//
// Convert the CTInfo in a 1D histo
//
TH1F *histo;
Int_t ntimes = fNumberBinCharge;
UShort_t *entries = cTInfo->GetEntries();
// Init histo
histo = new TH1F(name,name,fNumberBinCharge,0,300);
histo->Sumw2();
// Fill histo
for (Int_t k = 0; k < ntimes; k++) {
histo->SetBinContent(k+1,entries[k]);
histo->SetBinError(k+1,TMath::Sqrt(TMath::Abs(entries[k])));
}
return histo;
}
//_____________________________________________________________________________
TTree *AliTRDCalibra::ConvertVectorCTTreeHisto(TObjArray *vVectorCT
, TObjArray *pPlaCT
, const Char_t *name
, const Char_t *nametitle) const
{
//
// Convert the vector in a tree with two branchs: the group number
// and the TH1F histo reconstructed from the vector
//
// Size of the things
Int_t ntotal = (Int_t) pPlaCT->GetEntriesFast();
if (ntotal == 0) {
AliInfo("nothing to write!");
TTree *treeCT = new TTree(name,nametitle);
return treeCT;
}
// Variable of the tree
Int_t groupnumber = -1; // Group calibration
TH1F *histo = 0x0;
TObjArray vectorCT = *vVectorCT;
TObjArray plaCT = *pPlaCT;
// Init the tree
TTree *treeCT = new TTree(name,nametitle);
treeCT->Branch("groupnumber",&groupnumber,"groupnumber/I");
treeCT->Branch("histo","TH1F",&histo,32000,0);
// Fill
Int_t k = 0;
while (k < ntotal) {
TString nome(name);
groupnumber = ((AliTRDPlace *) plaCT.At(0))->GetPlace();
nome += groupnumber;
histo = ConvertVectorCTHisto(((AliTRDCTInfo *) vectorCT.At(0)),nome);
treeCT->Fill();
vectorCT.RemoveAt(0);
vectorCT.Compress();
plaCT.RemoveAt(0);
plaCT.Compress();
k++;
}
return treeCT;
}
//_____________________________________________________________________________
TTree *AliTRDCalibra::ConvertVectorPTreeHisto(TObjArray *vVectorP
, TObjArray *pPlaP
, const Char_t *name
, const Char_t *nametitle) const
{
//
// Convert the vector in a tree with two branchs: the group number
// and the TGraphErrors histo reconstructed from the vector.
// The name must contain "PRF" for PRF calibration and not "PRF"
// for Vdrift calibration
//
// Size of the things
Int_t ntotal = (Int_t) pPlaP->GetEntriesFast();
if (ntotal == 0) {
AliInfo("nothing to write!");
TTree *treeP = new TTree(name,nametitle);
return treeP;
}
// Variable of the tree
Int_t groupnumber = -1; // Group calibration
TGraphErrors *histo = 0x0;
TObjArray vectorP = *vVectorP;
TObjArray plaP = *pPlaP;
// Init the tree
TTree *treeP = new TTree(name,nametitle);
treeP->Branch("groupnumber",&groupnumber,"groupnumber/I");
treeP->Branch("histo","TGraphErrors",&histo,32000,0);
// Fill
Int_t k = 0;
while (k < ntotal) {
TString nome(name);
groupnumber = ((AliTRDPlace *) plaP.At(0))->GetPlace();
nome += groupnumber;
histo = ConvertVectorPHisto((AliTRDPInfo *) vectorP.At(0),nome);
treeP->Fill();
vectorP.RemoveAt(0);
vectorP.Compress();
plaP.RemoveAt(0);
plaP.Compress();
k++;
}
return treeP;
}
//_____________________________________________________________________________
TObjArray *AliTRDCalibra::ConvertTreeVector(TTree *tree) const
{
//
// Convert the branch groupnumber of the tree taken from
// TRD.calibration.root in case of vector method in a std::vector
// to be faster
//
// Initialise
TObjArray *vectorplace = new TObjArray();
// Variable of the tree
Int_t groupnumber = -1; // Group calibration
// Set the branch
tree->SetBranchAddress("groupnumber",&groupnumber);
// Fill
Int_t ntotal = tree->GetEntries();
for (Int_t k = 0; k < ntotal; k++) {
tree->GetEntry(k);
AliTRDPlace *placegroupnumber = new AliTRDPlace();
placegroupnumber->SetPlace(groupnumber);
vectorplace->Add((TObject *) placegroupnumber);
}
return vectorplace;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::MergeVectorCT(TObjArray *vVectorCT2, TObjArray *pPlaCT2)
{
//
// Add the two vectors and place the result in the first
//
if (((Int_t) pPlaCT2->GetEntriesFast()) != ((Int_t) vVectorCT2->GetEntriesFast())){
AliInfo("VectorCT2 doesn't correspond to PlaCT2!");
return kFALSE;
}
// CH case
for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
// Look if PlaCT1[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) pPlaCT2->GetEntriesFast(); j++) {
if (((AliTRDPlace *) pPlaCT2->At(j))->GetPlace() ==
((AliTRDPlace *) fPlaCH->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the second vector nothing to do
// If in the second vector
if (place != -1) {
AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();
UShort_t *entries = new UShort_t[fNumberBinCharge];
for (Int_t nu = 0; nu < fNumberBinCharge; nu++) {
entries[nu] = ((AliTRDCTInfo *) fVectorCH->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu]
+ ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu];
}
// Set
fCTInfo->SetEntries(entries);
// Nothing to do on PlaCT1
// Update the vector
fVectorCH->AddAt((TObject *) fCTInfo,((AliTRDPlace *) fPlaCH->At(k))->GetPlace());
}
}
// And at the end the vector in CT2 but not in CH1
for (Int_t k = 0; k < (Int_t) pPlaCT2->GetEntriesFast(); k++) {
// Look if pPlaCT2[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) fPlaCH->GetEntriesFast(); j++) {
if (((AliTRDPlace *) fPlaCH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the first vector
if (place == -1) {
AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();
fCTInfo = ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()));
// Add at the end
fPlaCH->Add((TObject *) (pPlaCT2->At(k)));
fVectorCH->Add((TObject *) fCTInfo);
}
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t AliTRDCalibra::MergeVectorP(TObjArray *vVectorP2
, TObjArray *pPlaP2
, Int_t i)
{
//
// Add the two vectors and place the result in the first
//
if (((Int_t) pPlaP2->GetEntriesFast()) != ((Int_t) vVectorP2->GetEntriesFast())) {
AliInfo("VectorP2 doesn't correspond to PlaP2!");
return kFALSE;
}
// PH case
if (i == 1) {
for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
// Look if fPlaPH[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPH->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the second vector nothing to do
// If in the second vector
if (place != -1) {
AliTRDPInfo *fPInfo = new AliTRDPInfo();
UShort_t *entries = new UShort_t[fTimeMax];
Float_t *sum = new Float_t[fTimeMax];
Float_t *sumsquare = new Float_t[fTimeMax];
for (Int_t nu = 0; nu < fTimeMax; nu++) {
entries[nu] = ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]
+ ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu];
Double_t calcul = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
* ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
+ (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
* ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
/ ((Double_t) fPInfo->GetEntries()[nu]);
sum[nu] = (Float_t) calcul;
Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
* ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
+ (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
* ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
/ ((Double_t) fPInfo->GetEntries()[nu]);
sumsquare[nu] = calculsquare;
}
// Set
fPInfo->SetSum(sum);
fPInfo->SetSumSquare(sumsquare);
fPInfo->SetEntries(entries);
// Nothing to do on PlaCT1
// Update the vector VectorCT1
fVectorPH->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPH->At(k))->GetPlace());
}
}
// And at the end the vector in P2 but not in CH1
for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
// Look if PlaCT2[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) fPlaPH->GetEntriesFast(); j++) {
if (((AliTRDPlace *) fPlaPH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the first vector
if (place == -1) {
AliTRDPInfo *fPInfo = new AliTRDPInfo();
fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
// Add at the end of CH1
fPlaPH->Add(((TObject *) pPlaP2->At(k)));
fVectorPH->Add((TObject *) fPInfo);
}
}
}
// PRF case
if (i == 1) {
for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
// Look if fPlaPRF[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the second vector nothing to do
// If in the second vector
if (place != -1) {
AliTRDPInfo *fPInfo = new AliTRDPInfo();
UShort_t *entries = new UShort_t[fNumberBinPRF];
Float_t *sum = new Float_t[fNumberBinPRF];
Float_t *sumsquare = new Float_t[fNumberBinPRF];
for (Int_t nu = 0; nu < fNumberBinPRF; nu++) {
entries[nu] = ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]
+ ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu];
Double_t calcul = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
* ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
+ (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
* ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
/ ((Double_t) fPInfo->GetEntries()[nu]);
sum[nu] = (Float_t) calcul;
Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
* ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
+ (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
* ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
/ ((Double_t) fPInfo->GetEntries()[nu]);
sumsquare[nu] = calculsquare;
}
// Set
fPInfo->SetSum(sum);
fPInfo->SetSumSquare(sumsquare);
fPInfo->SetEntries(entries);
// Nothing to do on PlaCT1
// Update the vector VectorCT1
fVectorPRF->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPRF->At(k))->GetPlace());
}
}
// And at the end the vector in P2 but not in CH1
for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
// Look if PlaCT2[k] it is also in the second vector
Int_t place = -1;
for (Int_t j = 0; j < (Int_t) fPlaPRF->GetEntriesFast(); j++) {
if (((AliTRDPlace *) fPlaPRF->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
place = j;
break;
}
}
// If not in the first vector
if (place == -1) {
AliTRDPInfo *fPInfo = new AliTRDPInfo();
fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
// Add at the end of CH1
fPlaPRF->Add(((TObject *) pPlaP2->At(k)));
fVectorPRF->Add((TObject *) fPInfo);
}
}
}
return kTRUE;
}
//____________Fit Methods______________________________________________________
//_____________________________________________________________________________
void AliTRDCalibra::FitPente(TH1* projPH, Int_t idect)
{
//
// Slope methode for the drift velocity
//
// Constants
const Float_t kDrWidth = AliTRDgeometry::DrThick();
Int_t binmax = 0;
Int_t binmin = 0;
fPhd[0] = 0.0;
fPhd[1] = 0.0;
fPhd[2] = 0.0;
Int_t ju = 0;
Double_t vdriftCoefE = 0.0;
Double_t t0CoefE = 0.0;
fVdriftCoef[1] = 0.0;
fT0Coef[1] = 0.0;
TLine *line = new TLine();
// Some variables
TAxis *xpph = projPH->GetXaxis();
Int_t nbins = xpph->GetNbins();
Double_t lowedge = xpph->GetBinLowEdge(1);
Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
Double_t widbins = (upedge - lowedge) / nbins;
Double_t limit = upedge + 0.5 * widbins;
// Beginning of the signal
TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
}
binmax = (Int_t) pentea->GetMaximumBin();
if (binmax == 1) {
binmax = 2;
AliInfo("Put the binmax from 1 to 2 to enable the fit");
}
pentea->Fit("pol2","0MR","",TMath::Max(pentea->GetBinCenter(binmax-1),0.0),pentea->GetBinCenter(binmax+1));
Float_t l3P1am = pentea->GetFunction("pol2")->GetParameter(1);
Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
if (l3P2am != 0) {
fPhd[0] = -(l3P1am / (2 * l3P2am));
}
if((l3P1am != 0.0) && (l3P2am != 0.0)){
t0CoefE = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
}
// Amplification region
binmax = 0;
ju = 0;
for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
if (((projPH->GetBinContent(kbin+1) - projPH->GetBinContent(kbin)) <= 0.0) && (ju == 0) && (kbin > (fPhd[0]/widbins))) {
binmax = kbin;
ju = 1;
}
}
if (binmax == 1) {
binmax = 2;
AliInfo("Put the binmax from 1 to 2 to enable the fit");
}
projPH->Fit("pol2","0MR","",TMath::Max(projPH->GetBinCenter(binmax-1),0.0),projPH->GetBinCenter(binmax+1));
Float_t l3P1amf = projPH->GetFunction("pol2")->GetParameter(1);
Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
if (l3P2amf != 0) {
fPhd[1] = -(l3P1amf / (2 * l3P2amf));
}
if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
vdriftCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
}
// Drift region
TH1D *pente = new TH1D("pente","pente",projPH->GetNbinsX(),0,(Float_t) limit);
for (Int_t k = binmax+4; k < projPH->GetNbinsX(); k++) {
pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
}
binmin = (Int_t) pente->GetMinimumBin();
if (binmin == 1) {
binmin = 2;
AliInfo("Put the binmax from 1 to 2 to enable the fit");
}
pente->Fit("pol2"
,"0MR"
,""
,TMath::Max(pente->GetBinCenter(binmin-1), 0.0)
,TMath::Min(pente->GetBinCenter(binmin+2),(Double_t) limit));
Float_t l3P1dr = pente->GetFunction("pol2")->GetParameter(1);
Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
if (l3P2dr != 0) {
fPhd[2] = -(l3P1dr / (2 * l3P2dr));
}
if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
vdriftCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
}
if ((fPhd[2] > fPhd[0]) &&
(fPhd[2] > fPhd[1]) &&
(fPhd[1] > fPhd[0])) {
fVdriftCoef[1] = (kDrWidth) / (fPhd[2]-fPhd[1]);
if (fPhd[0] >= 0.0) {
fT0Coef[1] = (fPhd[0] - fT0Shift) / widbins;
if (fT0Coef[1] < -1.0) {
fT0Coef[1] = fT0Coef[2];
}
}
else {
fT0Coef[1] = fT0Coef[2];
}
}
else {
fVdriftCoef[1] = -TMath::Abs(fVdriftCoef[2]);
fT0Coef[1] = fT0Coef[2];
}
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefVdrift[1][idect] = fVdriftCoef[1];
fCoefVdriftE[1] [idect] = vdriftCoefE;
fCoefT0[1][idect] = fT0Coef[1];
fCoefT0E[1][idect] = t0CoefE;
}
if (fDebug == 2) {
TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
cpentei->cd();
projPH->Draw();
line->SetLineColor(2);
line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
line->DrawLine(fPhd[2],0,fPhd[2],projPH->GetMaximum());
AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
AliInfo(Form("fVriftCoef[1] (with only the drift region(default)): %f",(Float_t) fVdriftCoef[1]));
}
if (fDebug != 2) {
delete pentea;
}
if (fDebug != 2) {
delete pente;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FitPH(TH1* projPH, Int_t idect)
{
//
// Fit methode for the drift velocity
//
// Constants
const Float_t kDrWidth = AliTRDgeometry::DrThick();
// Some variables
TAxis *xpph = projPH->GetXaxis();
Double_t upedge = xpph->GetBinUpEdge(xpph->GetNbins());
TF1 *fPH = new TF1("fPH",AliTRDCalibra::PH,-0.05,3.2,6);
fPH->SetParameter(0,0.469); // Scaling
fPH->SetParameter(1,0.18); // Start
fPH->SetParameter(2,0.0857325); // AR
fPH->SetParameter(3,1.89); // DR
fPH->SetParameter(4,0.08); // QA/QD
fPH->SetParameter(5,0.0); // Baseline
TLine *line = new TLine();
fVdriftCoef[0] = 0.0;
fT0Coef[0] = 0.0;
Double_t vdriftCoefE = 0.0;
Double_t t0CoefE = 0.0;
if (idect%fFitPHPeriode == 0) {
AliInfo(Form(" The detector %d will be fitted",idect));
fPH->SetParameter(0,(projPH->Integral()-(projPH->GetBinContent(1)*projPH->GetNbinsX())) * 0.00028); // Scaling
fPH->SetParameter(1,fPhd[0] - 0.1); // Start
fPH->SetParameter(2,fPhd[1] - fPhd[0]); // AR
fPH->SetParameter(3,fPhd[2] - fPhd[1]); // DR
fPH->SetParameter(4,0.225); // QA/QD
fPH->SetParameter(5,(Float_t) projPH->GetBinContent(1));
if (fDebug != 2) {
projPH->Fit(fPH,"0M","",0.0,upedge);
}
if (fDebug == 2) {
TCanvas *cpente = new TCanvas("cpente","cpente",50,50,600,800);
cpente->cd();
projPH->Fit(fPH,"M+","",0.0,upedge);
projPH->Draw("E0");
line->SetLineColor(4);
line->DrawLine(fPH->GetParameter(1)
,0
,fPH->GetParameter(1)
,projPH->GetMaximum());
line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)
,0
,fPH->GetParameter(1)+fPH->GetParameter(2)
,projPH->GetMaximum());
line->DrawLine(fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
,0
,fPH->GetParameter(1)+fPH->GetParameter(2)+fPH->GetParameter(3)
,projPH->GetMaximum());
}
if (fPH->GetParameter(3) != 0) {
fVdriftCoef[0] = kDrWidth / (fPH->GetParameter(3));
vdriftCoefE = (fPH->GetParError(3)/fPH->GetParameter(3))*fVdriftCoef[0];
fT0Coef[0] = fPH->GetParameter(1);
t0CoefE = fPH->GetParError(1);
}
else {
fVdriftCoef[0] = -TMath::Abs(fVdriftCoef[2]);
fT0Coef[0] = fT0Coef[2];
}
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefVdrift[0][idect] = fVdriftCoef[0];
fCoefVdriftE[0][idect] = vdriftCoefE;
fCoefT0[0][idect] = fT0Coef[0];
fCoefT0E[0][idect] = t0CoefE;
}
if (fDebug == 2) {
AliInfo(Form("fVdriftCoef[0]: %f",(Float_t) fVdriftCoef[0]));
}
}
else {
// Put the default value
if ((fDebug <= 1) ||
(fDebug == 4)) {
fCoefVdrift[0][idect] = -TMath::Abs(fVdriftCoef[2]);
fCoefT0[0][idect] = -TMath::Abs(fT0Coef[2]);
}
}
if (fDebug != 2) {
delete fPH;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FitPRF(TH1 *projPRF, Int_t idect)
{
//
// Fit methode for the sigma of the pad response function
//
fPRFCoef[0] = 0.0;
Double_t prfCoefE = 0.0;
if (fDebug != 2) {
projPRF->Fit("gaus","0M","",-fRangeFitPRF,fRangeFitPRF);
}
if (fDebug == 2) {
TCanvas *cfit = new TCanvas("cfit","cfit",50,50,600,800);
cfit->cd();
projPRF->Fit("gaus","M+","",-fRangeFitPRF,fRangeFitPRF);
projPRF->Draw();
}
fPRFCoef[0] = projPRF->GetFunction("gaus")->GetParameter(2);
prfCoefE = projPRF->GetFunction("gaus")->GetParError(2);
if(fPRFCoef[0] <= 0.0) fPRFCoef[0] = -fPRFCoef[1];
if ((fDebug == 1) ||
(fDebug == 4)) {
fCoefPRF[0][idect] = fPRFCoef[0];
fCoefPRFE[idect] = prfCoefE;
}
if (fDebug == 2) {
AliInfo(Form("fPRFCoef[0]: %f",(Float_t) fPRFCoef[0]));
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FitCH(TH1 *projch, Int_t idect)
{
//
// Fit methode for the gain factor
//
fChargeCoef[0] = 0.0;
fChargeCoef[1] = 0.0;
Double_t chargeCoefE0 = 0.0;
Double_t chargeCoefE1 = 0.0;
TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
fChargeCoef[1] = projch->GetMean();
chargeCoefE1 = projch->GetMeanError();
projch->Fit("landau","0",""
,(Float_t) fChargeCoef[1]/fBeginFitCharge
,projch->GetBinCenter(projch->GetNbinsX()));
fL3P0 = projch->GetFunction("landau")->GetParameter(0);
Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
fL3P2 = projch->GetFunction("landau")->GetParameter(2);
projch->Fit("gaus","0",""
,(Float_t) fChargeCoef[1]/fBeginFitCharge
,projch->GetBinCenter(projch->GetNbinsX()));
Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
fG3P2 = projch->GetFunction("gaus")->GetParameter(2);
fLandauGaus->SetParameters(fL3P0,l3P1,fL3P2,g3P0,fG3P2);
if ((fDebug <= 1) ||
(fDebug >= 3)) {
projch->Fit("fLandauGaus","0",""
,(Float_t) fChargeCoef[1]/fBeginFitCharge
,projch->GetBinCenter(projch->GetNbinsX()));
}
if (fDebug == 2) {
TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
cp->cd();
projch->Fit("fLandauGaus","+",""
,(Float_t) fChargeCoef[1]/fBeginFitCharge
,projch->GetBinCenter(projch->GetNbinsX()));
projch->Draw();
fLandauGaus->Draw("same");
}
if (projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) {
// Calcul of "real" coef
CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kTRUE);
fChargeCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
chargeCoefE0 = projch->GetFunction("fLandauGaus")->GetParError(1);
}
else {
// Calcul of "real" coef
CalculChargeCoefMean(fCountDet[0],(Int_t) idect,kFALSE);
fChargeCoef[0] = -TMath::Abs(fChargeCoef[3]);
}
if (fDebug == 2) {
AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fChargeCoef[0]));
AliInfo(Form("fChargeCoef[1]: %f",(Float_t) fChargeCoef[1]));
}
if ((fDebug == 1) ||
(fDebug == 4)) {
if (fChargeCoef[0] > 0.0) {
fCoefCharge[0][idect]= fChargeCoef[0];
fCoefChargeE[0][idect]= chargeCoefE0;
fCoefCharge[1][idect]= fChargeCoef[1];
fCoefChargeE[1][idect]= chargeCoefE1;
}
}
fL3P0 = fLandauGaus->Integral(0.3*projch->GetMean(),3*projch->GetMean());
fG3P2 = fLandauGaus->GetParameter(2);
fL3P2 = fLandauGaus->GetParameter(4);
if (fDebug != 2) {
delete fLandauGaus;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FitBisCH(TH1* projch, Int_t idect)
{
//
// Fit methode for the gain factor more time consuming
//
// Setting fit range and start values
Double_t fr[2];
//Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
Double_t pllo[4] = { 0.001, 0.001, 0.001, 0.001 };
Double_t plhi[4] = { 300.0, 300.0, 100000000.0, 300.0 };
Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
fr[0] = 0.3 * projch->GetMean();
fr[1] = 3.0 * projch->GetMean();
fChargeCoef[2] = 0.0;
Double_t chargeCoefE2 = 0.0;
Double_t chisqr;
Int_t ndf;
TF1 *fitsnr = LanGauFit(projch,&fr[0],&sv[0]
,&pllo[0],&plhi[0]
,&fp[0],&fpe[0]
,&chisqr,&ndf);
Double_t projchPeak;
Double_t projchFWHM;
LanGauPro(fp,projchPeak,projchFWHM);
if (fp[1] > 0) {
fChargeCoef[2] = fp[1];
chargeCoefE2 = fpe[1];
//chargeCoefE2 = chisqr;
}
else {
fChargeCoef[2] = -TMath::Abs(fChargeCoef[3]);
}
if (fDebug == 2) {
AliInfo(Form("fChargeCoef[2]: %f",(Float_t) fChargeCoef[2]));
TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
cpy->cd();
projch->Draw();
fitsnr->Draw("same");
}
if ((fDebug == 1) ||
(fDebug == 4)) {
if (fChargeCoef[2] > 0.0) {
fCoefCharge[2][idect]= fChargeCoef[2];
fCoefChargeE[2][idect]= chargeCoefE2;
}
}
if (fDebug != 2) {
delete fitsnr;
}
}
//_____________________________________________________________________________
void AliTRDCalibra::NormierungCharge()
{
//
// Normalisation of the gain factor resulting for the fits
//
// Calcul of the mean of the fit
Double_t sum = 0.0;
for (Int_t k = 0; k < (Int_t) fVectorFitCH->GetEntriesFast(); k++) {
Int_t total = 0;
Int_t detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
Float_t *coef = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef();
if (GetChamber(detector) == 2) {
total = 1728;
}
if (GetChamber(detector) != 2) {
total = 2304;
}
for (Int_t j = 0; j < total; j++) {
if (coef[j] >= 0) {
sum += coef[j];
}
}
}
if (sum > 0) {
fScaleFitFactor = fScaleFitFactor / sum;
}
else {
fScaleFitFactor = 1.0;
}
if ((fDebug == 3) ||
(fDebug == 4)) {
if ((fCoefChargeDB[0]->GetEntries() > 0.0) &&
(fCoefChargeDB[0]->GetSumOfWeights() > 0.0)) {
fCoefChargeDB[0]->Scale(fCoefChargeDB[0]->GetEntries() / fCoefChargeDB[0]->GetSumOfWeights());
}
if ((fMeanChargeOn) &&
(fCoefChargeDB[1]->GetEntries() > 0.0) &&
(fCoefChargeDB[1]->GetSumOfWeights() > 0.0)) {
fCoefChargeDB[1]->Scale(fCoefChargeDB[1]->GetEntries() / fCoefChargeDB[1]->GetSumOfWeights());
}
if ((fFitChargeBisOn) &&
(fCoefChargeDB[2]->GetEntries() > 0.0) &&
(fCoefChargeDB[2]->GetSumOfWeights() > 0.0)) {
fCoefChargeDB[2]->Scale(fCoefChargeDB[2]->GetEntries() / fCoefChargeDB[2]->GetSumOfWeights());
}
}
}
//_____________________________________________________________________________
TH1I *AliTRDCalibra::ReBin(TH1I *hist) const
{
//
// Rebin of the 1D histo for the gain calibration if needed.
// you have to choose fRebin, divider of fNumberBinCharge
//
TAxis *xhist = hist->GetXaxis();
TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
,xhist->GetBinLowEdge(1)
,xhist->GetBinUpEdge(xhist->GetNbins()));
AliInfo(Form("fRebin: %d",fRebin));
Int_t i = 1;
for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
Double_t sum = 0.0;
for (Int_t ji = i; ji < i+fRebin; ji++) {
sum += hist->GetBinContent(ji);
}
sum = sum / fRebin;
rehist->SetBinContent(k,sum);
i += fRebin;
}
if (fDebug == 2) {
TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
crebin->cd();
rehist->Draw();
}
return rehist;
}
//_____________________________________________________________________________
TH1F *AliTRDCalibra::ReBin(TH1F *hist) const
{
//
// Rebin of the 1D histo for the gain calibration if needed
// you have to choose fRebin divider of fNumberBinCharge
//
TAxis *xhist = hist->GetXaxis();
TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
,xhist->GetBinLowEdge(1)
,xhist->GetBinUpEdge(xhist->GetNbins()));
AliInfo(Form("fRebin: %d",fRebin));
Int_t i = 1;
for (Int_t k = 1; k <= (Int_t) xhist->GetNbins()/fRebin; k++) {
Double_t sum = 0.0;
for (Int_t ji = i; ji < i+fRebin; ji++) {
sum += hist->GetBinContent(ji);
}
sum = sum/fRebin;
rehist->SetBinContent(k,sum);
i += fRebin;
}
if (fDebug == 2) {
TCanvas *crebin = new TCanvas("crebin","",50,50,600,800);
crebin->cd();
rehist->Draw();
}
return rehist;
}
//_____________________________________________________________________________
TH1F *AliTRDCalibra::CorrectTheError(TGraphErrors *hist)
{
//
// In the case of the vectors method the trees contains TGraphErrors for PH and PRF
// to be able to add them after
// We convert it to a TH1F to be able to applied the same fit function method
// After having called this function you can not add the statistics anymore
//
TH1F *rehist = 0x0;
Int_t nbins = hist->GetN();
Double_t *x = hist->GetX();
Double_t *entries = hist->GetEX();
Double_t *mean = hist->GetY();
Double_t *square = hist->GetEY();
fEntriesCurrent = 0;
if (nbins < 2) {
return rehist;
}
Double_t step = x[1] - x[0];
Double_t minvalue = x[0] - step/2;
Double_t maxvalue = x[(nbins-1)] + step/2;
rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue);
for (Int_t k = 0; k < nbins; k++) {
rehist->SetBinContent(k+1,mean[k]);
if (entries[k] > 0.0) {
fEntriesCurrent += (Int_t) entries[k];
Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k]));
rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k]));
}
else {
rehist->SetBinError(k+1,0.0);
}
}
return rehist;
}
//_____________________________________________________________________________
TGraphErrors *AliTRDCalibra::AddProfiles(TGraphErrors *hist1
, TGraphErrors *hist2) const
{
//
// In the case of the vectors method we use TGraphErrors for PH and PRF
// to be able to add the them after
// Here we add the TGraphErrors
//
// First TGraphErrors
Int_t nbins1 = hist1->GetN();
Double_t *x1 = hist1->GetX();
Double_t *ex1 = hist1->GetEX();
Double_t *y1 = hist1->GetY();
Double_t *ey1 = hist1->GetEY();
TGraphErrors *rehist = new TGraphErrors(nbins1);
// Second TGraphErrors
Double_t *ex2 = hist2->GetEX();
Double_t *y2 = hist2->GetY();
Double_t *ey2 = hist2->GetEY();
// Define the Variables for the new TGraphErrors
Double_t x;
Double_t ex;
Double_t y;
Double_t ey;
for (Int_t k = 0; k < nbins1; k++) {
Double_t nentries = 0.0;
x = x1[k];
y = 0.0;
ey = 0.0;
ex = 0.0;
if ((ex2[k] == 0.0) &&
(ex1[k] == 0.0)) {
nentries = 0.0;
}
if ((ex2[k] == 0.0) &&
(ex1[k] > 0.0)) {
nentries = ex1[k];
y = y1[k];
ey = ey1[k];
ex = ex1[k];
}
if ((ex2[k] > 0.0) &&
(ex1[k] == 0.0)) {
nentries = ex2[k];
y = y2[k];
ey = ey2[k];
ex = ex2[k];
}
if ((ex2[k] > 0.0) &&
(ex1[k] > 0.0)) {
nentries = ex1[k] + ex2[k];
y = ( y1[k]*ex1[k]+ y2[k]*ex2[k]) / nentries;
ey = (ey1[k]*ex1[k]+ey2[k]*ex2[k]) / nentries;
ex = nentries;
}
rehist->SetPoint(k,x,y);
rehist->SetPointError(k,ex,ey);
}
return rehist;
}
//
//____________Some basic geometry function_____________________________________
//
//_____________________________________________________________________________
Int_t AliTRDCalibra::GetPlane(Int_t d) const
{
//
// Reconstruct the plane number from the detector number
//
return ((Int_t) (d % 6));
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::GetChamber(Int_t d) const
{
//
// Reconstruct the chamber number from the detector number
//
Int_t fgkNplan = 6;
return ((Int_t) (d % 30) / fgkNplan);
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::GetSector(Int_t d) const
{
//
// Reconstruct the sector number from the detector number
//
Int_t fg = 30;
return ((Int_t) (d / fg));
}
//
//____________Fill and Init tree Gain, PRF, Vdrift and T0______________________
//
//_____________________________________________________________________________
void AliTRDCalibra::InitTreePRF()
{
//
// Init the tree where the coefficients from the fit methods can be stored
//
gDirectory = gROOT;
fPRFPad = new Float_t[2304];
fPRF = new TTree("PRF","PRF");
fPRF->Branch("detector",&fPRFDetector,"detector/I");
fPRF->Branch("width" ,fPRFPad ,"width[2304]/F");
// Set to default value for the plane 0 supposed to be the first one
for (Int_t k = 0; k < 2304; k++) {
fPRFPad[k] = 0.515;
}
fPRFDetector = -1;
}
//_____________________________________________________________________________
void AliTRDCalibra::FillTreePRF(Int_t countdet)
{
//
// Fill the tree with the sigma of the pad response function for the detector countdet
//
Int_t numberofgroup = 0;
fPRFDetector = countdet;
fPRF->Fill();
if (GetChamber((Int_t)(countdet+1)) == 2) {
numberofgroup = 1728;
}
else {
numberofgroup = 2304;
}
// Reset to default value for the next
for (Int_t k = 0; k < numberofgroup; k++) {
if (GetPlane((Int_t) (countdet+1)) == 0) {
fPRFPad[k] = 0.515;
}
if (GetPlane((Int_t) (countdet+1)) == 1) {
fPRFPad[k] = 0.502;
}
if (GetPlane((Int_t) (countdet+1)) == 2) {
fPRFPad[k] = 0.491;
}
if (GetPlane((Int_t) (countdet+1)) == 3) {
fPRFPad[k] = 0.481;
}
if (GetPlane((Int_t) (countdet+1)) == 4) {
fPRFPad[k] = 0.471;
}
if (GetPlane((Int_t) (countdet+1)) == 5) {
fPRFPad[k] = 0.463;
}
}
fPRFDetector = -1;
}
//_____________________________________________________________________________
void AliTRDCalibra::ConvertVectorFitCHTree()
{
//
// Convert the vector stuff to a tree of 1D histos if the user
// want to write it after the fill functions
//
Int_t detector = -1;
Int_t numberofgroup = 1;
Float_t gainPad[2304];
fGain = new TTree("Gain","Gain");
fGain->Branch("detector",&detector,"detector/I");
fGain->Branch("gainPad" ,gainPad ,"gainPad[2304]/F");
Int_t loop = (Int_t) fVectorFitCH->GetEntriesFast();
for (Int_t k = 0; k < loop; k++) {
detector = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector();
if (GetChamber((Int_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetDetector()) == 2) {
numberofgroup = 1728;
}
else {
numberofgroup = 2304;
}
for (Int_t i = 0; i < numberofgroup; i++) {
if (((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] >= 0) {
gainPad[i] = ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i] * fScaleFitFactor;
}
else {
gainPad[i] = (Float_t) ((AliTRDFitCHInfo *) fVectorFitCH->At(k))->GetCoef()[i];
}
}
fGain->Fill();
}
}
//_____________________________________________________________________________
void AliTRDCalibra::FillTreeVdrift(Int_t countdet)
{
//
// Fill the tree with the drift velocities for the detector countdet
//
Int_t numberofgroup = 0;
fVdriftDetector = countdet;
fVdrift->Fill();
if (GetChamber((Int_t)(countdet+1)) == 2) {
numberofgroup = 1728;
}
else {
numberofgroup = 2304;
}
// Reset to default value the gain coef
for (Int_t k = 0; k < numberofgroup; k++) {
fVdriftPad[k] = -1.5;
}
fVdriftDetector = -1;
}
//_____________________________________________________________________________
void AliTRDCalibra::InitTreePH()
{
//
// Init the tree where the coefficients from the fit methods can be stored
//
gDirectory = gROOT;
fVdriftPad = new Float_t[2304];
fVdrift = new TTree("Vdrift","Vdrift");
fVdrift->Branch("detector",&fVdriftDetector,"detector/I");
fVdrift->Branch("vdrift" ,fVdriftPad ,"vdrift[2304]/F");
// Set to default value for the plane 0 supposed to be the first one
for (Int_t k = 0; k < 2304; k++) {
fVdriftPad[k] = -1.5;
}
fVdriftDetector = -1;
}
//_____________________________________________________________________________
void AliTRDCalibra::FillTreeT0(Int_t countdet)
{
//
// Fill the tree with the t0 value for the detector countdet
//
Int_t numberofgroup = 0;
fT0Detector = countdet;
fT0->Fill();
if (GetChamber((Int_t) (countdet+1)) == 2) {
numberofgroup = 1728;
}
else {
numberofgroup = 2304;
}
// Reset to default value
for (Int_t k = 0; k < numberofgroup; k++) {
fT0Pad[k] = 0.0;
}
fT0Detector = -1;
}
//_____________________________________________________________________________
void AliTRDCalibra::InitTreeT0()
{
//
// Init the tree where the coefficients from the fit methods can be stored
//
gDirectory = gROOT;
fT0Pad = new Float_t[2304];
fT0 = new TTree("T0","T0");
fT0->Branch("detector",&fT0Detector,"detector/I");
fT0->Branch("t0",fT0Pad,"t0[2304]/F");
//Set to default value for the plane 0 supposed to be the first one
for(Int_t k = 0; k < 2304; k++){
fT0Pad[k] = 0.0;
}
fT0Detector = -1;
}
//
//____________Private Functions________________________________________________
//
//_____________________________________________________________________________
Double_t AliTRDCalibra::PH(Double_t *x, Double_t *par)
{
//
// Function for the fit
//
//TF1 *fAsymmGauss = new TF1("fAsymmGauss",AsymmGauss,0,4,6);
//PARAMETERS FOR FIT PH
// PASAv.4
//fAsymmGauss->SetParameter(0,0.113755);
//fAsymmGauss->SetParameter(1,0.350706);
//fAsymmGauss->SetParameter(2,0.0604244);
//fAsymmGauss->SetParameter(3,7.65596);
//fAsymmGauss->SetParameter(4,1.00124);
//fAsymmGauss->SetParameter(5,0.870597); // No tail cancelation
Double_t xx = x[0];
if (xx < par[1]) {
return par[5];
}
Double_t dx = 0.005;
Double_t xs = par[1];
Double_t ss = 0.0;
Double_t paras[2] = { 0.0, 0.0 };
while (xs < xx) {
if ((xs >= par[1]) &&
(xs < (par[1]+par[2]))) {
//fAsymmGauss->SetParameter(0,par[0]);
//fAsymmGauss->SetParameter(1,xs);
//ss += fAsymmGauss->Eval(xx);
paras[0] = par[0];
paras[1] = xs;
ss += AsymmGauss(&xx,paras);
}
if ((xs >= (par[1]+par[2])) &&
(xs < (par[1]+par[2]+par[3]))) {
//fAsymmGauss->SetParameter(0,par[0]*par[4]);
//fAsymmGauss->SetParameter(1,xs);
//ss += fAsymmGauss->Eval(xx);
paras[0] = par[0]*par[4];
paras[1] = xs;
ss += AsymmGauss(&xx,paras);
}
xs += dx;
}
return ss + par[5];
}
//_____________________________________________________________________________
Double_t AliTRDCalibra::AsymmGauss(Double_t *x, Double_t *par)
{
//
// Function for the fit
//
//par[0] = normalization
//par[1] = mean
//par[2] = sigma
//norm0 = 1
//par[3] = lambda0
//par[4] = norm1
//par[5] = lambda1
Double_t par1save = par[1];
//Double_t par2save = par[2];
Double_t par2save = 0.0604244;
//Double_t par3save = par[3];
Double_t par3save = 7.65596;
//Double_t par5save = par[5];
Double_t par5save = 0.870597;
Double_t dx = x[0] - par1save;
Double_t sigma2 = par2save*par2save;
Double_t sqrt2 = TMath::Sqrt(2.0);
Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
* (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
* (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
//return par[0]*(exp1+par[4]*exp2);
return par[0] * (exp1 + 1.00124 * exp2);
}
//_____________________________________________________________________________
Double_t AliTRDCalibra::FuncLandauGaus(Double_t *x, Double_t *par)
{
//
// Sum Landau + Gaus with identical mean
//
Double_t valLandau = par[0] * TMath::Landau(x[0],par[1],par[2]);
//Double_t valGaus = par[3] * TMath::Gaus(x[0],par[4],par[5]);
Double_t valGaus = par[3] * TMath::Gaus(x[0],par[1],par[4]);
Double_t val = valLandau + valGaus;
return val;
}
//_____________________________________________________________________________
Double_t AliTRDCalibra::LanGauFun(Double_t *x, Double_t *par)
{
//
// Function for the fit
//
// Fit parameters:
// par[0]=Width (scale) parameter of Landau density
// par[1]=Most Probable (MP, location) parameter of Landau density
// par[2]=Total area (integral -inf to inf, normalization constant)
// par[3]=Width (sigma) of convoluted Gaussian function
//
// In the Landau distribution (represented by the CERNLIB approximation),
// the maximum is located at x=-0.22278298 with the location parameter=0.
// This shift is corrected within this function, so that the actual
// maximum is identical to the MP parameter.
//
// Numeric constants
Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
Double_t mpshift = -0.22278298; // Landau maximum location
// Control constants
Double_t np = 100.0; // Number of convolution steps
Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
// Variables
Double_t xx;
Double_t mpc;
Double_t fland;
Double_t sum = 0.0;
Double_t xlow;
Double_t xupp;
Double_t step;
Double_t i;
// MP shift correction
mpc = par[1] - mpshift * par[0];
// Range of convolution integral
xlow = x[0] - sc * par[3];
xupp = x[0] + sc * par[3];
step = (xupp - xlow) / np;
// Convolution integral of Landau and Gaussian by sum
for (i = 1.0; i <= np/2; i++) {
xx = xlow + (i-.5) * step;
fland = TMath::Landau(xx,mpc,par[0]) / par[0];
sum += fland * TMath::Gaus(x[0],xx,par[3]);
xx = xupp - (i-.5) * step;
fland = TMath::Landau(xx,mpc,par[0]) / par[0];
sum += fland * TMath::Gaus(x[0],xx,par[3]);
}
return (par[2] * step * sum * invsq2pi / par[3]);
}
//_____________________________________________________________________________
TF1 *AliTRDCalibra::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues
, Double_t *parlimitslo, Double_t *parlimitshi
, Double_t *fitparams, Double_t *fiterrors
, Double_t *chiSqr, Int_t *ndf)
{
//
// Function for the fit
//
Int_t i;
Char_t funname[100];
AliInfo(Form(funname,"Fitfcn_%s",his->GetName()));
TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
if (ffitold) {
delete ffitold;
}
TF1 *ffit = new TF1(funname,LanGauFun,fitrange[0],fitrange[1],4);
ffit->SetParameters(startvalues);
ffit->SetParNames("Width","MP","Area","GSigma");
for (i = 0; i < 4; i++) {
ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
}
his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
ffit->GetParameters(fitparams); // Obtain fit parameters
for (i = 0; i < 4; i++) {
fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
}
chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
ndf[0] = ffit->GetNDF(); // Obtain ndf
return (ffit); // Return fit function
}
//_____________________________________________________________________________
Int_t AliTRDCalibra::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm)
{
//
// Function for the fit
//
Double_t p;
Double_t x;
Double_t fy;
Double_t fxr;
Double_t fxl;
Double_t step;
Double_t l;
Double_t lold;
Int_t i = 0;
Int_t maxcalls = 10000;
// Search for maximum
p = params[1] - 0.1 * params[0];
step = 0.05 * params[0];
lold = -2.0;
l = -1.0;
while ((l != lold) && (i < maxcalls)) {
i++;
lold = l;
x = p + step;
l = LanGauFun(&x,params);
if (l < lold) {
step = -step / 10.0;
}
p += step;
}
if (i == maxcalls) {
return (-1);
}
maxx = x;
fy = l / 2.0;
// Search for right x location of fy
p = maxx + params[0];
step = params[0];
lold = -2.0;
l = -1e300;
i = 0;
while ( (l != lold) && (i < maxcalls) ) {
i++;
lold = l;
x = p + step;
l = TMath::Abs(LanGauFun(&x,params) - fy);
if (l > lold)
step = -step/10;
p += step;
}
if (i == maxcalls)
return (-2);
fxr = x;
// Search for left x location of fy
p = maxx - 0.5 * params[0];
step = -params[0];
lold = -2.0;
l = -1.0e300;
i = 0;
while ((l != lold) && (i < maxcalls)) {
i++;
lold = l;
x = p + step;
l = TMath::Abs(LanGauFun(&x,params) - fy);
if (l > lold) {
step = -step / 10.0;
}
p += step;
}
if (i == maxcalls) {
return (-3);
}
fxl = x;
fwhm = fxr - fxl;
return (0);
}