]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MFT/AliMuonForwardTrackAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
use new names for the PWG subdirectories
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackAnalysis.cxx
CommitLineData
e806e863 1//====================================================================================================================================================
2//
3// Class for the analysis of the ALICE muon forward tracks (MUON + MFT)
4//
5// Contact author: antonio.uras@cern.ch
6//
7//====================================================================================================================================================
8
9#include "TObject.h"
10#include "TClonesArray.h"
11#include "AliMuonForwardTrack.h"
12#include "AliMuonForwardTrackPair.h"
13#include "TMatrixD.h"
14#include "TTree.h"
15#include "TH1D.h"
16#include "AliLog.h"
17#include "TFile.h"
18#include "TParticle.h"
19#include "AliMUONTrackParam.h"
20#include "AliMUONTrackExtrap.h"
21#include "TGeoManager.h"
22#include "TRandom.h"
23#include "TLorentzVector.h"
24#include "TDatabasePDG.h"
25#include "AliMuonForwardTrackAnalysis.h"
26
27ClassImp(AliMuonForwardTrackAnalysis)
28
29//====================================================================================================================================================
30
31AliMuonForwardTrackAnalysis::AliMuonForwardTrackAnalysis():
32 TObject(),
33 fInputDir(0),
34 fOutputDir(0),
35 fInputTree(0),
36 fMuonForwardTracks(0),
37 fMuonForwardTrackPairs(0),
38 fMFTTrack(0),
39 fMFTTrackPair(0),
40 fMCRefTrack(0),
41 fEv(0),
42 fFirstEvent(-1),
43 fLastEvent(-1),
44 fNTracksOfEvent(0),
45 fNTracksAnalyzedOfEvent(0),
46 fNTracksAnalyzed(0),
47 fNPairsOfEvent(0),
48 fNPairsAnalyzedOfEvent(0),
49 fHistOffsetSingleMuonsX(0x0),
50 fHistOffsetSingleMuonsY(0x0),
51 fHistOffsetSingleMuons(0x0),
52 fHistWOffsetSingleMuons(0x0),
53 fHistErrorSingleMuonsX(0x0),
54 fHistErrorSingleMuonsY(0x0),
55 fHistOffsetSingleMuonsX_vsPtRapidity(0x0),
56 fHistOffsetSingleMuonsY_vsPtRapidity(0x0),
57 fHistSingleMuonsPtRapidity(0x0),
58 fHistWOffsetMuonPairs(0x0),
59 fHistMassMuonPairs(0x0),
60 fHistMassMuonPairsWithoutMFT(0x0),
61 fHistMassMuonPairsMC(0x0),
62 fHistRapidityPtMuonPairsMC(0x0),
63 fNMassBins(1000),
64 fMassMin(0),
65 fMassMax(10),
66 fSingleMuonAnalysis(1),
67 fMuonPairAnalysis(1),
68 fMatchTrigger(0),
69 fOption(0),
70 fXVertResMC(150.e-4),
71 fYVertResMC(150.e-4),
72 fZVertResMC(100.e-4),
73 fMaxNWrongClustersMC(999),
74 fPtMinSingleMuons(0)
75{
76
77 // default constructor
78
79 for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
80 for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
81 fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] = 0x0;
82 fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] = 0x0;
83 }
84 }
85
86}
87
88//====================================================================================================================================================
89
90Bool_t AliMuonForwardTrackAnalysis::Init(Char_t *inputFileName) {
91
92 BookHistos();
93
94 TFile *inputFile = new TFile(Form("%s/%s",fInputDir.Data(),inputFileName));
95 if (!inputFile || !inputFile->IsOpen()) {
96 AliError(Form("Error opening file %s", inputFileName));
97 return kFALSE;
98 }
99 fInputTree = (TTree*) inputFile->Get("AliMuonForwardTracks");
100 if (!fInputTree) {
101 AliError("Error reading input tree");
102 return kFALSE;
103 }
104
105 if (fFirstEvent<0 || fLastEvent<0 || fFirstEvent>fLastEvent || fFirstEvent>=fInputTree->GetEntries()) {
106 fFirstEvent = 0;
107 fLastEvent = fInputTree->GetEntries()-1;
108 }
109 else {
110 fLastEvent = TMath::Min(fLastEvent, Int_t(fInputTree->GetEntries()-1));
111 }
112
113 AliInfo(Form("Analysing events %d to %d", fFirstEvent, fLastEvent));
114
115 fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
116 fInputTree->SetBranchAddress("tracks", &fMuonForwardTracks);
117
118 TGeoManager::Import(Form("%s/geometry.root",fInputDir.Data()));
119
120 AliMUONTrackExtrap::SetField();
121
122 fMuonForwardTrackPairs = new TClonesArray("AliMuonForwardTrackPair");
123
124 return kTRUE;
125
126}
127
128//====================================================================================================================================================
129
130Bool_t AliMuonForwardTrackAnalysis::LoadNextEvent() {
131
132 if (fEv>fLastEvent) return kFALSE;
133 if (fEv<fFirstEvent) { fEv++; return kTRUE; }
134 fMuonForwardTracks -> Clear();
135 fInputTree->GetEvent(fEv);
136 AliInfo(Form("**** analyzing event # %4d (%3d tracks) ****", fEv, fMuonForwardTracks->GetEntries()));
137
138 if (fSingleMuonAnalysis) {
139 fNTracksAnalyzedOfEvent = 0;
140 fNTracksOfEvent = fMuonForwardTracks->GetEntries();
141 while (AnalyzeSingleMuon()) continue;
142 }
143
144 if (fMuonPairAnalysis) {
145 if (fMuonForwardTrackPairs) {
146 fMuonForwardTrackPairs->Clear();
147 }
148 BuildMuonPairs();
149 fNPairsAnalyzedOfEvent = 0;
150 fNPairsOfEvent = fMuonForwardTrackPairs->GetEntries();
151 while (AnalyzeMuonPair()) continue;
152 }
153
154 fEv++;
155
156 return kTRUE;
157
158}
159
160//====================================================================================================================================================
161
162Bool_t AliMuonForwardTrackAnalysis::AnalyzeSingleMuon() {
163
164 if (fNTracksAnalyzedOfEvent>=fNTracksOfEvent) return kFALSE;
165
166 fMFTTrack = (AliMuonForwardTrack*) fMuonForwardTracks->At(fNTracksAnalyzedOfEvent);
167 fNTracksAnalyzedOfEvent++;
168 if (fMatchTrigger && !fMFTTrack->GetMatchTrigger()) return kTRUE;
169 fMCRefTrack = fMFTTrack->GetMCTrackRef();
170
171 if (!fMCRefTrack) return kTRUE;
172 if (fMFTTrack->GetNWrongClustersMC()>fMaxNWrongClustersMC) return kTRUE;
173
174 Double_t xOrig=gRandom->Gaus(0., fXVertResMC);
175 Double_t yOrig=gRandom->Gaus(0., fXVertResMC);
176 Double_t zOrig=gRandom->Gaus(0., fZVertResMC);
177
178 AliMUONTrackParam *param = fMFTTrack->GetTrackParamAtMFTCluster(0);
179 AliMUONTrackExtrap::ExtrapToZCov(param, zOrig);
180
181 TLorentzVector pMu;
182 Double_t mMu = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
183 Double_t energy = TMath::Sqrt(param->P()*param->P() + mMu*mMu);
184 pMu.SetPxPyPzE(param->Px(), param->Py(), param->Pz(), energy);
185
186 if (fMFTTrack->Pt()<fPtMinSingleMuons) return kTRUE;
187
188 TMatrixD cov(5,5);
189 cov = param->GetCovariances();
190
191 fHistErrorSingleMuonsX -> Fill(1.e4*TMath::Sqrt(cov(0,0)));
192 fHistErrorSingleMuonsY -> Fill(1.e4*TMath::Sqrt(cov(2,2)));
193
194 Double_t dX = fMFTTrack->GetOffsetX(xOrig, zOrig);
195 Double_t dY = fMFTTrack->GetOffsetY(yOrig, zOrig);
196
197 Double_t offset = fMFTTrack->GetOffset(xOrig, yOrig, zOrig);
198 Double_t weightedOffset = fMFTTrack->GetWeightedOffset(xOrig, yOrig, zOrig);
199
200 // AliDebug(2, Form("pdg code = %d\n", fMCRefTrack->GetPdgCode()));
201
202 fHistOffsetSingleMuonsX -> Fill(1.e4*dX);
203 fHistOffsetSingleMuonsY -> Fill(1.e4*dY);
204 Int_t rapBin = fHistOffsetSingleMuonsX_vsPtRapidity->GetXaxis()->FindBin(pMu.Rapidity());
205 Int_t ptBin = fHistOffsetSingleMuonsX_vsPtRapidity->GetYaxis()->FindBin(pMu.Pt());
206 // AliDebug(1, Form("pt = %f (%d), rap = %f (%d)\n", pMu.Pt(), pMu.Rapidity(), ptBin, rapBin));
207 if (0<rapBin && rapBin<=fNRapBinsOffsetSingleMuons && 0<ptBin && ptBin<=fNPtBinsOffsetSingleMuons) {
208 fHistOffsetSingleMuonsX_tmp[rapBin-1][ptBin-1]->Fill(1.e4*dX);
209 fHistOffsetSingleMuonsY_tmp[rapBin-1][ptBin-1]->Fill(1.e4*dY);
210 }
211 fHistSingleMuonsPtRapidity -> Fill(pMu.Rapidity(), pMu.Pt());
212 fHistOffsetSingleMuons -> Fill(1.e4*offset);
213 fHistWOffsetSingleMuons -> Fill(weightedOffset);
214
215 fNTracksAnalyzed++;
216
217 return kTRUE;
218
219}
220
221//====================================================================================================================================================
222
223Bool_t AliMuonForwardTrackAnalysis::AnalyzeMuonPair() {
224
225 if (fNPairsAnalyzedOfEvent>=fNPairsOfEvent) return kFALSE;
226
227 fMFTTrackPair = (AliMuonForwardTrackPair*) fMuonForwardTrackPairs->At(fNPairsAnalyzedOfEvent);
228
229 if (fOption==kResonanceOnly && !fMFTTrackPair->IsResonance()) {
230 fNPairsAnalyzedOfEvent++;
231 return kTRUE;
232 }
233
234 Double_t xOrig=gRandom->Gaus(0., fXVertResMC);
235 Double_t yOrig=gRandom->Gaus(0., fXVertResMC);
236 Double_t zOrig=gRandom->Gaus(0., fZVertResMC);
237 AliDebug(1, Form("origin = (%f, %f, %f)", xOrig, yOrig, zOrig));
238
239 fHistMassMuonPairs -> Fill(fMFTTrackPair->GetMass(zOrig));
240 fHistWOffsetMuonPairs -> Fill(fMFTTrackPair->GetWeightedOffset(xOrig, yOrig, zOrig));
241 fHistMassMuonPairsWithoutMFT -> Fill(fMFTTrackPair->GetMassWithoutMFT(xOrig, yOrig, zOrig));
242 fHistMassMuonPairsMC -> Fill(fMFTTrackPair->GetMassMC());
243 fHistRapidityPtMuonPairsMC -> Fill(fMFTTrackPair->GetRapidityMC(), fMFTTrackPair->GetPtMC());
244
245 AliDebug(1, Form("mass = %f MC = %f", fMFTTrackPair->GetMass(zOrig), fMFTTrackPair->GetMassMC()));
246
247 fNPairsAnalyzedOfEvent++;
248
249 return kTRUE;
250
251}
252
253//====================================================================================================================================================
254
255void AliMuonForwardTrackAnalysis::BuildMuonPairs() {
256
257 for (Int_t iTrack=0; iTrack<fMuonForwardTracks->GetEntries(); iTrack++) {
258 for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
259
260 AliMuonForwardTrack *track0 = (AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack);
261 AliMuonForwardTrack *track1 = (AliMuonForwardTrack*) fMuonForwardTracks->At(jTrack);
262
263 if (fMatchTrigger) if (!track0->GetMatchTrigger() || !track1->GetMatchTrigger()) continue;
264 if (!track0->GetMCTrackRef() || !track1->GetMCTrackRef()) continue;
265 if (track0->GetNWrongClustersMC()>fMaxNWrongClustersMC || track1->GetNWrongClustersMC()>fMaxNWrongClustersMC) continue;
266 if (track0->Pt()<fPtMinSingleMuons || track1->Pt()<fPtMinSingleMuons) continue;
267
268 AliMuonForwardTrackPair *trackPair = new AliMuonForwardTrackPair(track0, track1);
269 if (fOption==kResonanceOnly && !trackPair->IsResonance()) {
270 delete trackPair;
271 continue;
272 }
273 new ((*fMuonForwardTrackPairs)[fMuonForwardTrackPairs->GetEntries()]) AliMuonForwardTrackPair(*trackPair);
274
275 }
276 }
277
278}
279
280//====================================================================================================================================================
281
282void AliMuonForwardTrackAnalysis::Terminate(Char_t *outputFileName) {
283
284 for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
285 for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
286 Int_t binMin_x = fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->FindBin(-3*fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
287 Int_t binMax_x = fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->FindBin(+3*fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
288 for (Int_t bin=1; bin<=fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetNbinsX(); bin++) {
289 if (bin<binMin_x || bin>binMax_x) fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->SetBinContent(bin, 0);
290 }
291 Int_t binMin_y = fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->FindBin(-3*fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
292 Int_t binMax_y = fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->FindBin(+3*fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
293 for (Int_t bin=1; bin<=fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetNbinsX(); bin++) {
294 if (bin<binMin_y || bin>binMax_y) fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->SetBinContent(bin, 0);
295 }
296 fHistOffsetSingleMuonsX_vsPtRapidity->SetBinContent(rapBin+1, ptBin+1, fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMS());
297 fHistOffsetSingleMuonsX_vsPtRapidity->SetBinError(rapBin+1, ptBin+1, fHistOffsetSingleMuonsX_tmp[rapBin][ptBin]->GetRMSError());
298 fHistOffsetSingleMuonsY_vsPtRapidity->SetBinContent(rapBin+1, ptBin+1, fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMS());
299 fHistOffsetSingleMuonsY_vsPtRapidity->SetBinError(rapBin+1, ptBin+1, fHistOffsetSingleMuonsY_tmp[rapBin][ptBin]->GetRMSError());
300 }
301 }
302
303 TFile *fileOut = new TFile(Form("%s/%s",fOutputDir.Data(),outputFileName), "recreate");
304
305 printf("Writing output objects to file %s\n", fileOut->GetName());
306
307 fHistOffsetSingleMuonsX -> Write();
308 fHistOffsetSingleMuonsY -> Write();
309 fHistOffsetSingleMuons -> Write();
310 fHistWOffsetSingleMuons -> Write();
311 fHistErrorSingleMuonsX -> Write();
312 fHistErrorSingleMuonsY -> Write();
313
314 fHistOffsetSingleMuonsX_vsPtRapidity -> Write();
315 fHistOffsetSingleMuonsY_vsPtRapidity -> Write();
316
317// for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
318// for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
319// fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] -> Write();
320// fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] -> Write();
321// }
322// }
323
324 fHistSingleMuonsPtRapidity -> Write();
325
326 fHistWOffsetMuonPairs -> Write();
327 fHistMassMuonPairs -> Write();
328 fHistMassMuonPairsWithoutMFT -> Write();
329 fHistMassMuonPairsMC -> Write();
330 fHistRapidityPtMuonPairsMC -> Write();
331
332 fileOut -> Close();
333
334}
335
336//====================================================================================================================================================
337
338void AliMuonForwardTrackAnalysis::BookHistos() {
339
340 fHistOffsetSingleMuonsX = new TH1D("fHistOffsetSingleMuonsX", "Offset for single muons along X", 200, -1000, 1000);
341 fHistOffsetSingleMuonsY = new TH1D("fHistOffsetSingleMuonsY", "Offset for single muons along Y", 200, -1000, 1000);
342 fHistErrorSingleMuonsX = new TH1D("fHistErrorSingleMuonsX", "Coordinate Error for single muons along X", 200, 0, 1000);
343 fHistErrorSingleMuonsY = new TH1D("fHistErrorSingleMuonsY", "Coordinate Error for single muons along Y", 200, 0, 1000);
344 fHistOffsetSingleMuons = new TH1D("fHistOffsetSingleMuons", "Offset for single muons", 100, 0, 2000);
345 fHistWOffsetSingleMuons = new TH1D("fHistWOffsetSingleMuons", "Weighted Offset for single muons", 300, 0, 15);
346
347 fHistOffsetSingleMuonsX_vsPtRapidity = new TH2D("fHistOffsetSingleMuonsX_vsPtRapidity", "Offset for single muons along X",
348 10, -4, -2.5, 10, 0.5, 5.5);
349 fHistOffsetSingleMuonsY_vsPtRapidity = new TH2D("fHistOffsetSingleMuonsY_vsPtRapidity", "Offset for single muons along Y",
350 10, -4, -2.5, 10, 0.5, 5.5);
351
352 for (Int_t rapBin=0; rapBin<fNRapBinsOffsetSingleMuons; rapBin++) {
353 for (Int_t ptBin=0; ptBin<fNPtBinsOffsetSingleMuons; ptBin++) {
354 fHistOffsetSingleMuonsX_tmp[rapBin][ptBin] = new TH1D(Form("fHistOffsetSingleMuonsX_tmp_%02d_%02d",rapBin,ptBin), "", 1000, -1000, 1000);
355 fHistOffsetSingleMuonsY_tmp[rapBin][ptBin] = new TH1D(Form("fHistOffsetSingleMuonsY_tmp_%02d_%02d",rapBin,ptBin), "", 1000, -1000, 1000);
356 }
357 }
358
359 fHistSingleMuonsPtRapidity = new TH2D("fHistSingleMuonsPtRapidity", "Phase Space for single muons", 10, -4, -2.5, 10, 0.5, 5.5);
360
361 fHistOffsetSingleMuonsX -> SetXTitle("Offset(X) [#mum]");
362 fHistOffsetSingleMuonsY -> SetXTitle("Offset(Y) [#mum]");
363 fHistErrorSingleMuonsX -> SetXTitle("Err. on track position at z_{vtx} (X) [#mum]");
364 fHistErrorSingleMuonsY -> SetXTitle("Err. on track position at z_{vtx} (Y) [#mum]");
365 fHistOffsetSingleMuons -> SetXTitle("Offset [#mum]");
366 fHistWOffsetSingleMuons -> SetXTitle("Weighted Offset");
367
368 fHistOffsetSingleMuonsX_vsPtRapidity -> SetXTitle("y^{#mu}");
369 fHistOffsetSingleMuonsY_vsPtRapidity -> SetXTitle("y^{#mu}");
370 fHistOffsetSingleMuonsX_vsPtRapidity -> SetYTitle("p_{T}^{#mu} [GeV/c]");
371 fHistOffsetSingleMuonsY_vsPtRapidity -> SetYTitle("p_{T}^{#mu} [GeV/c]");
372
373 fHistSingleMuonsPtRapidity -> SetXTitle("y^{#mu}");
374 fHistSingleMuonsPtRapidity -> SetYTitle("p_{T}^{#mu} [GeV/c]");
375
376 fHistOffsetSingleMuonsX -> Sumw2();
377 fHistOffsetSingleMuonsY -> Sumw2();
378 fHistErrorSingleMuonsX -> Sumw2();
379 fHistErrorSingleMuonsY -> Sumw2();
380 fHistOffsetSingleMuons -> Sumw2();
381 fHistWOffsetSingleMuons -> Sumw2();
382
383 fHistOffsetSingleMuonsX_vsPtRapidity -> Sumw2();
384 fHistOffsetSingleMuonsY_vsPtRapidity -> Sumw2();
385 fHistSingleMuonsPtRapidity -> Sumw2();
386
387 //--------------------------------------------
388
389 fHistWOffsetMuonPairs = new TH1D("fHistWOffsetMuonPairs", "Weighted Offset for Muon Pairs", 300, 0, 60);
390 fHistMassMuonPairs = new TH1D("fHistMassMuonPairs", "Dimuon Mass (MUON+MFT)", fNMassBins, fMassMin, fMassMax);
391 fHistMassMuonPairsWithoutMFT = new TH1D("fHistMassMuonPairsWithoutMFT", "Dimuon Mass (MUON only)", fNMassBins, fMassMin, fMassMax);
392 fHistMassMuonPairsMC = new TH1D("fHistMassMuonPairsMC", "Dimuon Mass (MC)", fNMassBins, fMassMin, fMassMax);
393 fHistRapidityPtMuonPairsMC = new TH2D("fHistRapidityPtMuonPairsMC", "Dimuon Phase Space (MC)", 100, -4.5, -2., 100, 0., 10.);
394
395 fHistWOffsetMuonPairs -> SetXTitle("Weighted Offset");
396 fHistMassMuonPairs -> SetXTitle("Mass [GeV/c^{2}]");
397 fHistMassMuonPairsWithoutMFT -> SetXTitle("Mass [GeV/c^{2}]");
398 fHistMassMuonPairsMC -> SetXTitle("Mass [GeV/c^{2}]");
399 fHistRapidityPtMuonPairsMC -> SetXTitle("Rapidity");
400 fHistRapidityPtMuonPairsMC -> SetYTitle("p_{T} [GeV/c]");
401
402 fHistWOffsetMuonPairs -> Sumw2();
403 fHistMassMuonPairs -> Sumw2();
404 fHistMassMuonPairsWithoutMFT -> Sumw2();
405 fHistMassMuonPairsMC -> Sumw2();
406 fHistRapidityPtMuonPairsMC -> Sumw2();
407
408}
409
410//====================================================================================================================================================
411