]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONTestAbso.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New class to make V2 clusters starting from digits or hits (fast simulation). Origin...
[u/mrichter/AliRoot.git] / MUON / MUONTestAbso.C
CommitLineData
b45fd22b 1void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0,Int_t testSingle=0)
b6340d07 2{
3 // Useful to test the absorber correction in the reconstruction procedure
4 // (mean energy loss + Branson correction)
5 // The AliMUONTrackParam class from the reconstruction is directly checked
6 // in this macro using the GEANT particle momentum downstream of the absorber.
7 // Histograms are saved in file MUONTestAbso.root.
8 // Use DrawTestAbso.C to display control histograms.
9
10
11 // Dynamically link some shared libs
12 if (gClassTable->GetID("AliRun") < 0) {
13 gROOT->LoadMacro("loadlibs.C");
14 loadlibs();
15 }
16
17 // Connect the Root Galice file containing Geometry, Kine and Hits
18
19 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
20
21
22 if (!file) {
23 printf("\n Creating galice.root \n");
24 file = new TFile("galice.root");
25 } else {
26 printf("\n galice.root found in file list");
27 }
28
29// Get AliRun object from file or create it if not on file
30 if (!gAlice) {
31 gAlice = (AliRun*)(file->Get("gAlice"));
32 if (gAlice) printf("AliRun object found on file\n");
33 if (!gAlice) {
34 printf("\n Create new gAlice object");
35 gAlice = new AliRun("gAlice","Alice test program");
36 }
37 }
38
39
40 // Create some histograms
41
b45fd22b 42 TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 1200, 0., 12.);
43// TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 100, 9., 10.5);
44 TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons theta < 3 deg", 100, -10., 10.);
45 TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons theta < 3 deg", 100, 0., 0.5);
b6340d07 46 //
b45fd22b 47 TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons theta > 3 deg", 100, -10., 10.);
48 TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons theta > 3 deg", 100, 0., 0.5);
49
50 //
51 TH1F *hRap = new TH1F("hRap"," Muon Rapidity gen.",20,2.5,4);
52 //
53
54 Int_t nBinProf=19;
55 Float_t xProfMin=5;
56 Float_t xProfMax=290;
57 char titleHist[50];
58 char numberHist[10];
59 TH1F *DeltaP1[100],*DeltaP2[100];
60 TH1F *DeltaAngleX1[100],*DeltaAngleX2[100];
61 TH1F *DeltaAngleY1[100],*DeltaAngleY2[100];
b6340d07 62
b45fd22b 63 TH1F *hP1 = new TH1F("hP1"," P theta < 3 deg ",nBinProf,xProfMin,xProfMax);
64 TProfile *hProfDeltaPvsP1 = new TProfile("hProfDeltaPvsP1"," delta P vs P theta < 3 deg ",nBinProf,xProfMin,xProfMax,-50,50,"s");
65 TH2F *h2DeltaPvsP1 = new TH2F("h2DeltaPvsP1"," delta P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
66 TH1F *hSigmaPvsP1 = new TH1F("hSigmaPvsP1"," deltaP/P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax);
67 for (Int_t iHist=0; iHist < nBinProf; iHist++){
68 sprintf(titleHist,"%s","hDelta1P");
69 sprintf(numberHist,"%d",iHist+1);
70 strcat(titleHist,numberHist);
71 DeltaP1[iHist] = new TH1F(titleHist," deltaP theta < 3 deg ",400,-50,50);
72
73 sprintf(titleHist,"%s","hDelta1AngleX");
74 sprintf(numberHist,"%d",iHist+1);
75 strcat(titleHist,numberHist);
76 DeltaAngleX1[iHist] = new TH1F(titleHist," delta Angle X theta < 3 deg ",400,-0.05,0.05);
77
78 sprintf(titleHist,"%s","hDelta1AngleY");
79 sprintf(numberHist,"%d",iHist+1);
80 strcat(titleHist,numberHist);
81 DeltaAngleY1[iHist] = new TH1F(titleHist," delta Angle Y theta < 3 deg ",400,-0.05,0.05);
82 }
83 TH1F *hSigGausPvsP1 = new TH1F("hSigGausPvsP1"," deltaP/P vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
84 TH1F *hSigGausAngleXvsP1 = new TH1F("hSigGausAngleXvsP1"," delta theta X vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
85 TH1F *hSigGausAngleYvsP1 = new TH1F("hSigGausAngleYvsP1"," delta theta Y vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
86
87 TH1F *hP2 = new TH1F("hP2"," P theta > 3 deg ",nBinProf,xProfMin,xProfMax);
88 TProfile *hProfDeltaPvsP2 = new TProfile("hProfDeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,-50,50,"s");
89 TH2F *h2DeltaPvsP2 = new TH2F("h2DeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
90 TH1F *hSigmaPvsP2 = new TH1F("hSigmaPvsP2"," deltaP/P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax);
91 for (Int_t iHist=0; iHist < nBinProf; iHist++){
92 sprintf(titleHist,"%s","hDelta2P");
93 sprintf(numberHist,"%d",iHist+1);
94 strcat(titleHist,numberHist);
95 DeltaP2[iHist] = new TH1F(titleHist," deltaP theta > 3 deg ",400,-50,50);
96
97 sprintf(titleHist,"%s","hDelta2AngleX");
98 sprintf(numberHist,"%d",iHist+1);
99 strcat(titleHist,numberHist);
100 DeltaAngleX2[iHist] = new TH1F(titleHist," delta Angle X theta > 3 deg ",400,-0.05,0.05);
101
102 sprintf(titleHist,"%s","hDelta2AngleY");
103 sprintf(numberHist,"%d",iHist+1);
104 strcat(titleHist,numberHist);
105 DeltaAngleY2[iHist] = new TH1F(titleHist," delta Angle Y theta > 3 deg ",400,-0.05,0.05);
106 }
107 TH1F *hSigGausPvsP2 = new TH1F("hSigGausPvsP2"," deltaP/P vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
108
109 TH1F *hSigGausAngleXvsP2 = new TH1F("hSigGausAngleXvsP2"," delta theta X vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
110 TH1F *hSigGausAngleYvsP2 = new TH1F("hSigGausAngleYvsP2"," delta theta Y vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
111
112 TProfile *hDeltaPxvsPhi = new TProfile("hDeltaPxvsPhi"," delta Px vs Phi",45,-180,180,-1,1,"s");
113 TProfile *hDeltaPyvsPhi = new TProfile("hDeltaPyvsPhi"," delta Py vs Phi",45,-180,180,-1,1,"s");
114 TProfile *hDeltaPhivsPz = new TProfile("hDeltaPhivsPz"," delta phi vs pZ",25,0,100,-4,4);
115
b6340d07 116 // Define constant parameters
117 Float_t massMuon = 0.105658389; // muon mass
118 Float_t massUpsilon = 9.46037; // Upsilon mass
b6340d07 119 Double_t zEndAbso = 503; // z position of the end of the absorber
b45fd22b 120 Double_t rLimit = zEndAbso * TMath::Tan(3*TMath::Pi()/180); // 3 deg. limit (different absorber composition)
b6340d07 121 Double_t printLevel = 0;
122
123 if (printLevel >= 1)
124 cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;
125
126 Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex
127
b45fd22b 128 for (Int_t nev=evNumber1; nev<= evNumber2; nev++) { // loop over events
129 if (printLevel >= 0 && nev%100 == 0) {
b6340d07 130 cout << " **** Event # " << nev <<endl;
131 }
132
133 Int_t nparticles = gAlice->GetEvent(nev);
134 if (printLevel >= 1) {
135 cout << " Total number of nparticles " << nparticles <<endl;
136 }
137 if (nev < evNumber1) continue;
138 if (nparticles <= 0) return;
139
b45fd22b 140 Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2],RapGen[2]; // parameters of generated particles at the vertex
141 ThetaGen[0] = ThetaGen[1] = 0;
142 RapGen[0] = RapGen[1] = 0;
143
b6340d07 144 for (int iPart = 0; iPart < nparticles ; iPart++) { // loop over particles
145
146 TParticle *particle;
147 particle = gAlice->Particle(iPart);
148
149 if ( particle->GetPdgCode() == kMuonPlus ) {
150 PxGen[0] = particle->Px();
151 PyGen[0] = particle->Py();
152 PzGen[0] = particle->Pz();
153 PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]);
b45fd22b 154 ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0]);
155 RapGen[0] = rapParticle(PxGen[0],PyGen[0],PzGen[0],massMuon);
b6340d07 156 if (printLevel >= 1) {
157 cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
158 }
159 }
160 if ( particle->GetPdgCode() == kMuonMinus ) {
161 PxGen[1] = particle->Px();
162 PyGen[1] = particle->Py();
163 PzGen[1] = particle->Pz();
164 PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]);
b45fd22b 165 ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1]);
166 RapGen[1] = rapParticle(PxGen[1],PyGen[1],PzGen[1],massMuon);
b6340d07 167 if (printLevel >= 1)
168 cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
169 }
b45fd22b 170 } // end loop over particles
b6340d07 171
172
173 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
174
175
176 TTree *TH = gAlice->TreeH();
177 Int_t ntracks =(Int_t) TH->GetEntries();
b45fd22b 178
179 Bool_t hitChamberMUPlus[10],hitChamberMUMinus[10];
180 Bool_t hitStationMUPlus[5],hitStationMUMinus[5];
181 for (Int_t iCh=0; iCh < 10; iCh++){
182 hitChamberMUPlus[iCh] = kFALSE;
183 hitChamberMUMinus[iCh] = kFALSE;
184 }
185 for (Int_t iSt=0; iSt < 5; iSt++){
186 hitStationMUPlus[iSt] = kFALSE;
187 hitStationMUMinus[iSt] = kFALSE;
188 }
189
190 Int_t nHitChamber=0;
191
b6340d07 192 for (Int_t i = 0; i < 2; i++){
193 pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
194 }
195
196 for (Int_t track=0; track<ntracks;track++) { // loop over tracks
197 gAlice->ResetHits();
198 Int_t nbytes += TH->GetEvent(track);
199 if (MUON) {
200 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
201 mHit;
202 mHit=(AliMUONHit*)MUON->NextHit()) // loop over GEANT muon hits
203 {
204 Int_t nch = mHit->Chamber(); // chamber number
205
b6340d07 206
207 Int_t ipart = mHit->Particle();
208 Double_t x = mHit->X(); // x-pos of hit in chamber #1
209 Double_t y = mHit->Y(); // y-pos of hit in chamber #1
210 Double_t z = mHit->Z(); // z-pos of hit in chamber #1
b45fd22b 211
212 Int_t iSt = (nch-1)/2;
213 if ((nch-1) < 10) {
214 if (ipart == kMuonPlus){
215 hitChamberMUPlus[nch-1] = kTRUE;
216 hitStationMUPlus[iSt] = kTRUE;
217 }
218 if (ipart == kMuonMinus){
219 hitChamberMUMinus[nch-1] = kTRUE;
220 hitStationMUMinus[iSt] = kTRUE;
221 }
222 nHitChamber++;
223 }
b6340d07 224
b45fd22b 225 if (nch != 1) continue;
b6340d07 226
227 if (ipart == kMuonPlus || ipart == kMuonMinus) {
228 Double_t px0=mHit->Px();
229 Double_t py0=mHit->Py();
230 Double_t pz0=mHit->Pz();
231 Double_t nonBendingSlope=px0/pz0;
232 Double_t bendingSlope=py0/pz0;
233
234 // create an AliMUONTrackParam object in chamber #1
235 AliMUONTrackParam *trackParam = new AliMUONTrackParam();
b45fd22b 236 Int_t sign = -TMath::Sign(1.0,ipart);
237 Double_t bendingMometum = sign* TMath::Sqrt(pz0*pz0+py0*py0); // Bug px -> py !!!
b6340d07 238 Double_t inverseBendingMomentum = 1/bendingMometum;
239
240 trackParam->SetBendingCoor(y);
241 trackParam->SetNonBendingCoor(x);
242 trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
243 trackParam->SetBendingSlope(bendingSlope);
244 trackParam->SetNonBendingSlope(nonBendingSlope);
245 trackParam->SetZ(z);
246
247 if (printLevel >= 1) {
248 cout <<" **** before extrap " <<endl;
249 cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
250 trackParam->Dump();
251 }
b45fd22b 252 if (!testSingle)
b6340d07 253 trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
254 if (printLevel >= 1) {
255 cout <<" **** after extrap " <<endl;
256 trackParam->Dump();
257 }
258
259 Int_t iMuon;
260 if (ipart == kMuonPlus) iMuon = 0;
261 else iMuon = 1;
262
263 // calculate track parameters extrapolated to vertex.
264 bendingSlope = trackParam->GetBendingSlope();
265 nonBendingSlope = trackParam->GetNonBendingSlope();
266 Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
267 pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
268 pX[iMuon] = pZ[iMuon] * nonBendingSlope;
269 pY[iMuon] = pZ[iMuon] * bendingSlope;
270 pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]);
b45fd22b 271 theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon]);
b6340d07 272 radius[iMuon] = TMath::Sqrt(x*x+y*y);
273
274 if (printLevel >= 1)
275 cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
276 delete trackParam;
277
278 } // if MuonPlus or MuonMinus
279 } // loop over MUON hits
280 } // if MUON
b45fd22b 281 } // loop over tracks
282
283
284 Bool_t goodMUPlus = kTRUE;
285 Bool_t goodMUMinus = kTRUE;
286 for (Int_t iCh=0; iCh < 1; iCh++) { // !!!! 10 -> 1
287 if (!hitChamberMUPlus[iCh]) {
288 goodMUPlus = kFALSE;
289 printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
290 }
291 if (!hitChamberMUMinus[iCh]) {
292 goodMUMinus = kFALSE;
293 printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
294 }
295 }
296
297// for (Int_t iSt=0; iSt < 5; iSt++) {
298// if (!hitStationMUPlus[iSt]) {
299// goodMUPlus = kFALSE;
300// // printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
301// }
302// if (!hitStationMUMinus[iSt]) {
303// goodMUMinus = kFALSE;
304// // printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
305// }
306// }
b6340d07 307
308 if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist
309 Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
310 hInvMassRes->Fill(massResonance);
b45fd22b 311 }
312
313
314 for (Int_t i = 0; i<2; i++){ // !!! 1 -> 2
315 if (i == 0 && !goodMUPlus) continue;
316 if (i == 1 && !goodMUMinus) continue;
317 if (pX[i] == 0) continue;
318
319 hRap->Fill(RapGen[i]);
320
321 Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
322 cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
323 Double_t deltaAngle = TMath::ACos(cosAngle)*180/TMath::Pi();
324 Double_t deltaPNorm = 0;
325 if (testSingle) {
326 deltaPNorm = (PTotGen[i]-pTot[i])* TMath::Cos(ThetaGen[i]);
327 } else {
328 deltaPNorm = PTotGen[i]-pTot[i];
329 }
330
331 Float_t thetaGX = TMath::ATan2(PxGen[i],PzGen[i]);
332 Float_t thetaRX = TMath::ATan2(pX[i],pZ[i]);
333 Float_t deltaAngleX = thetaRX - thetaGX;
334
335 Float_t thetaGY = TMath::ATan2(PyGen[i],PzGen[i]);
336 Float_t thetaRY = TMath::ATan2(pY[i],pZ[i]);
337 Float_t deltaAngleY = thetaRY - thetaGY;
338
339 Float_t phiG = TMath::ATan2(PyGen[i],PxGen[i])*180/TMath::Pi();
340 Float_t phiR = TMath::ATan2(pY[i],pX[i])*180/TMath::Pi();
341
342 hDeltaPxvsPhi->Fill(phiG,PxGen[i]-pX[i]);
343 hDeltaPyvsPhi->Fill(phiG,PyGen[i]-pY[i]);
344 hDeltaPhivsPz->Fill(PzGen[i],phiR-phiG);
345
346 Int_t iHist=0;
347 if (ThetaGen[i] < (3*TMath::Pi()/180)) {
348 hDeltaP1->Fill(pTot[i]-PTotGen[i]);
349 hDeltaAngle1->Fill(deltaAngle);
350 hP1->Fill(PTotGen[i]);
351 hProfDeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
352 h2DeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
353 iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
354 if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) {
355 DeltaP1[iHist]->Fill(deltaPNorm);
356 DeltaAngleX1[iHist]->Fill(deltaAngleX);
357 DeltaAngleY1[iHist]->Fill(deltaAngleY);
358 }
359 } else {
360 hDeltaP2->Fill(pTot[i]-PTotGen[i]);
361 hDeltaAngle2->Fill(deltaAngle);
362 hP2->Fill(PTotGen[i]);
363 hProfDeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
364 h2DeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
365 iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
366 if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) {
367 DeltaP2[iHist]->Fill(deltaPNorm);
368 DeltaAngleX2[iHist]->Fill(deltaAngleX);
369 DeltaAngleY2[iHist]->Fill(deltaAngleY);
b6340d07 370 }
371 }
372 }
373 } // loop over event
b45fd22b 374
375
376 Float_t sigmaP = 0;
377 Float_t xBin = 0;
378 for (Int_t iBin = 1; iBin <= nBinProf; iBin++){
379 sigmaP = hProfDeltaPvsP1->GetBinError(iBin);
380 xBin = hProfDeltaPvsP1->GetBinCenter(iBin);
381 hSigmaPvsP1->SetBinContent(iBin,sigmaP/xBin);
382 sigmaP = hProfDeltaPvsP2->GetBinError(iBin);
383 xBin = hProfDeltaPvsP2->GetBinCenter(iBin);
384 hSigmaPvsP2->SetBinContent(iBin,sigmaP/xBin);
385 }
386
387 TF1 *g1= new TF1("g1","gaus",-25,25) ;
388 TF1 *g2= new TF1("g2","gaus",-25,25) ;
389 Float_t sigmaPGaus;
390 Float_t xBinGaus;
391 for (Int_t iHist=0; iHist < nBinProf; iHist++){
392 DeltaP1[iHist]->Fit("g1","RQ");
393 sigmaPGaus = g1->GetParameter(2);
394 printf(" ** 1 ** iHist= %d , sigmaPGaus = %f \n",iHist,sigmaPGaus);
395 xBinGaus = hSigGausPvsP1->GetBinCenter(iHist+1);
396 hSigGausPvsP1->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
397
398 DeltaP2[iHist]->Fit("g2","RQ");
399 sigmaPGaus = g2->GetParameter(2);
400 printf(" ** 2 ** iHist= %d , sigmaGaus = %f \n",iHist,sigmaPGaus);
401 xBinGaus = hSigGausPvsP2->GetBinCenter(iHist+1);
402 hSigGausPvsP2->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
403 }
b6340d07 404
b45fd22b 405
406 // Angles
407 TF1 *g3= new TF1("g3","gaus") ;
408 TF1 *g4= new TF1("g4","gaus") ;
409 TF1 *g5= new TF1("g5","gaus") ;
410 TF1 *g6= new TF1("g6","gaus") ;
411 Float_t sigmaAngleGaus,limitGaus,errorSigma;
412 Float_t nSigFit = 3;
413 for (Int_t iHist=0; iHist < nBinProf; iHist++){
414 limitGaus = nSigFit * (DeltaAngleX1[iHist]->GetRMS());
415 g3->SetRange(-limitGaus,limitGaus);
416 DeltaAngleX1[iHist]->Fit("g3","RQ");
417 sigmaAngleGaus = g3->GetParameter(2);
418 errorSigma = g3->GetParError(2);
419 printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
420 hSigGausAngleXvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
421 hSigGausAngleXvsP1->SetBinError(iHist+1,errorSigma);
422
423 limitGaus = nSigFit * (DeltaAngleY1[iHist]->GetRMS());
424 g4->SetRange(-limitGaus,limitGaus);
425 DeltaAngleY1[iHist]->Fit("g4","RQ");
426 sigmaAngleGaus = g4->GetParameter(2);
427 errorSigma = g4->GetParError(2);
428 printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
429 hSigGausAngleYvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
430 hSigGausAngleYvsP1->SetBinError(iHist+1,errorSigma);
431
432 limitGaus = nSigFit * (DeltaAngleX2[iHist]->GetRMS());
433 g5->SetRange(-limitGaus,limitGaus);
434 DeltaAngleX2[iHist]->Fit("g5","RQ");
435 sigmaAngleGaus = g5->GetParameter(2);
436 errorSigma = g5->GetParError(2);
437 printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
438 hSigGausAngleXvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
439 hSigGausAngleXvsP2->SetBinError(iHist+1,errorSigma);
440
441 limitGaus = nSigFit * (DeltaAngleY2[iHist]->GetRMS());
442 g6->SetRange(-limitGaus,limitGaus);
443 DeltaAngleY2[iHist]->Fit("g6","RQ");
444 sigmaAngleGaus = g6->GetParameter(2);
445 errorSigma = g6->GetParError(2);
446 printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
447 hSigGausAngleYvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
448 hSigGausAngleYvsP2->SetBinError(iHist+1,errorSigma);
449
450 }
b6340d07 451 // save histograms in MUONTestAbso.root
452 TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
b45fd22b 453
b6340d07 454 hInvMassRes->Write();
b45fd22b 455 hRap->Write();
456
b6340d07 457 hDeltaP1->Write();
458 hDeltaAngle1->Write();
b45fd22b 459 hP1->Write();
460 hProfDeltaPvsP1->Write();
461 h2DeltaPvsP1->Write();
462 hSigmaPvsP1->Write();
463 hSigGausPvsP1->Write();
464 hSigGausAngleXvsP1->Write();
465 hSigGausAngleYvsP1->Write();
466
b6340d07 467 hDeltaP2->Write();
468 hDeltaAngle2->Write();
b45fd22b 469 hP2->Write();
470 hProfDeltaPvsP2->Write();
471 h2DeltaPvsP2->Write();
472 hSigmaPvsP2->Write();
473 hSigGausPvsP2->Write();
474 hSigGausAngleXvsP2->Write();
475 hSigGausAngleYvsP2->Write();
476
477 for (Int_t iHist=0; iHist < nBinProf; iHist++){
478 DeltaP1[iHist]->Write();
479 DeltaP2[iHist]->Write();
480 DeltaAngleX1[iHist]->Write();
481 DeltaAngleY1[iHist]->Write();
482 DeltaAngleX2[iHist]->Write();
483 DeltaAngleY2[iHist]->Write();
484 }
485
486 hDeltaPxvsPhi->Write();
487 hDeltaPyvsPhi->Write();
488 hDeltaPhivsPz->Write();
489
b6340d07 490 histoFile->Close();
491
492}
493
494
495Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
496{
497 // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass
498
499 Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
500 Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
501 return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
502}
b45fd22b 503
504Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
505{
506 // return rapidity for particle
507
508 // Particle energy
509 Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
510 // Rapidity
511 Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
512 return Rapidity;
513}