Converting PWG/TRD to native cmake
[u/mrichter/AliRoot.git] / MUON / MUONRecoCheck.C
CommitLineData
b8dc484b 1/**************************************************************************
48217459 2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
b8dc484b 15
e54bf126 16// $Id$
17
18/// \ingroup macros
19/// \file MUONRecoCheck.C
20/// \brief Utility macro to check the muon reconstruction.
21///
22/// Reconstructed tracks are compared to reference tracks. The reference tracks
23/// are built from AliTrackReference for the hit in chamber (0..9) and from
24/// kinematics (TreeK) for the vertex parameters.
25///
c687aa97 26/// \author Jean-Pierre Cussonneau, Philippe Pillot, Subatech
e54bf126 27
b8dc484b 28// ROOT includes
e1fe81b8 29#include <Riostream.h>
f202486b 30#include "TMath.h"
ce8cd162 31#include "TObjArray.h"
b8dc484b 32#include "TH1.h"
f202486b 33#include "TH2.h"
c687aa97 34#include "TH3.h"
f202486b 35#include "TGraphErrors.h"
e1fe81b8 36#include "TGraphAsymmErrors.h"
f202486b 37#include "TF1.h"
b8dc484b 38#include "TFile.h"
c687aa97 39#include "TCanvas.h"
40#include "TLegend.h"
e1fe81b8 41#include "TGeoManager.h"
b8dc484b 42
43// STEER includes
a99c3449 44#include "AliCDBManager.h"
e1fe81b8 45#include "AliGeomManager.h"
c687aa97 46#include "AliLog.h"
b8dc484b 47
48// MUON includes
a99c3449 49#include "AliMUONCDB.h"
b8dc484b 50#include "AliMUONConstants.h"
51#include "AliMUONTrack.h"
52#include "AliMUONRecoCheck.h"
53#include "AliMUONTrackParam.h"
a99c3449 54#include "AliMUONRecoParam.h"
48217459 55#include "AliMUONVTrackStore.h"
f202486b 56#include "AliMUONVCluster.h"
c687aa97 57#include "AliMUONTrackExtrap.h"
e1fe81b8 58#include "AliMUONESDInterface.h"
9f164762 59#include "AliMUONVTriggerTrackStore.h"
60#include "AliMUONTriggerTrack.h"
dd6a48d2 61#include "AliMpDEIterator.h"
b8dc484b 62
c687aa97 63Double_t langaufun(Double_t *x, Double_t *par);
e1fe81b8 64void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
c231afb9 65void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma);
dd6a48d2 66void FillResidual(TH1* h, Int_t i, Double_t& sigma, TGraphErrors* gMean, TGraphErrors* gSigma, Bool_t correctForSystematics, Bool_t fitResiduals);
e1fe81b8 67TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2);
68TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3);
69TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2 = 0x0, const char* fitting = "");
dd6a48d2 70void Zoom(TH1* h, Double_t fractionCut = 0.01);
c687aa97 71
72//------------------------------------------------------------------------------------
c231afb9 73void MUONRecoCheck(Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
74 const char* ocdbPath = "local://$ALICE_ROOT/OCDB", const Double_t pMin = 0., const Double_t pMax = 300.,
dd6a48d2 75 const Int_t pNBins = 30, Int_t absorberRegion = -1, Bool_t fitResiduals = kTRUE)
e54bf126 76{
e1fe81b8 77 /// Associate the reconstructed tracks with the simulated ones and check the quality of the reconstruction
78 /// (tracking/trigger efficiency; momentum, slope,... resolutions at first cluster and at vertex; cluster resolution).
c231afb9 79 /// - You can choose the momentum range and number of bins used to study the track resolution versus momentum.
80 /// - You can limit the calculation of track resolution at vertex to the tracks crossing the absorber in a given region
e1fe81b8 81 /// with the flag "absorberRegion": -1=all, 1=[2,3]deg, 2=[3,10]deg.
82
83 Double_t aAbsLimits[2];
84 if (absorberRegion > -1) {
85 if (absorberRegion == 1) {
86 aAbsLimits[0] = 2.;
87 aAbsLimits[1] = 3.;
88 } else if (absorberRegion == 2) {
89 aAbsLimits[0] = 3.;
90 aAbsLimits[1] = 10.;
91 } else {
92 cout<<"Unknown absorber region. Valid choices are: -1=all, 1=[2,3]deg, 2=[3,10]deg"<<endl;
93 return;
94 }
95 } else {
96 aAbsLimits[0] = 0.;
97 aAbsLimits[1] = 90.;
98 }
b8dc484b 99
c231afb9 100 cout<<endl<<"Momentum range for track resolution studies: "<<pNBins<<" bins in ["<<pMin<<","<<pMax<<"] GeV/c"<<endl;
101 if (pMin >= pMax || pNBins <= 0) {
102 cout<<"--> incorrect momentum range"<<endl<<endl;
103 return;
104 } else cout<<endl;
105
c687aa97 106 AliLog::SetClassDebugLevel("AliMCEvent",-1);
107
dd6a48d2 108 // ###################################### initialize ###################################### //
109 AliMUONRecoCheck rc(esdFileName, pathSim);
110
111 // load necessary data from OCDB
112 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
113 AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data", Form("local://%s",gSystem->pwd()));
114 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
115 if (!AliMUONCDB::LoadField()) return;
116 AliMUONTrackExtrap::SetField();
117 if (!AliMUONCDB::LoadMapping()) return;
118// AliGeomManager::LoadGeometry();
119 AliGeomManager::LoadGeometry("geometry.root");
120 if (!AliGeomManager::GetGeometry()) return;
121 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
122 if (!recoParam) return;
123 AliMUONESDInterface::ResetTracker(recoParam);
124
125 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
126 Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
127 // compute the mask of requested stations from recoParam
128 UInt_t requestedStationMask = 0;
129 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
130 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
131 Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
132
c687aa97 133 // ###################################### define histograms ###################################### //
b8dc484b 134 // File for histograms and histogram booking
135 TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
48217459 136
9f164762 137 TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks / evt",15,-0.5,14.5);
b8dc484b 138 TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
96ebe67e 139 TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
b8dc484b 140 TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
9f164762 141 TH1F *hTriggerable = new TH1F("hTriggerable"," Nb of triggerable tracks / evt",15,-0.5,14.5);
142 TH1F *hTriggered = new TH1F("hTriggered"," Nb of triggered tracks / evt",15,-0.5,14.5);
b8dc484b 143
c687aa97 144 // momentum resolution at vertex
145 histoFile->mkdir("momentumAtVertex","momentumAtVertex");
146 histoFile->cd("momentumAtVertex");
147
c231afb9 148 const Double_t pEdges[2] = {pMin, pMax};
c687aa97 149 const Int_t deltaPAtVtxNBins = 250;
c231afb9 150 Double_t deltaPAtVtxEdges[2];
151 if (pMax < 50.) { deltaPAtVtxEdges[0] = -20.; deltaPAtVtxEdges[1] = 5.; }
dd6a48d2 152 else { deltaPAtVtxEdges[0] = -50.; deltaPAtVtxEdges[1] = 50.; }
c687aa97 153
154 TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
155
156 TH2D *hResMomVertexVsMom = new TH2D("hResMomVertexVsMom","#Delta_{p} at vertex versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
157 TH2D *hResMomVertexVsMom_2_3_Deg = new TH2D("hResMomVertexVsMom_2_3_Deg","#Delta_{p} at vertex versus p for tracks between 2 and 3 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
158 TH2D *hResMomVertexVsMom_3_10_Deg = new TH2D("hResMomVertexVsMom_3_10_Deg","#Delta_{p} at vertex versus p for tracks between 3 and 10 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
159 TH2D *hResMomVertexVsMom_0_2_DegMC = new TH2D("hResMomVertexVsMom_0_2_DegMC","#Delta_{p} at vertex versus p for tracks with MC angle below 2 degrees;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtVtxNBins/10,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
160
e1fe81b8 161 TH2D *hResMomVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_0_2_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
162 TH2D *hResMomVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_2_3_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
163 TH2D *hResMomVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResMomVertexVsPosAbsEnd_3_10_DegMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
c687aa97 164
165 TH2D *hResMomVertexVsAngle = new TH2D("hResMomVertexVsAngle","#Delta_{p} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
166 TH2D *hResMomVertexVsMCAngle = new TH2D("hResMomVertexVsMCAngle","#Delta_{p} at vertex versus MC angle;MC angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
167 TH3D *hResMomVertexVsAngleVsMom = new TH3D("hResMomVertexVsAngleVsMom","#Delta_{p} at vertex versus track position at absorber end converted to degrees versus momentum;p (GeV/c);angle (Deg);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],100,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
168
e1fe81b8 169 TGraphAsymmErrors* gMeanResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
c687aa97 170 gMeanResMomVertexVsMom->SetName("gMeanResMomVertexVsMom");
171 gMeanResMomVertexVsMom->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
e1fe81b8 172 TGraphAsymmErrors* gMostProbResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
c687aa97 173 gMostProbResMomVertexVsMom->SetName("gMostProbResMomVertexVsMom");
174 gMostProbResMomVertexVsMom->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
e1fe81b8 175 TGraphAsymmErrors* gSigmaResMomVertexVsMom = new TGraphAsymmErrors(pNBins);
c687aa97 176 gSigmaResMomVertexVsMom->SetName("gSigmaResMomVertexVsMom");
e1fe81b8 177 gSigmaResMomVertexVsMom->SetTitle("#sigma_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
c687aa97 178
dd6a48d2 179 // transverse momentum resolution at vertex
180 TH2D *hResPtVertexVsPt = new TH2D("hResPtVertexVsPt","#Delta_{p_{t}} at vertex versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*pNBins,pEdges[0]/10.,pEdges[1]/10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0]/10.,deltaPAtVtxEdges[1]/10.);
181
c687aa97 182 // momentum resolution at first cluster
183 histoFile->mkdir("momentumAtFirstCluster","momentumAtFirstCluster");
184 histoFile->cd("momentumAtFirstCluster");
185
186 const Int_t deltaPAtFirstClNBins = 500;
c231afb9 187 Double_t deltaPAtFirstClEdges[2];
188 if (pMax < 25.) { deltaPAtFirstClEdges[0] = -5.; deltaPAtFirstClEdges[1] = 5.; }
189 else if (pMax < 50.) { deltaPAtFirstClEdges[0] = -10.; deltaPAtFirstClEdges[1] = 10.; }
190 else { deltaPAtFirstClEdges[0] = -25.; deltaPAtFirstClEdges[1] = 25.; }
c687aa97 191
192 TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
193 TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*pNBins,pEdges[0],pEdges[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
194
e1fe81b8 195 TGraphAsymmErrors* gMeanResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
c687aa97 196 gMeanResMomFirstClusterVsMom->SetName("gMeanResMomFirstClusterVsMom");
197 gMeanResMomFirstClusterVsMom->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
e1fe81b8 198 TGraphAsymmErrors* gSigmaResMomFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
c687aa97 199 gSigmaResMomFirstClusterVsMom->SetName("gSigmaResMomFirstClusterVsMom");
e1fe81b8 200 gSigmaResMomFirstClusterVsMom->SetTitle("#sigma_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
201
dd6a48d2 202 // transverse momentum resolution at vertex
203 TH2D *hResPtFirstClusterVsPt = new TH2D("hResPtFirstClusterVsPt","#Delta_{p_{t}} at first cluster versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*pNBins,pEdges[0]/10.,pEdges[1]/10.,deltaPAtFirstClNBins,deltaPAtFirstClEdges[0]/10.,deltaPAtFirstClEdges[1]/10.);
204
e1fe81b8 205 // angular resolution at vertex
206 histoFile->mkdir("slopesAtVertex","slopesAtVertex");
207 histoFile->cd("slopesAtVertex");
208
209 const Int_t deltaSlopeAtVtxNBins = 500;
210 const Double_t deltaSlopeAtVtxEdges[2] = {-0.05, 0.05};
211
212 TH1F *hResSlopeXVertex = new TH1F("hResSlopeXVertex","#Delta_{slope_{X}} at vertex;#Delta_{slope_{X}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
213 TH1F *hResSlopeYVertex = new TH1F("hResSlopeYVertex","#Delta_{slope_{Y}} at vertex;#Delta_{slope_{Y}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
214 TH2D *hResSlopeXVertexVsMom = new TH2D("hResSlopeXVertexVsMom","#Delta_{slope_{X}} at vertex versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
215 TH2D *hResSlopeYVertexVsMom = new TH2D("hResSlopeYVertexVsMom","#Delta_{slope_{Y}} at vertex versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
216
217 TH2D *hResSlopeXVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
218 TH2D *hResSlopeYVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_0_2_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
219 TH2D *hResSlopeXVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
220 TH2D *hResSlopeYVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_2_3_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
221 TH2D *hResSlopeXVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeXVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
222 TH2D *hResSlopeYVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResSlopeYVertexVsPosAbsEnd_3_10_DegMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
223
224 TH2D *hResSlopeXVertexVsAngle = new TH2D("hResSlopeXVertexVsAngle","#Delta_{slope_{X}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
225 TH2D *hResSlopeYVertexVsAngle = new TH2D("hResSlopeYVertexVsAngle","#Delta_{slope_{Y}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
226 TH2D *hResSlopeXVertexVsMCAngle = new TH2D("hResSlopeXVertexVsMCAngle","#Delta_{slope_{X}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
227 TH2D *hResSlopeYVertexVsMCAngle = new TH2D("hResSlopeYVertexVsMCAngle","#Delta_{slope_{Y}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
228
229 TGraphAsymmErrors* gMeanResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
230 gMeanResSlopeXVertexVsMom->SetName("gMeanResSlopeXVertexVsMom");
231 gMeanResSlopeXVertexVsMom->SetTitle("<#Delta_{slope_{X}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{X}}>");
232 TGraphAsymmErrors* gSigmaResSlopeXVertexVsMom = new TGraphAsymmErrors(pNBins);
233 gSigmaResSlopeXVertexVsMom->SetName("gSigmaResSlopeXVertexVsMom");
234 gSigmaResSlopeXVertexVsMom->SetTitle("#sigma_{slope_{X}} at vertex versus p;p (GeV/c);#sigma_{slope_{X}}");
235 TGraphAsymmErrors* gMeanResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
236 gMeanResSlopeYVertexVsMom->SetName("gMeanResSlopeYVertexVsMom");
237 gMeanResSlopeYVertexVsMom->SetTitle("<#Delta_{slope_{Y}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
238 TGraphAsymmErrors* gSigmaResSlopeYVertexVsMom = new TGraphAsymmErrors(pNBins);
239 gSigmaResSlopeYVertexVsMom->SetName("gSigmaResSlopeYVertexVsMom");
240 gSigmaResSlopeYVertexVsMom->SetTitle("#sigma_{slope_{Y}} at vertex versus p;p (GeV/c);#sigma_{slope_{Y}}");
241
242 // angular resolution at first cluster
243 histoFile->mkdir("slopesAtFirstCluster","slopesAtFirstCluster");
244 histoFile->cd("slopesAtFirstCluster");
245
246 const Int_t deltaSlopeAtFirstClNBins = 500;
247 const Double_t deltaSlopeAtFirstClEdges[2] = {-0.01, 0.01};
248
249 TH1F *hResSlopeXFirstCluster = new TH1F("hResSlopeXFirstCluster","#Delta_{slope_{X}} at first cluster;#Delta_{slope_{X}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
250 TH2D *hResSlopeXFirstClusterVsMom = new TH2D("hResSlopeXFirstClusterVsMom","#Delta_{slope_{X}} at first cluster versus p;p (GeV/c);#Delta_{slope_{X}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
251 TH1F *hResSlopeYFirstCluster = new TH1F("hResSlopeYFirstCluster","#Delta_{slope_{Y}} at first cluster;#Delta_{slope_{Y}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
252 TH2D *hResSlopeYFirstClusterVsMom = new TH2D("hResSlopeYFirstClusterVsMom","#Delta_{slope_{Y}} at first cluster versus p;p (GeV/c);#Delta_{slope_{Y}}",2*pNBins,pEdges[0],pEdges[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
253
254 TGraphAsymmErrors* gMeanResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
255 gMeanResSlopeXFirstClusterVsMom->SetName("gMeanResSlopeXFirstClusterVsMom");
256 gMeanResSlopeXFirstClusterVsMom->SetTitle("<#Delta_{slope_{X}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{X}}>");
257 TGraphAsymmErrors* gSigmaResSlopeXFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
258 gSigmaResSlopeXFirstClusterVsMom->SetName("gSigmaResSlopeXFirstClusterVsMom");
259 gSigmaResSlopeXFirstClusterVsMom->SetTitle("#sigma_{slope_{X}} at first cluster versus p;p (GeV/c);#sigma_{slope_{X}}");
260 TGraphAsymmErrors* gMeanResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
261 gMeanResSlopeYFirstClusterVsMom->SetName("gMeanResSlopeYFirstClusterVsMom");
262 gMeanResSlopeYFirstClusterVsMom->SetTitle("<#Delta_{slope_{Y}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
263 TGraphAsymmErrors* gSigmaResSlopeYFirstClusterVsMom = new TGraphAsymmErrors(pNBins);
264 gSigmaResSlopeYFirstClusterVsMom->SetName("gSigmaResSlopeYFirstClusterVsMom");
265 gSigmaResSlopeYFirstClusterVsMom->SetTitle("#sigma_{slope_{Y}} at first cluster versus p;p (GeV/c);#sigma_{slope_{Y}}");
266
267 // DCA resolution and MCS angular dispersion
268 histoFile->mkdir("DCA","DCA");
269 histoFile->cd("DCA");
270
271 const Int_t deltaPDCANBins = 500;
272 const Double_t deltaPDCAEdges[2] = {0., 1000.};
c231afb9 273 const Double_t deltaPMCSAngEdges[2] = {-1.5, 1.5};
e1fe81b8 274
275 TH1F *hPDCA = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
276 TH2D *hPDCAVsMom_2_3_Deg = new TH2D("hPDCAVsMom_2_3_Deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
277 TH2D *hPDCAVsMom_3_10_Deg = new TH2D("hPDCAVsMom_3_10_Deg","p #times DCA versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
278 TH2D *hPMCSAngVsMom_2_3_Deg = new TH2D("hPMCSAngVsMom_2_3_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
279 TH2D *hPMCSAngVsMom_3_10_Deg = new TH2D("hPMCSAngVsMom_3_10_Deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*pNBins,pEdges[0],pEdges[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
280
281 TH2D *hPDCAVsPosAbsEnd_0_2_DegMC = new TH2D("hPDCAVsPosAbsEnd_0_2_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
282 TH2D *hPDCAVsPosAbsEnd_2_3_DegMC = new TH2D("hPDCAVsPosAbsEnd_2_3_DegMC","p #times DCA}versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
283 TH2D *hPDCAVsPosAbsEnd_3_10_DegMC = new TH2D("hPDCAVsPosAbsEnd_3_10_DegMC","p #times DCA versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
284
285 TH2D *hPDCAVsAngle = new TH2D("hPDCAVsAngle","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
286 TH2D *hPDCAVsMCAngle = new TH2D("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
287
c231afb9 288 TGraphAsymmErrors* gMeanPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
289 gMeanPDCAVsMom_2_3_Deg->SetName("gMeanPDCAVsMom_2_3_Deg");
290 gMeanPDCAVsMom_2_3_Deg->SetTitle("<p #times DCA> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
e1fe81b8 291 TGraphAsymmErrors* gSigmaPDCAVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
292 gSigmaPDCAVsMom_2_3_Deg->SetName("gSigmaPDCAVsMom_2_3_Deg");
293 gSigmaPDCAVsMom_2_3_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
c231afb9 294 TGraphAsymmErrors* gMeanPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
295 gMeanPDCAVsMom_3_10_Deg->SetName("gMeanPDCAVsMom_3_10_Deg");
296 gMeanPDCAVsMom_3_10_Deg->SetTitle("<p #times DCA> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
e1fe81b8 297 TGraphAsymmErrors* gSigmaPDCAVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
298 gSigmaPDCAVsMom_3_10_Deg->SetName("gSigmaPDCAVsMom_3_10_Deg");
299 gSigmaPDCAVsMom_3_10_Deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
300 TGraphAsymmErrors* gMeanPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
301 gMeanPMCSAngVsMom_2_3_Deg->SetName("gMeanPMCSAngVsMom_2_3_Deg");
302 gMeanPMCSAngVsMom_2_3_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
303 TGraphAsymmErrors* gSigmaPMCSAngVsMom_2_3_Deg = new TGraphAsymmErrors(pNBins);
304 gSigmaPMCSAngVsMom_2_3_Deg->SetName("gSigmaPMCSAngVsMom_2_3_Deg");
305 gSigmaPMCSAngVsMom_2_3_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
306 TGraphAsymmErrors* gMeanPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
307 gMeanPMCSAngVsMom_3_10_Deg->SetName("gMeanPMCSAngVsMom_3_10_Deg");
308 gMeanPMCSAngVsMom_3_10_Deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
309 TGraphAsymmErrors* gSigmaPMCSAngVsMom_3_10_Deg = new TGraphAsymmErrors(pNBins);
310 gSigmaPMCSAngVsMom_3_10_Deg->SetName("gSigmaPMCSAngVsMom_3_10_Deg");
311 gSigmaPMCSAngVsMom_3_10_Deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
312
313 // eta resolution at vertex
314 histoFile->mkdir("etaAtVertex","etaAtVertex");
315 histoFile->cd("etaAtVertex");
316
317 const Int_t deltaEtaAtVtxNBins = 500;
318 const Double_t deltaEtaAtVtxEdges[2] = {-0.5, 0.5};
319
320 TH1F *hResEtaVertex = new TH1F("hResEtaVertex","#Delta_{eta} at vertex;#Delta_{eta}", deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
321 TH2D *hResEtaVertexVsMom = new TH2D("hResEtaVertexVsMom","#Delta_{eta} at vertex versus p;p (GeV/c);#Delta_{eta}",2*pNBins,pEdges[0],pEdges[1], deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
322
323 TH2D *hResEtaVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_0_2_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
324 TH2D *hResEtaVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_2_3_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
325 TH2D *hResEtaVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResEtaVertexVsPosAbsEnd_3_10_DegMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
326
327 TH2D *hResEtaVertexVsAngle = new TH2D("hResEtaVertexVsAngle","#Delta_{eta} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
328 TH2D *hResEtaVertexVsMCAngle = new TH2D("hResEtaVertexVsMCAngle","#Delta_{eta} at vertex versus MC angle;MC angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
329
330 TGraphAsymmErrors* gMeanResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
331 gMeanResEtaVertexVsMom->SetName("gMeanResEtaVertexVsMom");
332 gMeanResEtaVertexVsMom->SetTitle("<#Delta_{eta}> at vertex versus p;p (GeV/c);<#Delta_{eta}>");
333 TGraphAsymmErrors* gSigmaResEtaVertexVsMom = new TGraphAsymmErrors(pNBins);
334 gSigmaResEtaVertexVsMom->SetName("gSigmaResEtaVertexVsMom");
335 gSigmaResEtaVertexVsMom->SetTitle("#sigma_{eta} at vertex versus p;p (GeV/c);#sigma_{eta}");
336
337 // phi resolution at vertex
338 histoFile->mkdir("phiAtVertex","phiAtVertex");
339 histoFile->cd("phiAtVertex");
340
341 const Int_t deltaPhiAtVtxNBins = 500;
342 const Double_t deltaPhiAtVtxEdges[2] = {-0.5, 0.5};
343
344 TH1F *hResPhiVertex = new TH1F("hResPhiVertex","#Delta_{phi} at vertex;#Delta_{phi}", deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
345 TH2D *hResPhiVertexVsMom = new TH2D("hResPhiVertexVsMom","#Delta_{phi} at vertex versus p;p (GeV/c);#Delta_{phi}",2*pNBins,pEdges[0],pEdges[1], deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
346
347 TH2D *hResPhiVertexVsPosAbsEnd_0_2_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_0_2_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
348 TH2D *hResPhiVertexVsPosAbsEnd_2_3_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_2_3_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
349 TH2D *hResPhiVertexVsPosAbsEnd_3_10_DegMC = new TH2D("hResPhiVertexVsPosAbsEnd_3_10_DegMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
350
351 TH2D *hResPhiVertexVsAngle = new TH2D("hResPhiVertexVsAngle","#Delta_{phi} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
352 TH2D *hResPhiVertexVsMCAngle = new TH2D("hResPhiVertexVsMCAngle","#Delta_{phi} at vertex versus MC angle;MC angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
353
354 TGraphAsymmErrors* gMeanResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
355 gMeanResPhiVertexVsMom->SetName("gMeanResPhiVertexVsMom");
356 gMeanResPhiVertexVsMom->SetTitle("<#Delta_{phi}> at vertex versus p;p (GeV/c);<#Delta_{phi}>");
357 TGraphAsymmErrors* gSigmaResPhiVertexVsMom = new TGraphAsymmErrors(pNBins);
358 gSigmaResPhiVertexVsMom->SetName("gSigmaResPhiVertexVsMom");
359 gSigmaResPhiVertexVsMom->SetTitle("#sigma_{phi} at vertex versus p;p (GeV/c);#sigma_{phi}");
c687aa97 360
361 // cluster resolution
362 histoFile->mkdir("clusters","clusters");
363 histoFile->cd("clusters");
364
dd6a48d2 365 // find the highest chamber resolution and set histogram bins
366 const Int_t clusterResNBins = 5000;
367 Double_t maxRes[2] = {-1., -1.};
368 for (Int_t i = 0; i < 10; i++) {
369 if (recoParam->GetDefaultNonBendingReso(i) > maxRes[0]) maxRes[0] = recoParam->GetDefaultNonBendingReso(i);
370 if (recoParam->GetDefaultBendingReso(i) > maxRes[1]) maxRes[1] = recoParam->GetDefaultBendingReso(i);
371 }
372 Double_t clusterResMaxX = 10.*maxRes[0];
373 Double_t clusterResMaxY = 10.*maxRes[1];
374
f202486b 375 TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
376 TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
377 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
dd6a48d2 378 hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), clusterResNBins, -clusterResMaxX, clusterResMaxX);
379 hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), clusterResNBins, -clusterResMaxY, clusterResMaxY);
380 }
381
382 // fill correspondence between DEId and indices for histo (starting from 1)
383 Int_t deNBins = 0;
384 Int_t deIndices[1100];
385 Int_t deIds[200];
386 AliMpDEIterator it;
387 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
388 it.First(i);
389 while (!it.IsDone()) {
390 deNBins++;
391 deIndices[it.CurrentDEId()] = deNBins;
392 deIds[deNBins] = it.CurrentDEId();
393 it.Next();
394 }
f202486b 395 }
dd6a48d2 396 TH2F* hResidualXPerDE = new TH2F("hResidualXPerDE", "cluster-track residual-X distribution per DE;DE ID;#Delta_{X} (cm)", deNBins, 0.5, deNBins+0.5, clusterResNBins, -clusterResMaxX, clusterResMaxX);
397 for (Int_t i = 1; i <= deNBins; i++) hResidualXPerDE->GetXaxis()->SetBinLabel(i, Form("%d",deIds[i]));
398 TH2F* hResidualYPerDE = new TH2F("hResidualYPerDE", "cluster-track residual-Y distribution per DE;DE ID;#Delta_{Y} (cm)", deNBins, 0.5, deNBins+0.5, clusterResNBins, -clusterResMaxY, clusterResMaxY);
399 for (Int_t i = 1; i <= deNBins; i++) hResidualYPerDE->GetXaxis()->SetBinLabel(i, Form("%d",deIds[i]));
c687aa97 400
f202486b 401 TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
402 gResidualXPerChMean->SetName("gResidualXPerChMean");
e1fe81b8 403 gResidualXPerChMean->SetTitle("cluster-trackRef residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
f202486b 404 gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
405 TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
406 gResidualYPerChMean->SetName("gResidualYPerChMean");
e1fe81b8 407 gResidualYPerChMean->SetTitle("cluster-trackRef residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
f202486b 408 gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
409 TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
410 gResidualXPerChSigma->SetName("gResidualXPerChSigma");
e1fe81b8 411 gResidualXPerChSigma->SetTitle("cluster-trackRef residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
f202486b 412 gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
413 TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
414 gResidualYPerChSigma->SetName("gResidualYPerChSigma");
e1fe81b8 415 gResidualYPerChSigma->SetTitle("cluster-trackRef residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
f202486b 416 gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
dd6a48d2 417
418 TGraphErrors* gResidualXPerDEMean = new TGraphErrors(deNBins);
419 gResidualXPerDEMean->SetName("gResidualXPerDEMean");
420 gResidualXPerDEMean->SetTitle("cluster-trackRef residual-X per DE: mean;DE ID;<#Delta_{X}> (cm)");
421 gResidualXPerDEMean->SetMarkerStyle(kFullDotLarge);
422 TGraphErrors* gResidualYPerDEMean = new TGraphErrors(deNBins);
423 gResidualYPerDEMean->SetName("gResidualYPerDEMean");
424 gResidualYPerDEMean->SetTitle("cluster-trackRef residual-Y per dE: mean;DE ID;<#Delta_{Y}> (cm)");
425 gResidualYPerDEMean->SetMarkerStyle(kFullDotLarge);
426 TGraphErrors* gResidualXPerDESigma = new TGraphErrors(deNBins);
427 gResidualXPerDESigma->SetName("gResidualXPerDESigma");
428 gResidualXPerDESigma->SetTitle("cluster-trackRef residual-X per DE: sigma;DE ID;#sigma_{X} (cm)");
429 gResidualXPerDESigma->SetMarkerStyle(kFullDotLarge);
430 TGraphErrors* gResidualYPerDESigma = new TGraphErrors(deNBins);
431 gResidualYPerDESigma->SetName("gResidualYPerDESigma");
432 gResidualYPerDESigma->SetTitle("cluster-trackRef residual-Y per DE: sigma;DE ID;#sigma_{Y} (cm)");
433 gResidualYPerDESigma->SetMarkerStyle(kFullDotLarge);
434
9f164762 435 histoFile->mkdir("trigger");
436 histoFile->cd("trigger");
437 TH1F* hResidualTrigX11 = new TH1F("hResiudalTrigX11", "Residual X11", 100, -10., 10.);
438 TH1F* hResidualTrigY11 = new TH1F("hResiudalTrigY11", "Residual Y11", 100, -10., 10.);
439 TH1F* hResidualTrigSlopeY = new TH1F("hResiudalTrigSlopeY", "Residual Y slope", 100, -0.1, 0.1);
440 TH1F* hTriggerableMatchFailed = new TH1F("hTriggerableMatchFailed", "Triggerable multiplicity for events with no match", 15, -0.5, 14.5);
441
dd6a48d2 442 // ###################################### fill histograms ###################################### //
b8dc484b 443 Int_t ievent;
444 Int_t nReconstructibleTracks = 0;
445 Int_t nReconstructedTracks = 0;
446 Int_t nReconstructibleTracksCheck = 0;
f202486b 447 AliMUONTrackParam *trackParam;
e1fe81b8 448 Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
449 Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
248ee365 450 Double_t dPhi;
e1fe81b8 451 Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
452 Double_t xDCA,yDCA,dca,pU;
453 Double_t aMCSMoy = 0., aMCS2Moy = 0., dMCSMoy = 0., dMCS2Moy = 0., adMCSMoy = 0.;
454 Int_t nMCS = 0;
48217459 455
dd6a48d2 456 Int_t nevents = rc.NumberOfEvents();
457 if (nevents < nEvent || nEvent < 0) nEvent = nevents;
458
459 // loop over events
48217459 460 for (ievent=0; ievent<nEvent; ievent++)
461 {
c687aa97 462 if ((ievent+1)%100 == 0) cout<<"\rEvent processing... "<<ievent+1<<flush;
b8dc484b 463
a99c3449 464 AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
f202486b 465 AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
b8dc484b 466
48217459 467 hReconstructible->Fill(trackRefStore->GetSize());
468 hReco->Fill(trackStore->GetSize());
469
470 nReconstructibleTracks += trackRefStore->GetSize();
471 nReconstructedTracks += trackStore->GetSize();
9f164762 472
473 AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(ievent);
474 AliMUONVTriggerTrackStore* triggerTrackStore = rc.TriggeredTracks(ievent);
475
476 hTriggerable->Fill(triggerTrackRefStore->GetSize());
477 hTriggered->Fill(triggerTrackStore->GetSize());
478
479 // loop over trigger trackRef
480 TIter nextTrig(triggerTrackRefStore->CreateIterator());
481 AliMUONTriggerTrack* triggerTrackRef;
482 Int_t nTriggerMatches = 0;
483 while ( ( triggerTrackRef = static_cast<AliMUONTriggerTrack*>(nextTrig()) ) )
484 {
485
486 AliMUONTriggerTrack* triggerTrackMatched = 0x0;
487
488 // loop over trackReco and look for compatible track
489 TIter nextTrig2(triggerTrackStore->CreateIterator());
490 AliMUONTriggerTrack* triggerTrackReco;
491 while ( ( triggerTrackReco = static_cast<AliMUONTriggerTrack*>(nextTrig2()) ) )
492 {
493
494 // check if trackReco is compatible with trackRef
495 if (triggerTrackReco->Match(*triggerTrackRef, sigmaCut)) {
496 triggerTrackMatched = triggerTrackReco;
497 nTriggerMatches++;
498 break;
499 }
500 }
501
502 if (triggerTrackMatched) { // tracking requirements verified, track is found
503 hResidualTrigX11->Fill( triggerTrackMatched->GetX11() - triggerTrackRef->GetX11() );
504 hResidualTrigY11->Fill( triggerTrackMatched->GetY11() - triggerTrackRef->GetY11() );
505 hResidualTrigSlopeY->Fill( triggerTrackMatched->GetSlopeY() - triggerTrackRef->GetSlopeY() );
506 }
507 } // loop on trigger track ref
508
509 if ( nTriggerMatches != triggerTrackStore->GetSize() )
510 hTriggerableMatchFailed->Fill(triggerTrackRefStore->GetSize());
48217459 511
f202486b 512 // loop over trackRef
48217459 513 TIter next(trackRefStore->CreateIterator());
514 AliMUONTrack* trackRef;
48217459 515 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
516 {
48217459 517
9f164762 518 hTrackRefID->Fill(trackRef->GetUniqueID());
f202486b 519
520 AliMUONTrack* trackMatched = 0x0;
521 Int_t nMatchClusters = 0;
522
523 // loop over trackReco and look for compatible track
48217459 524 TIter next2(trackStore->CreateIterator());
525 AliMUONTrack* trackReco;
48217459 526 while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
527 {
f202486b 528
529 // check if trackReco is compatible with trackRef
f202486b 530 if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
531 trackMatched = trackReco;
f202486b 532 break;
61fed964 533 }
f202486b 534
48217459 535 }
536
f202486b 537 if (trackMatched) { // tracking requirements verified, track is found
48217459 538 nReconstructibleTracksCheck++;
f202486b 539 hNClusterComp->Fill(nMatchClusters);
c687aa97 540
541 // compute track position at the end of the absorber
542 AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)trackMatched->GetTrackParamAtCluster()->First()));
543 AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
544 xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
545 yAbs = trackParamAtAbsEnd.GetBendingCoor();
546 dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
547 aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
e1fe81b8 548 pX2 = trackParamAtAbsEnd.Px();
549 pY2 = trackParamAtAbsEnd.Py();
550 pZ2 = trackParamAtAbsEnd.Pz();
551 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
552 aMCS = TMath::ATan(-pT2/pZ2) * TMath::RadToDeg();
c687aa97 553
48217459 554 trackParam = trackRef->GetTrackParamAtVertex();
555 x1 = trackParam->GetNonBendingCoor();
556 y1 = trackParam->GetBendingCoor();
557 z1 = trackParam->GetZ();
e1fe81b8 558 slopex1 = trackParam->GetNonBendingSlope();
559 slopey1 = trackParam->GetBendingSlope();
48217459 560 pX1 = trackParam->Px();
561 pY1 = trackParam->Py();
562 pZ1 = trackParam->Pz();
563 p1 = trackParam->P();
f202486b 564 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
c687aa97 565 aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
e1fe81b8 566 eta1 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT1/pZ1)));
567 phi1 = TMath::Pi()+TMath::ATan2(-pY1, -pX1);
c687aa97 568
f202486b 569 trackParam = trackMatched->GetTrackParamAtVertex();
48217459 570 x2 = trackParam->GetNonBendingCoor();
571 y2 = trackParam->GetBendingCoor();
572 z2 = trackParam->GetZ();
e1fe81b8 573 slopex2 = trackParam->GetNonBendingSlope();
574 slopey2 = trackParam->GetBendingSlope();
48217459 575 pX2 = trackParam->Px();
576 pY2 = trackParam->Py();
577 pZ2 = trackParam->Pz();
578 p2 = trackParam->P();
f202486b 579 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
e1fe81b8 580 eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
581 phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
48217459 582
248ee365 583 dPhi = phi2-phi1;
584 if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
585 else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
586
e1fe81b8 587 AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First()));
588 pU = trackParamAtDCA.P();
589 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
590 xDCA = trackParamAtDCA.GetNonBendingCoor();
591 yDCA = trackParamAtDCA.GetBendingCoor();
592 dca = TMath::Sqrt(xDCA*xDCA + yDCA*yDCA);
593
48217459 594 hResMomVertex->Fill(p2-p1);
e1fe81b8 595 hResSlopeXVertex->Fill(slopex2-slopex1);
596 hResSlopeYVertex->Fill(slopey2-slopey1);
597 hPDCA->Fill(0.5*(p2+pU)*dca);
598 hResEtaVertex->Fill(eta2-eta1);
248ee365 599 hResPhiVertex->Fill(dPhi);
e1fe81b8 600 if (aMC >= aAbsLimits[0] && aMC <= aAbsLimits[1]) {
601 hResMomVertexVsMom->Fill(p1,p2-p1);
602 hResSlopeXVertexVsMom->Fill(p1,slopex2-slopex1);
603 hResSlopeYVertexVsMom->Fill(p1,slopey2-slopey1);
604 hResEtaVertexVsMom->Fill(p1,eta2-eta1);
248ee365 605 hResPhiVertexVsMom->Fill(p1,dPhi);
dd6a48d2 606 hResPtVertexVsPt->Fill(pT1,pT2-pT1);
e1fe81b8 607 }
c687aa97 608 hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
e1fe81b8 609 if (aAbs > 2. && aAbs < 3.) {
610 hResMomVertexVsMom_2_3_Deg->Fill(p1,p2-p1);
611 hPDCAVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*dca);
612 hPMCSAngVsMom_2_3_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
613 }
614 else if (aAbs >= 3. && aAbs < 10.) {
615 hResMomVertexVsMom_3_10_Deg->Fill(p1,p2-p1);
616 hPDCAVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*dca);
617 hPMCSAngVsMom_3_10_Deg->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
618 aMCSMoy += 0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad();
619 aMCS2Moy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
620 dMCSMoy += 0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd());
621 dMCS2Moy += (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd())) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
622 adMCSMoy += (0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad()) * (0.5*(p2+pU)*(dAbs-pT1/pZ1*AliMUONConstants::AbsZEnd()));
623 nMCS++;
624 }
c687aa97 625 if (aMC < 2.) {
626 hResMomVertexVsMom_0_2_DegMC->Fill(p1,p2-p1);
627 hResMomVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,p2-p1);
e1fe81b8 628 hResSlopeXVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopex2-slopex1);
629 hResSlopeYVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,slopey2-slopey1);
630 hPDCAVsPosAbsEnd_0_2_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
631 hResEtaVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,eta2-eta1);
248ee365 632 hResPhiVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,dPhi);
e1fe81b8 633 }
634 else if (aMC >= 2. && aMC < 3) {
635 hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
636 hResSlopeXVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopex2-slopex1);
637 hResSlopeYVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,slopey2-slopey1);
638 hPDCAVsPosAbsEnd_2_3_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
639 hResEtaVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,eta2-eta1);
248ee365 640 hResPhiVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,dPhi);
e1fe81b8 641 }
642 else if (aMC >= 3. && aMC < 10.) {
643 hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
644 hResSlopeXVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopex2-slopex1);
645 hResSlopeYVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,slopey2-slopey1);
646 hPDCAVsPosAbsEnd_3_10_DegMC->Fill(dAbs,0.5*(p2+pU)*dca);
647 hResEtaVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,eta2-eta1);
248ee365 648 hResPhiVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,dPhi);
c687aa97 649 }
c687aa97 650 hResMomVertexVsAngle->Fill(aAbs,p2-p1);
e1fe81b8 651 hResSlopeXVertexVsAngle->Fill(aAbs,slopex2-slopex1);
652 hResSlopeYVertexVsAngle->Fill(aAbs,slopey2-slopey1);
653 hPDCAVsAngle->Fill(aAbs,0.5*(p2+pU)*dca);
654 hResEtaVertexVsAngle->Fill(aAbs,eta2-eta1);
248ee365 655 hResPhiVertexVsAngle->Fill(aAbs,dPhi);
c687aa97 656 hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
e1fe81b8 657 hResSlopeXVertexVsMCAngle->Fill(aMC,slopex2-slopex1);
658 hResSlopeYVertexVsMCAngle->Fill(aMC,slopey2-slopey1);
659 hPDCAVsMCAngle->Fill(aMC,0.5*(p2+pU)*dca);
660 hResEtaVertexVsMCAngle->Fill(aMC,eta2-eta1);
248ee365 661 hResPhiVertexVsMCAngle->Fill(aMC,dPhi);
f202486b 662
96ebe67e 663 trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
48217459 664 x1 = trackParam->GetNonBendingCoor();
665 y1 = trackParam->GetBendingCoor();
666 z1 = trackParam->GetZ();
e1fe81b8 667 slopex1 = trackParam->GetNonBendingSlope();
668 slopey1 = trackParam->GetBendingSlope();
48217459 669 pX1 = trackParam->Px();
670 pY1 = trackParam->Py();
671 pZ1 = trackParam->Pz();
672 p1 = trackParam->P();
f202486b 673 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
674
f202486b 675 trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
48217459 676 x2 = trackParam->GetNonBendingCoor();
677 y2 = trackParam->GetBendingCoor();
678 z2 = trackParam->GetZ();
e1fe81b8 679 slopex2 = trackParam->GetNonBendingSlope();
680 slopey2 = trackParam->GetBendingSlope();
48217459 681 pX2 = trackParam->Px();
682 pY2 = trackParam->Py();
683 pZ2 = trackParam->Pz();
684 p2 = trackParam->P();
f202486b 685 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
48217459 686
96ebe67e 687 hResMomFirstCluster->Fill(p2-p1);
f202486b 688 hResMomFirstClusterVsMom->Fill(p1,p2-p1);
dd6a48d2 689 hResPtFirstClusterVsPt->Fill(pT1,pT2-pT1);
f202486b 690
e1fe81b8 691 hResSlopeXFirstCluster->Fill(slopex2-slopex1);
692 hResSlopeYFirstCluster->Fill(slopey2-slopey1);
693 hResSlopeXFirstClusterVsMom->Fill(p1,slopex2-slopex1);
694 hResSlopeYFirstClusterVsMom->Fill(p1,slopey2-slopey1);
695
f202486b 696 // Fill residuals
f202486b 697 AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
698 while (trackParamAtCluster1) {
699 AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
700 AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
701 while (trackParamAtCluster2) {
702 AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
703 if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
704 hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
705 hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
dd6a48d2 706 hResidualXPerDE->Fill(deIndices[cluster1->GetDetElemId()], cluster1->GetX() - cluster2->GetX());
707 hResidualYPerDE->Fill(deIndices[cluster1->GetDetElemId()], cluster1->GetY() - cluster2->GetY());
f202486b 708 break;
709 }
710 trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
711 }
712 trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
713 }
714
b8dc484b 715 }
f202486b 716
b8dc484b 717 } // end loop track ref.
b8dc484b 718
48217459 719 } // end loop on event
c687aa97 720 cout<<"\rEvent processing... "<<nevents<<" done"<<endl;
48217459 721
c687aa97 722 // ###################################### compute stuff ###################################### //
e1fe81b8 723 cout<<"\nWhen not specified, resolution at vertex is computed for ";
724 if (absorberRegion == 1) cout<<"tracks in the absorber region [2,3] deg."<<endl;
725 else if (absorberRegion == 2) cout<<"tracks in the absorber region [3,10] deg."<<endl;
726 else cout<<"all tracks"<<endl;
727
c687aa97 728 // compute momentum resolution at vertex versus p
e1fe81b8 729 TF1 *f2 = new TF1("f2",langaufun,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1],4);
c687aa97 730 Int_t rebinFactorX = TMath::Max(hResMomVertexVsMom->GetNbinsX()/pNBins, 1);
731 for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom->GetNbinsX(); i+=rebinFactorX) {
732 cout<<"\rFitting momentum residuals at vertex... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
733 TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
c687aa97 734 f2->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
735 tmp->Fit("f2","WWNQ");
736 Double_t fwhm = f2->GetParameter(0);
737 Double_t sigma = f2->GetParameter(3);
738 Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
739 Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/tmp->GetBinWidth(1)),1);
740 while (deltaPAtVtxNBins%rebin!=0) rebin--;
741 tmp->Rebin(rebin);
742 tmp->Fit("f2","NQ");
743 fwhm = f2->GetParameter(0);
744 sigma = f2->GetParameter(3);
745 sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
746 Double_t fwhmErr = f2->GetParError(0);
747 Double_t sigmaErr = f2->GetParError(3);
748 Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
e1fe81b8 749 hResMomVertexVsMom->GetXaxis()->SetRange(i-rebinFactorX+1,i);
248ee365 750 Double_t p = (tmp->GetEntries() > 0) ? hResMomVertexVsMom->GetMean() : 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
e1fe81b8 751 hResMomVertexVsMom->GetXaxis()->SetRange();
752 Double_t pErr[2] = {p-hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1), hResMomVertexVsMom->GetBinLowEdge(i+1)-p};
753 gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
754 gMeanResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], tmp->GetMeanError(), tmp->GetMeanError());
755 gMostProbResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, -f2->GetParameter(1));
756 gMostProbResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], f2->GetParError(1), f2->GetParError(1));
757 gSigmaResMomVertexVsMom->SetPoint(i/rebinFactorX-1, p, 100.*sigmaP/p);
758 gSigmaResMomVertexVsMom->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], 100.*sigmaPErr/p, 100.*sigmaPErr/p);
f202486b 759 delete tmp;
760 }
c687aa97 761 cout<<"\rFitting momentum residuals at vertex... "<<pNBins<<"/"<<pNBins<<endl;
762
e1fe81b8 763 // compute momentum relative resolution at first cluster versus p
764 FitGausResVsMom(hResMomFirstClusterVsMom, pNBins, 0., 1., "momentum residuals at first cluster", gMeanResMomFirstClusterVsMom, gSigmaResMomFirstClusterVsMom);
c687aa97 765 rebinFactorX = TMath::Max(hResMomFirstClusterVsMom->GetNbinsX()/pNBins, 1);
766 for (Int_t i = rebinFactorX; i <= hResMomFirstClusterVsMom->GetNbinsX(); i+=rebinFactorX) {
e1fe81b8 767 Double_t x,y;
768 gSigmaResMomFirstClusterVsMom->GetPoint(i/rebinFactorX-1, x, y);
769 gSigmaResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
770 gSigmaResMomFirstClusterVsMom->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYlow(i/rebinFactorX-1)/x);
771 gSigmaResMomFirstClusterVsMom->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResMomFirstClusterVsMom->GetErrorYhigh(i/rebinFactorX-1)/x);
f202486b 772 }
773
e1fe81b8 774 // compute slopeX resolution at vertex versus p
775 FitGausResVsMom(hResSlopeXVertexVsMom, pNBins, 0., 2.e-3, "slopeX residuals at vertex", gMeanResSlopeXVertexVsMom, gSigmaResSlopeXVertexVsMom);
776
777 // compute slopeY resolution at vertex versus p
778 FitGausResVsMom(hResSlopeYVertexVsMom, pNBins, 0., 2.e-3, "slopeY residuals at vertex", gMeanResSlopeYVertexVsMom, gSigmaResSlopeYVertexVsMom);
779
780 // compute slopeX resolution at first cluster versus p
781 FitGausResVsMom(hResSlopeXFirstClusterVsMom, pNBins, 0., 3.e-4, "slopeX residuals at first cluster", gMeanResSlopeXFirstClusterVsMom, gSigmaResSlopeXFirstClusterVsMom);
782
783 // compute slopeY resolution at first cluster versus p
784 FitGausResVsMom(hResSlopeYFirstClusterVsMom, pNBins, 0., 2.e-4, "slopeY residuals at first cluster", gMeanResSlopeYFirstClusterVsMom, gSigmaResSlopeYFirstClusterVsMom);
785
786 // compute p*DCA resolution in the region [2,3] deg at absorber end
c231afb9 787 FitPDCAVsMom(hPDCAVsMom_2_3_Deg, pNBins, "p*DCA (tracks in [2,3] deg.)", gMeanPDCAVsMom_2_3_Deg, gSigmaPDCAVsMom_2_3_Deg);
e1fe81b8 788
789 // compute p*DCA resolution in the region [3,10] deg at absorber end
c231afb9 790 FitPDCAVsMom(hPDCAVsMom_3_10_Deg, pNBins, "p*DCA (tracks in [3,10] deg.)", gMeanPDCAVsMom_3_10_Deg, gSigmaPDCAVsMom_3_10_Deg);
e1fe81b8 791
792 // compute MCS angular dispersion in the region [2,3] deg at absorber end
793 FitGausResVsMom(hPMCSAngVsMom_2_3_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsMom_2_3_Deg, gSigmaPMCSAngVsMom_2_3_Deg);
794
795 // compute MCS angular dispersion in the region [3,10] deg at absorber end
796 FitGausResVsMom(hPMCSAngVsMom_3_10_Deg, pNBins, 0., 2.e-3, "p*MCSAngle (tracks in [3,10] deg.)", gMeanPMCSAngVsMom_3_10_Deg, gSigmaPMCSAngVsMom_3_10_Deg);
797
798 // compute eta resolution at vertex versus p
799 FitGausResVsMom(hResEtaVertexVsMom, pNBins, 0., 0.1, "eta residuals at vertex", gMeanResEtaVertexVsMom, gSigmaResEtaVertexVsMom);
800
801 // compute phi resolution at vertex versus p
802 FitGausResVsMom(hResPhiVertexVsMom, pNBins, 0., 0.01, "phi residuals at vertex", gMeanResPhiVertexVsMom, gSigmaResPhiVertexVsMom);
803
dd6a48d2 804 // compute cluster-track residual mean and dispersion per chamber
805 Double_t clusterResPerCh[10][2];
f202486b 806 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
dd6a48d2 807 FillResidual(hResidualXInCh[i], i, clusterResPerCh[i][0], gResidualXPerChMean, gResidualXPerChSigma, kTRUE, fitResiduals);
808 FillResidual(hResidualYInCh[i], i, clusterResPerCh[i][1], gResidualYPerChMean, gResidualYPerChSigma, kTRUE, fitResiduals);
809 }
810
811 // compute cluster-track residual mean and dispersion per DE
812 Double_t clusterResPerDE[200][2];
813 for (Int_t i = 0; i < deNBins; i++) {
814 TH1D *tmp = hResidualXPerDE->ProjectionY("tmp",i+1,i+1,"e");
815 FillResidual(tmp, i, clusterResPerDE[i][0], gResidualXPerDEMean, gResidualXPerDESigma, kTRUE, fitResiduals);
816 delete tmp;
817 tmp = hResidualYPerDE->ProjectionY("tmp",i+1,i+1,"e");
818 FillResidual(tmp, i, clusterResPerDE[i][1], gResidualYPerDEMean, gResidualYPerDESigma, kTRUE, fitResiduals);
819 delete tmp;
f202486b 820 }
821
c687aa97 822 // ###################################### display histograms ###################################### //
e1fe81b8 823 // diplay momentum residuals
824 TCanvas* cResMom = DrawVsAng("cResMom", "momentum residual at vertex in 3 angular regions", hResMomVertex, hResMomVertexVsAngle);
825 TCanvas* cResMomMC = DrawVsAng("cResMomMC", "momentum residual at vertex in 3 MC angular regions", hResMomVertex, hResMomVertexVsMCAngle);
826 TCanvas* cResMomVsPos = DrawVsPos("cResMomVsPos", "momentum residual at vertex versus position at absorber end in 3 MC angular regions",
827 hResMomVertexVsPosAbsEnd_0_2_DegMC, hResMomVertexVsPosAbsEnd_2_3_DegMC, hResMomVertexVsPosAbsEnd_3_10_DegMC);
828 TCanvas* cResMom_2_3_Deg = DrawResMomVsMom("cResMom_2_3_Deg", "momentum residual for tracks between 2 and 3 degrees",
829 hResMomVertexVsMom_2_3_Deg, 10, f2, "momentum residuals at vertex (tracks in [2,3] deg.)");
830 TCanvas* cResMom_3_10_Deg = DrawResMomVsMom("cResMom_3_10_Deg", "momentum residual for tracks between 3 and 10 degrees",
831 hResMomVertexVsMom_3_10_Deg, 10, f2, "momentum residuals at vertex (tracks in [3,10] deg.)");
832 TCanvas* cResMom_0_2_DegMC = DrawResMomVsMom("cResMom_0_2_DegMC", "momentum residuals for tracks with MC angle < 2 degrees", hResMomVertexVsMom_0_2_DegMC, 5);
833
834 // diplay slopeX residuals
835 TCanvas* cResSlopeX = DrawVsAng("cResSlopeX", "slope_{X} residual at vertex in 3 angular regions", hResSlopeXVertex, hResSlopeXVertexVsAngle);
836 TCanvas* cResSlopeXMC = DrawVsAng("cResSlopeXMC", "slope_{X} residual at vertex in 3 MC angular regions", hResSlopeXVertex, hResSlopeXVertexVsMCAngle);
837 TCanvas* cResSlopeXVsPos = DrawVsPos("cResSlopeXVsPos", "slope_{X} residual at vertex versus position at absorber end in 3 MC angular regions",
838 hResSlopeXVertexVsPosAbsEnd_0_2_DegMC, hResSlopeXVertexVsPosAbsEnd_2_3_DegMC, hResSlopeXVertexVsPosAbsEnd_3_10_DegMC);
839
840 // diplay slopeY residuals
841 TCanvas* cResSlopeY = DrawVsAng("cResSlopeY", "slope_{Y} residual at vertex in 3 angular regions", hResSlopeYVertex, hResSlopeYVertexVsAngle);
842 TCanvas* cResSlopeYMC = DrawVsAng("cResSlopeYMC", "slope_{Y} residual at vertex in 3 MC angular regions", hResSlopeYVertex, hResSlopeYVertexVsMCAngle);
843 TCanvas* cResSlopeYVsPos = DrawVsPos("cResSlopeYVsPos", "slope_{Y} residual at vertex versus position at absorber end in 3 MC angular regions",
844 hResSlopeYVertexVsPosAbsEnd_0_2_DegMC, hResSlopeYVertexVsPosAbsEnd_2_3_DegMC, hResSlopeYVertexVsPosAbsEnd_3_10_DegMC);
845
846 // diplay P*DCA
847 TCanvas* cPDCA = DrawVsAng("cPDCA", "p #times DCA in 3 angular regions", hPDCA, hPDCAVsAngle);
848 TCanvas* cPDCAMC = DrawVsAng("cPDCAMC", "p #times DCA in 3 MC angular regions", hPDCA, hPDCAVsMCAngle);
849 TCanvas* cPDCAVsPos = DrawVsPos("cPDCAVsPos", "p #times DCA versus position at absorber end in 3 MC angular regions",
850 hPDCAVsPosAbsEnd_0_2_DegMC, hPDCAVsPosAbsEnd_2_3_DegMC, hPDCAVsPosAbsEnd_3_10_DegMC);
851
852 // diplay eta residuals
853 TCanvas* cResEta = DrawVsAng("cResEta", "eta residual at vertex in 3 angular regions", hResEtaVertex, hResEtaVertexVsAngle);
854 TCanvas* cResEtaMC = DrawVsAng("cResEtaMC", "eta residual at vertex in 3 MC angular regions", hResEtaVertex, hResEtaVertexVsMCAngle);
855 TCanvas* cResEtaVsPos = DrawVsPos("cResEtaVsPos", "eta residual at vertex versus position at absorber end in 3 MC angular regions",
856 hResEtaVertexVsPosAbsEnd_0_2_DegMC, hResEtaVertexVsPosAbsEnd_2_3_DegMC, hResEtaVertexVsPosAbsEnd_3_10_DegMC);
857
858 // diplay phi residuals
859 TCanvas* cResPhi = DrawVsAng("cResPhi", "phi residual at vertex in 3 angular regions", hResPhiVertex, hResPhiVertexVsAngle);
860 TCanvas* cResPhiMC = DrawVsAng("cResPhiMC", "phi residual at vertex in 3 MC angular regions", hResPhiVertex, hResPhiVertexVsMCAngle);
861 TCanvas* cResPhiVsPos = DrawVsPos("cResPhiVsPos", "phi residual at vertex versus position at absorber end in 3 MC angular regions",
862 hResPhiVertexVsPosAbsEnd_0_2_DegMC, hResPhiVertexVsPosAbsEnd_2_3_DegMC, hResPhiVertexVsPosAbsEnd_3_10_DegMC);
c687aa97 863
864 // ###################################### save histogram ###################################### //
b8dc484b 865 histoFile->Write();
c687aa97 866
867 histoFile->cd("momentumAtVertex");
868 gMeanResMomVertexVsMom->Write();
869 gMostProbResMomVertexVsMom->Write();
870 gSigmaResMomVertexVsMom->Write();
e1fe81b8 871 cResMom->Write();
872 cResMomMC->Write();
873 cResMomVsPos->Write();
874 cResMom_2_3_Deg->Write();
875 cResMom_3_10_Deg->Write();
876 cResMom_0_2_DegMC->Write();
877
878 histoFile->cd("slopesAtVertex");
879 gMeanResSlopeXVertexVsMom->Write();
880 gMeanResSlopeYVertexVsMom->Write();
881 gSigmaResSlopeXVertexVsMom->Write();
882 gSigmaResSlopeYVertexVsMom->Write();
883 cResSlopeX->Write();
884 cResSlopeY->Write();
885 cResSlopeXMC->Write();
886 cResSlopeYMC->Write();
887 cResSlopeXVsPos->Write();
888 cResSlopeYVsPos->Write();
889
890 histoFile->cd("DCA");
c231afb9 891 gMeanPDCAVsMom_2_3_Deg->Write();
e1fe81b8 892 gSigmaPDCAVsMom_2_3_Deg->Write();
c231afb9 893 gMeanPDCAVsMom_3_10_Deg->Write();
e1fe81b8 894 gSigmaPDCAVsMom_3_10_Deg->Write();
895 gMeanPMCSAngVsMom_2_3_Deg->Write();
896 gSigmaPMCSAngVsMom_2_3_Deg->Write();
897 gMeanPMCSAngVsMom_3_10_Deg->Write();
898 gSigmaPMCSAngVsMom_3_10_Deg->Write();
899 cPDCA->Write();
900 cPDCAMC->Write();
901 cPDCAVsPos->Write();
902
903 histoFile->cd("etaAtVertex");
904 gMeanResEtaVertexVsMom->Write();
905 gSigmaResEtaVertexVsMom->Write();
906 cResEta->Write();
907 cResEtaMC->Write();
908 cResEtaVsPos->Write();
909
910 histoFile->cd("phiAtVertex");
911 gMeanResPhiVertexVsMom->Write();
912 gSigmaResPhiVertexVsMom->Write();
913 cResPhi->Write();
914 cResPhiMC->Write();
915 cResPhiVsPos->Write();
c687aa97 916
917 histoFile->cd("momentumAtFirstCluster");
918 gMeanResMomFirstClusterVsMom->Write();
919 gSigmaResMomFirstClusterVsMom->Write();
920
e1fe81b8 921 histoFile->cd("slopesAtFirstCluster");
922 gMeanResSlopeXFirstClusterVsMom->Write();
923 gMeanResSlopeYFirstClusterVsMom->Write();
924 gSigmaResSlopeXFirstClusterVsMom->Write();
925 gSigmaResSlopeYFirstClusterVsMom->Write();
926
c687aa97 927 histoFile->cd("clusters");
f202486b 928 gResidualXPerChMean->Write();
929 gResidualXPerChSigma->Write();
930 gResidualYPerChMean->Write();
931 gResidualYPerChSigma->Write();
dd6a48d2 932 gResidualXPerDEMean->Write();
933 gResidualXPerDESigma->Write();
934 gResidualYPerDEMean->Write();
935 gResidualYPerDESigma->Write();
c687aa97 936
b8dc484b 937 histoFile->Close();
c687aa97 938
e1fe81b8 939 // ###################################### clean memory ###################################### //
940 delete cResMom;
941 delete cResMomMC;
942 delete cResMomVsPos;
943 delete cResMom_2_3_Deg;
944 delete cResMom_3_10_Deg;
945 delete cResMom_0_2_DegMC;
946 delete cResSlopeX;
947 delete cResSlopeY;
948 delete cResSlopeXMC;
949 delete cResSlopeYMC;
950 delete cResSlopeXVsPos;
951 delete cResSlopeYVsPos;
952 delete cPDCA;
953 delete cPDCAMC;
954 delete cPDCAVsPos;
955 delete cResEta;
956 delete cResEtaMC;
957 delete cResEtaVsPos;
958 delete cResPhi;
959 delete cResPhiMC;
960 delete cResPhiVsPos;
c687aa97 961
e1fe81b8 962 // ###################################### print statistics ###################################### //
963 printf("\n");
964 printf("nb of reconstructible tracks: %d \n", nReconstructibleTracks);
965 printf("nb of reconstructed tracks: %d \n", nReconstructedTracks);
966 printf("nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
967
968 aMCSMoy /= (Double_t) nMCS;
969 aMCS2Moy /= (Double_t) nMCS;
970 dMCSMoy /= (Double_t) nMCS;
971 dMCS2Moy /= (Double_t) nMCS;
972 adMCSMoy /= (Double_t) nMCS;
973 Double_t sigma2_ThetaMCS = aMCS2Moy - aMCSMoy*aMCSMoy;
974 Double_t sigma2_PosMCS = dMCS2Moy - dMCSMoy*dMCSMoy;
975 Double_t cov_ThetaPosMCS = - (adMCSMoy - aMCSMoy*dMCSMoy);
976 printf("\nmultiple scattering of tracks between 3 and 10 deg. at absorber end:\n");
977 printf(" sigma_ThetaMCS = %f\n", TMath::Sqrt(sigma2_ThetaMCS));
978 printf(" sigma_PosMCS = %f\n", TMath::Sqrt(sigma2_PosMCS));
979 printf(" cov_ThetaPosMCS = %f\n", cov_ThetaPosMCS);
980 printf(" --> sigma_DCA = %f\n", TMath::Sqrt(AliMUONConstants::AbsZEnd()*AliMUONConstants::AbsZEnd()*sigma2_ThetaMCS
981 - 2.*AliMUONConstants::AbsZEnd()*cov_ThetaPosMCS + sigma2_PosMCS));
dd6a48d2 982
983 printf("\nchamber resolution:\n");
984 printf(" - non-bending:");
985 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %5.3f":", %5.3f",clusterResPerCh[i][0]);
986 printf("\n - bending:");
987 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) printf((i==0)?" %6.4f":", %6.4f",clusterResPerCh[i][1]);
988 printf("\n\n");
989
990 printf("\nDE resolution:\n");
991 printf(" - non-bending:");
992 for (Int_t i = 0; i < deNBins; i++) printf((i==0)?" %5.3f":", %5.3f",clusterResPerDE[i][0]);
993 printf("\n - bending:");
994 for (Int_t i = 0; i < deNBins; i++) printf((i==0)?" %6.4f":", %6.4f",clusterResPerDE[i][1]);
995 printf("\n\n");
c687aa97 996}
997
998//------------------------------------------------------------------------------------
999Double_t langaufun(Double_t *x, Double_t *par) {
1000
1001 //Fit parameters:
1002 //par[0]=Width (scale) parameter of Landau density
1003 //par[1]=Most Probable (MP, location) parameter of Landau density
1004 //par[2]=Total area (integral -inf to inf, normalization constant)
1005 //par[3]=Width (sigma) of convoluted Gaussian function
1006 //
1007 //In the Landau distribution (represented by the CERNLIB approximation),
1008 //the maximum is located at x=-0.22278298 with the location parameter=0.
1009 //This shift is corrected within this function, so that the actual
1010 //maximum is identical to the MP parameter.
1011
1012 // Numeric constants
1013 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
1014 Double_t mpshift = -0.22278298; // Landau maximum location
1015
1016 // Control constants
1017 Double_t np = 100.0; // number of convolution steps
1018 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
1019
1020 // Variables
1021 Double_t xx;
1022 Double_t mpc;
1023 Double_t fland;
1024 Double_t sum = 0.0;
1025 Double_t xlow,xupp;
1026 Double_t step;
1027 Double_t i;
1028
1029
1030 // MP shift correction
1031 mpc = par[1] - mpshift * par[0];
1032
1033 // Range of convolution integral
1034 xlow = x[0] - sc * par[3];
1035 xupp = x[0] + sc * par[3];
1036
1037 step = (xupp-xlow) / np;
1038
1039 // Convolution integral of Landau and Gaussian by sum
1040 for(i=1.0; i<=np/2; i++) {
1041 xx = xlow + (i-.5) * step;
1042 //change x -> -x because the tail of the Landau is at the left here...
1043 fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
1044 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1045
1046 //change x -> -x because the tail of the Landau is at the left here...
1047 xx = xupp - (i-.5) * step;
1048 fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
1049 sum += fland * TMath::Gaus(x[0],xx,par[3]);
1050 }
1051
1052 return (par[2] * step * sum * invsq2pi / par[3]);
b8dc484b 1053}
1054
e1fe81b8 1055//------------------------------------------------------------------------------------
1056void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0,
1057 const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
1058{
1059 /// generic function to fit residuals versus momentum with a gaussian
1060 static TF1* fGaus = 0x0;
1061 if (!fGaus) fGaus = new TF1("fGaus","gaus");
1062
1063 Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1064 for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1065 cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1066 TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
1067 fGaus->SetParameters(tmp->GetEntries(), mean0, sigma0);
1068 tmp->Fit("fGaus","WWNQ");
1069 Int_t rebin = TMath::Max(Int_t(0.5*fGaus->GetParameter(2)/tmp->GetBinWidth(1)),1);
1070 while (tmp->GetNbinsX()%rebin!=0) rebin--;
1071 tmp->Rebin(rebin);
1072 tmp->Fit("fGaus","NQ");
1073 h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
248ee365 1074 Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
e1fe81b8 1075 h->GetXaxis()->SetRange();
1076 Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
1077 gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
1078 gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(1), fGaus->GetParError(1));
1079 gSigma->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(2));
1080 gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(2), fGaus->GetParError(2));
1081 delete tmp;
1082 }
1083 cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1084}
1085
1086//------------------------------------------------------------------------------------
c231afb9 1087void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
e1fe81b8 1088{
1089 /// generic function to fit p*DCA distributions
1090 static TF1* fPGaus = 0x0;
1091 if (!fPGaus) fPGaus = new TF1("fPGaus","x*gaus");
1092
1093 Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1094 for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1095 cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1096 TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
c231afb9 1097 fPGaus->SetParameters(1.,0.,80.);
1098 Int_t rebin = 25.*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1)));
e1fe81b8 1099 while (tmp->GetNbinsX()%rebin!=0) rebin--;
1100 tmp->Rebin(rebin);
1101 tmp->Fit("fPGaus","NQ");
1102 h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
248ee365 1103 Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
e1fe81b8 1104 h->GetXaxis()->SetRange();
1105 Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
c231afb9 1106 gMean->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(1));
1107 gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(1), fPGaus->GetParError(1));
e1fe81b8 1108 gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
1109 gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
1110 delete tmp;
1111 }
1112 cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1113}
1114
1115//------------------------------------------------------------------------------------
dd6a48d2 1116void FillResidual(TH1* h, Int_t i, Double_t& sigma, TGraphErrors* gMean, TGraphErrors* gSigma, Bool_t correctForSystematics, Bool_t fitResiduals)
1117{
1118 /// fill graphs with residual mean and sigma
1119 static TF1* fRGaus = 0x0;
1120 Double_t mean, meanErr, sigmaErr;
1121
1122 if (fitResiduals) {
1123
1124 if (!fRGaus) fRGaus = new TF1("fRGaus","gaus");
1125 fRGaus->SetParameters(h->GetEntries(), 0., 0.1);
1126 h->Fit("fRGaus", "WWNQ");
1127 mean = fRGaus->GetParameter(1);
1128 meanErr = fRGaus->GetParError(1);
1129 sigma = fRGaus->GetParameter(2);
1130 sigmaErr = fRGaus->GetParError(2);
1131
1132 } else {
1133
1134 Zoom(h);
1135 mean = h->GetMean();
1136 meanErr = h->GetMeanError();
1137 sigma = h->GetRMS();
1138 sigmaErr = h->GetRMSError();
1139 h->GetXaxis()->SetRange(0,0);
1140
1141 }
1142
1143 gMean->SetPoint(i, i+1, mean);
1144 gMean->SetPointError(i, 0., meanErr);
1145 if (correctForSystematics) {
1146 Double_t s = TMath::Sqrt(mean*mean + sigma*sigma);
1147 sigmaErr = (s>0.) ? TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + mean*mean*meanErr*meanErr) / s : 0.;
1148 sigma = s;
1149 }
1150 gSigma->SetPoint(i, i+1, sigma);
1151 gSigma->SetPointError(i, 0., sigmaErr);
1152}
1153
1154//------------------------------------------------------------------------------------
e1fe81b8 1155TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2)
1156{
1157 /// generic function to draw histograms versus absorber angular region
1158 TCanvas* c = new TCanvas(name, title);
1159 c->cd();
1160 h1->Draw();
1161 TH1D *proj1 = h2->ProjectionY(Form("%s_proj_0_2",h2->GetName()),1,2);
1162 proj1->Draw("sames");
1163 proj1->SetLineColor(2);
1164 TH1D *proj2 = h2->ProjectionY(Form("%s_proj_2_3",h2->GetName()),3,3);
1165 proj2->Draw("sames");
1166 proj2->SetLineColor(4);
1167 TH1D *proj3 = h2->ProjectionY(Form("%s__proj_3_10",h2->GetName()),4,10);
1168 proj3->Draw("sames");
1169 proj3->SetLineColor(3);
1170 return c;
1171}
1172
1173//------------------------------------------------------------------------------------
1174TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3)
1175{
1176 /// generic function to draw histograms versus position at absorber end
1177 TCanvas* c = new TCanvas(name, title);
1178 c->cd();
1179 h1->Draw();
1180 h1->SetMarkerColor(2);
1181 h2->Draw("sames");
1182 h2->SetMarkerColor(4);
1183 h3->Draw("sames");
1184 h3->SetMarkerColor(3);
1185 return c;
1186}
1187
1188//------------------------------------------------------------------------------------
1189TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2, const char* fitting)
1190{
1191 /// generic function to draw and eventually fit momentum residuals versus momentum
1192 TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
1193 TCanvas* c = new TCanvas(name, title);
1194 c->cd();
1195 TH1D* proj = 0x0;
1196 h->Sumw2();
1197 Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
1198 for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1199 if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
1200 proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
1201 if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
1202 proj->Draw((i==rebinFactorX)?"hist":"histsames");
1203 proj->SetLineColor(i/rebinFactorX);
1204 if (f2) {
1205 f2->SetParameters(0.2,0.,1.,1.);
1206 f2->SetLineColor(i/rebinFactorX);
1207 proj->Fit("f2","WWNQ","sames");
1208 Double_t fwhm = f2->GetParameter(0);
1209 Double_t sigma = f2->GetParameter(3);
1210 Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
1211 Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
1212 while (proj->GetNbinsX()%rebin!=0) rebin--;
1213 proj->Rebin(rebin);
1214 proj->Scale(1./rebin);
1215 proj->Fit("f2","Q","sames");
1216 } else proj->SetLineWidth(2);
1217 Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
1218 l->AddEntry(proj,Form("%5.1f GeV",p));
1219 }
1220 if (f2) cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
1221 l->Draw("same");
1222 return c;
1223}
1224
dd6a48d2 1225//------------------------------------------------------------------------------------
1226void Zoom(TH1* h, Double_t fractionCut)
1227{
1228 /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
1229 Int_t maxEventsCut = fractionCut * h->GetEntries();
1230 Int_t nBins = h->GetNbinsX();
1231
1232 // set low edge
1233 Int_t minBin;
1234 Int_t eventsCut = 0;
1235 for (minBin = 1; minBin <= nBins; minBin++) {
1236 eventsCut += h->GetBinContent(minBin);
1237 if (eventsCut > maxEventsCut) break;
1238 }
1239
1240 // set high edge
1241 Int_t maxBin;
1242 eventsCut = 0;
1243 for (maxBin = nBins; maxBin >= 1; maxBin--) {
1244 eventsCut += h->GetBinContent(maxBin);
1245 if (eventsCut > maxEventsCut) break;
1246 }
1247
1248 // set new axis range
1249 h->GetXaxis()->SetRange(--minBin, ++maxBin);
1250}
1251