]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONRecoCheck.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Useful macros for the shifter
[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
f202486b 29#include "TMath.h"
b8dc484b 30#include "TClonesArray.h"
31#include "TH1.h"
f202486b 32#include "TH2.h"
c687aa97 33#include "TH3.h"
f202486b 34#include "TGraphErrors.h"
35#include "TF1.h"
b8dc484b 36#include "TFile.h"
c687aa97 37#include "TCanvas.h"
38#include "TLegend.h"
b8dc484b 39
40// STEER includes
a99c3449 41#include "AliCDBManager.h"
c687aa97 42#include "AliLog.h"
b8dc484b 43
44// MUON includes
a99c3449 45#include "AliMUONCDB.h"
b8dc484b 46#include "AliMUONConstants.h"
47#include "AliMUONTrack.h"
48#include "AliMUONRecoCheck.h"
49#include "AliMUONTrackParam.h"
a99c3449 50#include "AliMUONRecoParam.h"
48217459 51#include "AliMUONVTrackStore.h"
f202486b 52#include "AliMUONVCluster.h"
c687aa97 53#include "AliMUONTrackExtrap.h"
b8dc484b 54
c687aa97 55Double_t langaufun(Double_t *x, Double_t *par);
56
57//------------------------------------------------------------------------------------
58void MUONRecoCheck (Int_t nEvent = -1, const char* pathSim="./generated/", const char* esdFileName="AliESDs.root",
f202486b 59 const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
e54bf126 60{
b8dc484b 61
c687aa97 62 AliLog::SetClassDebugLevel("AliMCEvent",-1);
63
64 // ###################################### define histograms ###################################### //
b8dc484b 65 // File for histograms and histogram booking
66 TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
48217459 67
b8dc484b 68 TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks ",15,-0.5,14.5);
69 TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
96ebe67e 70 TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5);
b8dc484b 71 TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
72
c687aa97 73 // momentum resolution at vertex
74 histoFile->mkdir("momentumAtVertex","momentumAtVertex");
75 histoFile->cd("momentumAtVertex");
76
77 const Int_t pNBins = 30;
78 const Double_t pEdges[2] = {0., 300.};
79 const Int_t deltaPAtVtxNBins = 250;
80 const Double_t deltaPAtVtxEdges[2] = {-35., 15.};
81
82 TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
83
84 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]);
85 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]);
86 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]);
87 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]);
88
89 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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
90 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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
91 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)",100,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
92
93 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]);
94 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]);
95 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]);
96
97 TGraphErrors* gMeanResMomVertexVsMom = new TGraphErrors(pNBins);
98 gMeanResMomVertexVsMom->SetName("gMeanResMomVertexVsMom");
99 gMeanResMomVertexVsMom->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
100 TGraphErrors* gMostProbResMomVertexVsMom = new TGraphErrors(pNBins);
101 gMostProbResMomVertexVsMom->SetName("gMostProbResMomVertexVsMom");
102 gMostProbResMomVertexVsMom->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
103 TGraphErrors* gSigmaResMomVertexVsMom = new TGraphErrors(pNBins);
104 gSigmaResMomVertexVsMom->SetName("gSigmaResMomVertexVsMom");
105 gSigmaResMomVertexVsMom->SetTitle("#Delta_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
106 TF1 *f2 = new TF1("f2",langaufun,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1],4);
107
108 // momentum resolution at first cluster
109 histoFile->mkdir("momentumAtFirstCluster","momentumAtFirstCluster");
110 histoFile->cd("momentumAtFirstCluster");
111
112 const Int_t deltaPAtFirstClNBins = 500;
113 const Double_t deltaPAtFirstClEdges[2] = {-25., 25.};
114
115 TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
116 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]);
117
118 TGraphErrors* gMeanResMomFirstClusterVsMom = new TGraphErrors(pNBins);
119 gMeanResMomFirstClusterVsMom->SetName("gMeanResMomFirstClusterVsMom");
120 gMeanResMomFirstClusterVsMom->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
121 TGraphErrors* gSigmaResMomFirstClusterVsMom = new TGraphErrors(pNBins);
122 gSigmaResMomFirstClusterVsMom->SetName("gSigmaResMomFirstClusterVsMom");
123 gSigmaResMomFirstClusterVsMom->SetTitle("#Delta_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
124 TF1* f = new TF1("f","gausn");
125
126 // cluster resolution
127 histoFile->mkdir("clusters","clusters");
128 histoFile->cd("clusters");
129
f202486b 130 TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()];
131 TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()];
132 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
133 hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), 1000, -1., 1.);
134 hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), 1000, -0.5, 0.5);
135 }
c687aa97 136
f202486b 137 TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
138 gResidualXPerChMean->SetName("gResidualXPerChMean");
c687aa97 139 gResidualXPerChMean->SetTitle("cluster-track residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
f202486b 140 gResidualXPerChMean->SetMarkerStyle(kFullDotLarge);
141 TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh());
142 gResidualYPerChMean->SetName("gResidualYPerChMean");
143 gResidualYPerChMean->SetTitle("cluster-track residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
144 gResidualYPerChMean->SetMarkerStyle(kFullDotLarge);
145 TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
146 gResidualXPerChSigma->SetName("gResidualXPerChSigma");
147 gResidualXPerChSigma->SetTitle("cluster-track residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
148 gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge);
149 TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh());
150 gResidualYPerChSigma->SetName("gResidualYPerChSigma");
c687aa97 151 gResidualYPerChSigma->SetTitle("cluster-track residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
f202486b 152 gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge);
48217459 153
c687aa97 154 // ###################################### initialize ###################################### //
a99c3449 155 AliMUONRecoCheck rc(esdFileName, pathSim);
156
157 // load necessary data from OCDB
158 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
159 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
c687aa97 160 if (!AliMUONCDB::LoadField()) return;
161 AliMUONTrackExtrap::SetField();
a99c3449 162 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
163 if (!recoParam) return;
a0bfad0e 164
f202486b 165 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
a99c3449 166 Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
f202486b 167 // compute the mask of requested stations from recoParam
168 UInt_t requestedStationMask = 0;
169 for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
170 // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
171 Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
48217459 172
173 Int_t nevents = rc.NumberOfEvents();
b8dc484b 174
edcc5b56 175 if (nevents < nEvent || nEvent < 0) nEvent = nevents;
b8dc484b 176
177 Int_t ievent;
178 Int_t nReconstructibleTracks = 0;
179 Int_t nReconstructedTracks = 0;
180 Int_t nReconstructibleTracksCheck = 0;
f202486b 181 AliMUONTrackParam *trackParam;
182 Double_t x1,y1,z1,pX1,pY1,pZ1,p1,pT1;
183 Double_t x2,y2,z2,pX2,pY2,pZ2,p2,pT2;
c687aa97 184 Double_t xAbs,yAbs,dAbs,aAbs,aMC;
48217459 185
c687aa97 186 // ###################################### fill histograms ###################################### //
48217459 187 for (ievent=0; ievent<nEvent; ievent++)
188 {
c687aa97 189 if ((ievent+1)%100 == 0) cout<<"\rEvent processing... "<<ievent+1<<flush;
b8dc484b 190
a99c3449 191 AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
f202486b 192 AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
b8dc484b 193
48217459 194 hReconstructible->Fill(trackRefStore->GetSize());
195 hReco->Fill(trackStore->GetSize());
196
197 nReconstructibleTracks += trackRefStore->GetSize();
198 nReconstructedTracks += trackStore->GetSize();
199
f202486b 200 // loop over trackRef
48217459 201 TIter next(trackRefStore->CreateIterator());
202 AliMUONTrack* trackRef;
48217459 203 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) )
204 {
48217459 205
f202486b 206 hTrackRefID->Fill(trackRef->GetMCLabel());
207
208 AliMUONTrack* trackMatched = 0x0;
209 Int_t nMatchClusters = 0;
210
211 // loop over trackReco and look for compatible track
48217459 212 TIter next2(trackStore->CreateIterator());
213 AliMUONTrack* trackReco;
48217459 214 while ( ( trackReco = static_cast<AliMUONTrack*>(next2()) ) )
215 {
f202486b 216
217 // check if trackReco is compatible with trackRef
218 Int_t n = 0;
219 if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) {
220 trackMatched = trackReco;
221 nMatchClusters = n;
222 break;
61fed964 223 }
f202486b 224
48217459 225 }
226
f202486b 227 if (trackMatched) { // tracking requirements verified, track is found
48217459 228 nReconstructibleTracksCheck++;
f202486b 229 hNClusterComp->Fill(nMatchClusters);
c687aa97 230
231 // compute track position at the end of the absorber
232 AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)trackMatched->GetTrackParamAtCluster()->First()));
233 AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
234 xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
235 yAbs = trackParamAtAbsEnd.GetBendingCoor();
236 dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
237 aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
238
48217459 239 trackParam = trackRef->GetTrackParamAtVertex();
240 x1 = trackParam->GetNonBendingCoor();
241 y1 = trackParam->GetBendingCoor();
242 z1 = trackParam->GetZ();
243 pX1 = trackParam->Px();
244 pY1 = trackParam->Py();
245 pZ1 = trackParam->Pz();
246 p1 = trackParam->P();
f202486b 247 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
c687aa97 248 aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
249
f202486b 250 trackParam = trackMatched->GetTrackParamAtVertex();
48217459 251 x2 = trackParam->GetNonBendingCoor();
252 y2 = trackParam->GetBendingCoor();
253 z2 = trackParam->GetZ();
254 pX2 = trackParam->Px();
255 pY2 = trackParam->Py();
256 pZ2 = trackParam->Pz();
257 p2 = trackParam->P();
f202486b 258 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
48217459 259
260 hResMomVertex->Fill(p2-p1);
f202486b 261 hResMomVertexVsMom->Fill(p1,p2-p1);
c687aa97 262 hResMomVertexVsAngleVsMom->Fill(p1,aAbs,p2-p1);
263 if (aAbs > 2. && aAbs < 3.) hResMomVertexVsMom_2_3_Deg->Fill(p1,p2-p1);
264 else if (aAbs >= 3. && aAbs < 10.) hResMomVertexVsMom_3_10_Deg->Fill(p1,p2-p1);
265 if (aMC < 2.) {
266 hResMomVertexVsMom_0_2_DegMC->Fill(p1,p2-p1);
267 hResMomVertexVsPosAbsEnd_0_2_DegMC->Fill(dAbs,p2-p1);
268 }
269 else if (aMC >= 2. && aMC < 3) hResMomVertexVsPosAbsEnd_2_3_DegMC->Fill(dAbs,p2-p1);
270 else if (aMC >= 3. && aMC < 10.) hResMomVertexVsPosAbsEnd_3_10_DegMC->Fill(dAbs,p2-p1);
271 hResMomVertexVsAngle->Fill(aAbs,p2-p1);
272 hResMomVertexVsMCAngle->Fill(aMC,p2-p1);
f202486b 273
96ebe67e 274 trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
48217459 275 x1 = trackParam->GetNonBendingCoor();
276 y1 = trackParam->GetBendingCoor();
277 z1 = trackParam->GetZ();
278 pX1 = trackParam->Px();
279 pY1 = trackParam->Py();
280 pZ1 = trackParam->Pz();
281 p1 = trackParam->P();
f202486b 282 pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
283
f202486b 284 trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
48217459 285 x2 = trackParam->GetNonBendingCoor();
286 y2 = trackParam->GetBendingCoor();
287 z2 = trackParam->GetZ();
288 pX2 = trackParam->Px();
289 pY2 = trackParam->Py();
290 pZ2 = trackParam->Pz();
291 p2 = trackParam->P();
f202486b 292 pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
48217459 293
96ebe67e 294 hResMomFirstCluster->Fill(p2-p1);
f202486b 295 hResMomFirstClusterVsMom->Fill(p1,p2-p1);
296
297 // Fill residuals
298 // Loop over clusters of first track
299 AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First();
300 while (trackParamAtCluster1) {
301 AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr();
302 AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First();
303 while (trackParamAtCluster2) {
304 AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr();
305 if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) {
306 hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX());
307 hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY());
308 break;
309 }
310 trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2);
311 }
312 trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1);
313 }
314
b8dc484b 315 }
f202486b 316
b8dc484b 317 } // end loop track ref.
b8dc484b 318
48217459 319 } // end loop on event
c687aa97 320 cout<<"\rEvent processing... "<<nevents<<" done"<<endl;
48217459 321
c687aa97 322 // ###################################### compute stuff ###################################### //
323 // compute momentum resolution at vertex versus p
324 Int_t rebinFactorX = TMath::Max(hResMomVertexVsMom->GetNbinsX()/pNBins, 1);
325 for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom->GetNbinsX(); i+=rebinFactorX) {
326 cout<<"\rFitting momentum residuals at vertex... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
327 TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
328 Double_t p = 0.5 * (hResMomVertexVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom->GetBinLowEdge(i+1));
329 f2->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
330 tmp->Fit("f2","WWNQ");
331 Double_t fwhm = f2->GetParameter(0);
332 Double_t sigma = f2->GetParameter(3);
333 Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
334 Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/tmp->GetBinWidth(1)),1);
335 while (deltaPAtVtxNBins%rebin!=0) rebin--;
336 tmp->Rebin(rebin);
337 tmp->Fit("f2","NQ");
338 fwhm = f2->GetParameter(0);
339 sigma = f2->GetParameter(3);
340 sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
341 Double_t fwhmErr = f2->GetParError(0);
342 Double_t sigmaErr = f2->GetParError(3);
343 Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
344 gMeanResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,tmp->GetMean());
345 gMeanResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),tmp->GetMeanError());
346 gMostProbResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,-f2->GetParameter(1));
347 gMostProbResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),f2->GetParError(1));
348 gSigmaResMomVertexVsMom->SetPoint(i/rebinFactorX-1,p,100.*sigmaP/p);
349 gSigmaResMomVertexVsMom->SetPointError(i/rebinFactorX-1,hResMomVertexVsMom->GetBinWidth(i),100.*sigmaPErr/p);
f202486b 350 delete tmp;
351 }
c687aa97 352 cout<<"\rFitting momentum residuals at vertex... "<<pNBins<<"/"<<pNBins<<endl;
353
354 // compute momentum resolution at first cluster versus p
355 rebinFactorX = TMath::Max(hResMomFirstClusterVsMom->GetNbinsX()/pNBins, 1);
356 for (Int_t i = rebinFactorX; i <= hResMomFirstClusterVsMom->GetNbinsX(); i+=rebinFactorX) {
357 cout<<"\rFitting momentum residuals at first cluster... "<<i/rebinFactorX<<"/"<<pNBins<<flush;
358 TH1D *tmp = hResMomFirstClusterVsMom->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
359 Double_t p = 0.5 * (hResMomFirstClusterVsMom->GetBinLowEdge(i-rebinFactorX+1) + hResMomFirstClusterVsMom->GetBinLowEdge(i+1));
360 f->SetParameters(tmp->GetEntries(),0.,1.);
361 tmp->Fit("f","WWNQ");
362 Int_t rebin = TMath::Max(Int_t(0.5*f->GetParameter(2)/tmp->GetBinWidth(1)),1);
363 while (deltaPAtFirstClNBins%rebin!=0) rebin--;
364 tmp->Rebin(rebin);
365 tmp->Fit("f","NQ");
366 gMeanResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1,p,f->GetParameter(1));
367 gMeanResMomFirstClusterVsMom->SetPointError(i/rebinFactorX-1,hResMomFirstClusterVsMom->GetBinWidth(i),f->GetParError(1));
368 gSigmaResMomFirstClusterVsMom->SetPoint(i/rebinFactorX-1,p,100.*f->GetParameter(2)/p);
369 gSigmaResMomFirstClusterVsMom->SetPointError(i/rebinFactorX-1,hResMomFirstClusterVsMom->GetBinWidth(i),100.*f->GetParError(2)/p);
f202486b 370 delete tmp;
371 }
c687aa97 372 cout<<"\rFitting momentum residuals at first cluster... "<<pNBins<<"/"<<pNBins<<endl;
f202486b 373
374 // compute residual mean and dispersion
375 for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
376 hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS());
377 gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean());
378 gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError());
379 gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS());
380 gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError());
381 hResidualXInCh[i]->GetXaxis()->SetRange(0,0);
382 hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS());
383 gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean());
384 gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError());
385 gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS());
386 gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError());
387 hResidualYInCh[i]->GetXaxis()->SetRange(0,0);
388 }
389
c687aa97 390 // ###################################### display histograms ###################################### //
391 // diplay momentum residual for different angular region
392 TCanvas cResMom("cResMom", "momentum residual at vertex in 3 angular regions");
393 cResMom.cd();
394 hResMomVertex->Draw();
395 TH1D *hResMomVertex_0_2_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_0_2_Deg",1,2);
396 hResMomVertex_0_2_Deg->Draw("sames");
397 hResMomVertex_0_2_Deg->SetLineColor(2);
398 TH1D *hResMomVertex_2_3_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_2_3_Deg",3,3);
399 hResMomVertex_2_3_Deg->Draw("sames");
400 hResMomVertex_2_3_Deg->SetLineColor(4);
401 TH1D *hResMomVertex_3_10_Deg = hResMomVertexVsAngle->ProjectionY("hResMomVertex_3_10_Deg",4,10);
402 hResMomVertex_3_10_Deg->Draw("sames");
403 hResMomVertex_3_10_Deg->SetLineColor(3);
404
405 // diplay momentum residual for different angular region
406 TCanvas cResMomMC("cResMomMC", "momentum residual at vertex in 3 MC angular regions");
407 cResMomMC.cd();
408 hResMomVertex->Draw();
409 TH1D *hResMomVertex_0_2_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_0_2_DegMC",1,2);
410 hResMomVertex_0_2_DegMC->Draw("sames");
411 hResMomVertex_0_2_DegMC->SetLineColor(2);
412 TH1D *hResMomVertex_2_3_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_2_3_DegMC",3,3);
413 hResMomVertex_2_3_DegMC->Draw("sames");
414 hResMomVertex_2_3_DegMC->SetLineColor(4);
415 TH1D *hResMomVertex_3_10_DegMC = hResMomVertexVsMCAngle->ProjectionY("hResMomVertex_3_10_DegMC",4,10);
416 hResMomVertex_3_10_DegMC->Draw("sames");
417 hResMomVertex_3_10_DegMC->SetLineColor(3);
418
419 // diplay momentum residual versus position at absorber end for different MC angular region
420 TCanvas cResMomVsPos("cResMomVsPos", "momentum residual at vertex versus position at absorber end in 3 MC angular regions");
421 cResMomVsPos.cd();
422 hResMomVertexVsPosAbsEnd_0_2_DegMC->Draw();
423 hResMomVertexVsPosAbsEnd_0_2_DegMC->SetMarkerColor(2);
424 hResMomVertexVsPosAbsEnd_2_3_DegMC->Draw("sames");
425 hResMomVertexVsPosAbsEnd_2_3_DegMC->SetMarkerColor(4);
426 hResMomVertexVsPosAbsEnd_3_10_DegMC->Draw("sames");
427 hResMomVertexVsPosAbsEnd_3_10_DegMC->SetMarkerColor(3);
428
429 // diplay momentum residual of tracks between 2 and 3 deg. for different momentum values
430 Int_t pNBinsShown = 10;
431 TLegend lResMom_2_3_Deg(0.15,0.25,0.3,0.85);
432 TCanvas cResMom_2_3_Deg("cResMom_2_3_Deg", "momentum residual for tracks between 2 and 3 degrees");
433 cResMom_2_3_Deg.cd();
434 TH1D* proj = 0x0;
435 hResMomVertexVsMom_2_3_Deg->Sumw2();
436 rebinFactorX = TMath::Max(hResMomVertexVsMom_2_3_Deg->GetNbinsX()/pNBinsShown, 1);
437 for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_2_3_Deg->GetNbinsX(); i+=rebinFactorX) {
438 cout<<"\rFitting momentum residuals at vertex (tracks in [2,3] deg.)... "<<i/rebinFactorX<<"/"<<pNBinsShown<<flush;
439 proj = hResMomVertexVsMom_2_3_Deg->ProjectionY(Form("hRes23_%d",i/rebinFactorX),i-rebinFactorX+1,i);
440 if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
441 proj->Draw((i==rebinFactorX)?"hist":"histsames");
442 proj->SetLineColor(i/rebinFactorX);
443 f2->SetParameters(0.2,0.,1.,1.);
444 f2->SetLineColor(i/rebinFactorX);
445 proj->Fit("f2","WWNQ","sames");
446 Double_t fwhm = f2->GetParameter(0);
447 Double_t sigma = f2->GetParameter(3);
448 Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
449 Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
450 while (deltaPAtVtxNBins%rebin!=0) rebin--;
451 proj->Rebin(rebin);
452 proj->Scale(1./rebin);
453 proj->Fit("f2","Q","sames");
454 Double_t p = 0.5 * (hResMomVertexVsMom_2_3_Deg->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_2_3_Deg->GetBinLowEdge(i+1));
455 lResMom_2_3_Deg.AddEntry(proj,Form("%5.1f GeV",p));
456 }
457 cout<<"\rFitting momentum residuals at vertex (tracks in [2,3] deg.)... "<<pNBinsShown<<"/"<<pNBinsShown<<endl;
458 lResMom_2_3_Deg.Draw("same");
b8dc484b 459
c687aa97 460 // diplay momentum residual of tracks between 3 and 10 deg. for different momentum values
461 pNBinsShown = 10;
462 TLegend lResMom_3_10_Deg(0.15,0.25,0.3,0.85);
463 TCanvas cResMom_3_10_Deg("cResMom_3_10_Deg", "momentum residual for tracks between 3 and 10 degrees");
464 cResMom_3_10_Deg.cd();
465 proj = 0x0;
466 hResMomVertexVsMom_3_10_Deg->Sumw2();
467 rebinFactorX = TMath::Max(hResMomVertexVsMom_3_10_Deg->GetNbinsX()/pNBinsShown, 1);
468 for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_3_10_Deg->GetNbinsX(); i+=rebinFactorX) {
469 cout<<"\rFitting momentum residuals at vertex (tracks in [3,10] deg.)... "<<i/rebinFactorX<<"/"<<pNBinsShown<<flush;
470 proj = hResMomVertexVsMom_3_10_Deg->ProjectionY(Form("hRes310_%d",i/rebinFactorX),i-rebinFactorX+1,i);
471 if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
472 proj->Draw((i==rebinFactorX)?"hist":"histsames");
473 proj->SetLineColor(i/rebinFactorX);
474 f2->SetParameters(0.2,0.,1.,1.);
475 f2->SetLineColor(i/rebinFactorX);
476 proj->Fit("f2","WWNQ","sames");
477 Double_t fwhm = f2->GetParameter(0);
478 Double_t sigma = f2->GetParameter(3);
479 Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
480 Int_t rebin = TMath::Max(Int_t(0.5*sigmaP/proj->GetBinWidth(1)),1);
481 while (deltaPAtVtxNBins%rebin!=0) rebin--;
482 proj->Rebin(rebin);
483 proj->Scale(1./rebin);
484 proj->Fit("f2","Q","sames");
485 Double_t p = 0.5 * (hResMomVertexVsMom_3_10_Deg->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_3_10_Deg->GetBinLowEdge(i+1));
486 lResMom_3_10_Deg.AddEntry(proj,Form("%5.1f GeV",p));
487 }
488 cout<<"\rFitting momentum residuals at vertex (tracks in [3,10] deg.)... "<<pNBinsShown<<"/"<<pNBinsShown<<endl;
489 lResMom_3_10_Deg.Draw("same");
490
491 // diplay momentum residuals of tracks with MC angle < 2 deg. for different momentum values
492 pNBinsShown = 5;
493 TLegend lResMom_0_2_DegMC(0.15,0.25,0.3,0.85);
494 TCanvas cResMom_0_2_DegMC("cResMom_0_2_DegMC", "momentum residuals for tracks with MC angle < 2 degrees");
495 cResMom_0_2_DegMC.cd();
496 proj = 0x0;
497 hResMomVertexVsMom_0_2_DegMC->Sumw2();
498 rebinFactorX = TMath::Max(hResMomVertexVsMom_0_2_DegMC->GetNbinsX()/pNBinsShown, 1);
499 for (Int_t i = rebinFactorX; i <= hResMomVertexVsMom_0_2_DegMC->GetNbinsX(); i+=rebinFactorX) {
500 proj = hResMomVertexVsMom_0_2_DegMC->ProjectionY(Form("hRes02_%d",i/rebinFactorX),i-rebinFactorX+1,i);
501 if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
502 proj->Draw((i==rebinFactorX)?"hist":"histsames");
503 proj->SetLineColor(i/rebinFactorX);
504 proj->SetLineWidth(2);
505 Double_t p = 0.5 * (hResMomVertexVsMom_0_2_DegMC->GetBinLowEdge(i-rebinFactorX+1) + hResMomVertexVsMom_0_2_DegMC->GetBinLowEdge(i+1));
506 lResMom_0_2_DegMC.AddEntry(proj,Form("%5.1f GeV",p));
507 }
508 lResMom_0_2_DegMC.Draw("same");
509
510 // ###################################### save histogram ###################################### //
b8dc484b 511 histoFile->Write();
c687aa97 512
513 histoFile->cd("momentumAtVertex");
514 gMeanResMomVertexVsMom->Write();
515 gMostProbResMomVertexVsMom->Write();
516 gSigmaResMomVertexVsMom->Write();
517 cResMom.Write();
518 cResMomMC.Write();
519 cResMomVsPos.Write();
520 cResMom_2_3_Deg.Write();
521 cResMom_3_10_Deg.Write();
522 cResMom_0_2_DegMC.Write();
523
524 histoFile->cd("momentumAtFirstCluster");
525 gMeanResMomFirstClusterVsMom->Write();
526 gSigmaResMomFirstClusterVsMom->Write();
527
528 histoFile->cd("clusters");
f202486b 529 gResidualXPerChMean->Write();
530 gResidualXPerChSigma->Write();
531 gResidualYPerChMean->Write();
532 gResidualYPerChSigma->Write();
c687aa97 533
b8dc484b 534 histoFile->Close();
c687aa97 535
536 printf(" nb of reconstructible tracks: %d \n", nReconstructibleTracks);
537 printf(" nb of reconstructed tracks: %d \n", nReconstructedTracks);
538 printf(" nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
539
540}
541
542//------------------------------------------------------------------------------------
543Double_t langaufun(Double_t *x, Double_t *par) {
544
545 //Fit parameters:
546 //par[0]=Width (scale) parameter of Landau density
547 //par[1]=Most Probable (MP, location) parameter of Landau density
548 //par[2]=Total area (integral -inf to inf, normalization constant)
549 //par[3]=Width (sigma) of convoluted Gaussian function
550 //
551 //In the Landau distribution (represented by the CERNLIB approximation),
552 //the maximum is located at x=-0.22278298 with the location parameter=0.
553 //This shift is corrected within this function, so that the actual
554 //maximum is identical to the MP parameter.
555
556 // Numeric constants
557 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
558 Double_t mpshift = -0.22278298; // Landau maximum location
559
560 // Control constants
561 Double_t np = 100.0; // number of convolution steps
562 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
563
564 // Variables
565 Double_t xx;
566 Double_t mpc;
567 Double_t fland;
568 Double_t sum = 0.0;
569 Double_t xlow,xupp;
570 Double_t step;
571 Double_t i;
572
573
574 // MP shift correction
575 mpc = par[1] - mpshift * par[0];
576
577 // Range of convolution integral
578 xlow = x[0] - sc * par[3];
579 xupp = x[0] + sc * par[3];
580
581 step = (xupp-xlow) / np;
582
583 // Convolution integral of Landau and Gaussian by sum
584 for(i=1.0; i<=np/2; i++) {
585 xx = xlow + (i-.5) * step;
586 //change x -> -x because the tail of the Landau is at the left here...
587 fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
588 sum += fland * TMath::Gaus(x[0],xx,par[3]);
589
590 //change x -> -x because the tail of the Landau is at the left here...
591 xx = xupp - (i-.5) * step;
592 fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
593 sum += fland * TMath::Gaus(x[0],xx,par[3]);
594 }
595
596 return (par[2] * step * sum * invsq2pi / par[3]);
b8dc484b 597}
598