]>
Commit | Line | Data |
---|---|---|
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 | 63 | Double_t langaufun(Double_t *x, Double_t *par); |
e1fe81b8 | 64 | void FitGausResVsMom(TH2* h, Int_t nBins, const Double_t mean0, const Double_t sigma0, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma); |
c231afb9 | 65 | void FitPDCAVsMom(TH2* h, Int_t nBins, const char* fitting, TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma); |
dd6a48d2 | 66 | void FillResidual(TH1* h, Int_t i, Double_t& sigma, TGraphErrors* gMean, TGraphErrors* gSigma, Bool_t correctForSystematics, Bool_t fitResiduals); |
e1fe81b8 | 67 | TCanvas* DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2); |
68 | TCanvas* DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3); | |
69 | TCanvas* DrawResMomVsMom(const char* name, const char* title, TH2* h, Int_t nBins, TF1* f2 = 0x0, const char* fitting = ""); | |
dd6a48d2 | 70 | void Zoom(TH1* h, Double_t fractionCut = 0.01); |
c687aa97 | 71 | |
72 | //------------------------------------------------------------------------------------ | |
c231afb9 | 73 | void 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 | //------------------------------------------------------------------------------------ | |
999 | Double_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 | //------------------------------------------------------------------------------------ |
1056 | void 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 | 1087 | void 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 | ||
dd6a48d2 | 1115 | //------------------------------------------------------------------------------------ |
1116 | void 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 | ||
e1fe81b8 | 1154 | //------------------------------------------------------------------------------------ |
1155 | TCanvas* 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 | //------------------------------------------------------------------------------------ | |
1174 | TCanvas* 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 | //------------------------------------------------------------------------------------ | |
1189 | TCanvas* 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 | //------------------------------------------------------------------------------------ |
1226 | void 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 |