]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONResoEffChamber.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
READMErec.txt updated.
[u/mrichter/AliRoot.git] / MUON / MUONResoEffChamber.C
CommitLineData
8ca62886 1/**************************************************************************
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 **************************************************************************/
15
16/* $Id$ */
17
e54bf126 18/// \ingroup macros
19/// \file MUONResoEffChamber.C
20/// \brief Macro to calculate the resolution and the efficiency of chamber(s)
21///
22/// Macro to calculate the resolution and the efficiency of chamber(s) chosen in function
23/// of Phi (angle on the chamber between the X axis and the straight line created by the
24/// center of the chamber and the impact of particles on this chamber, the azimuthal angle)
25/// and Theta (the polar angle), or in function of ThetaI (angle of incidence of particles
26/// on the chamber)
27///
28/// \author Nicolas Le Bris, Subatech
8ca62886 29
bd022fc7 30#if !defined(__CINT__) || defined(__MAKECINT__)
31
32// ROOT includes
33#include "TBranch.h"
34#include "TClonesArray.h"
35#include "TTree.h"
36#include "TNtuple.h"
37#include "TParticle.h"
38#include "TFile.h"
39
40#include "TH1.h"
41#include "TH1F.h"
42#include "TH2F.h"
43#include "TF1.h"
44#include "TMath.h"
45
46#include "TCanvas.h"
47#include "TGraph.h"
48#include "TGraphErrors.h"
49#include "TGraph2D.h"
50#include "TStyle.h"
51#include "TFitter.h"
52#include "TRandom.h"
53
54// STEER includes
55#include "AliRun.h"
56#include "AliRunLoader.h"
57#include "AliHeader.h"
58#include "AliLoader.h"
59#include "AliTracker.h"
60#include "AliStack.h"
f7a1cc68 61#include "AliMagF.h"
bd022fc7 62
63
64// MUON includes
65#include "AliMUON.h"
66#include "AliMUONData.h"
67#include "AliMUONConstants.h"
68
69#include "AliMUONHit.h"
70#include "AliMUONHitForRec.h"
bd022fc7 71#include "AliMUONTrack.h"
72#include "AliMUONTrackParam.h"
73#include "AliMUONTrackExtrap.h"
74
75#include "AliMpVSegmentation.h"
bd022fc7 76#include "AliMpDEManager.h"
77#endif
78
79
bd022fc7 80
81 const Double_t kInvPi = 1./3.14159265;
82
83
84//Chamber number:
85 Int_t chamberNbr;
86//Number of events:
87 Int_t nEvents, iEvent;
88//Number of tracks:
89 Int_t nTracks, iTrack;
90//Number of hits:
91 Int_t nHits,iHit;
92//Chamber(s) chosen
93 Int_t firstChamber, lastChamber;
94
95
96 AliMUONTrack * track ;
97 AliMUONHitForRec * hit = 0;
98 AliMUONTrackParam * trackParam = 0;
99
100 TClonesArray * tracks ;
101 TClonesArray * trackParams;
102 TClonesArray * hits ;
103
104
105
106
107
108
109/*****************************************************************************************************************/
110/*****************************************************************************************************************/
111 /*EFFICIENCY*/
112
113void efficiency( Int_t event2Check=0, char * filename="galice.root" )
114{
115 cout<<"\nChamber(s) chosen;\nFirst chamber:";
116 cin>>firstChamber;
117 cout<<"Last chamber:";
118 cin>>lastChamber;
119 cout<<"\n\n";
120
121//Creating Run Loader and openning file containing Hits
122//--------------------------------------------------------------------------------------------------------------
123 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
124 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
125
126
127//Getting MUONLoader
128//--------------------------------------------------------------------------------------------------------------
129 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
130 MUONLoader -> LoadTracks("READ");
131 MUONLoader -> LoadRecPoints("READ");
132
133
134//Creating a MUON data container
135//--------------------------------------------------------------------------------------------------------------
136 AliMUONData muondata(MUONLoader,"MUON","MUON");
137
138
139 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
140
141 //Loop on events
142 Int_t trackNb = 0;
143 Int_t chamber[10] = {0};
144 Int_t detEltNew, detElt;
145 Int_t detEltOld = 0;
146
147 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
148 {
149 if ( event2Check!=0 )
150 iEvent=event2Check;
151 printf("\r>>> Event %d ",iEvent);
152 RunLoader -> GetEvent(iEvent);
153 //Addressing
154 muondata.SetTreeAddress("RT");
155 muondata.GetRecTracks();
156 tracks = muondata.RecTracks();
157
158 //Loop on track
159 nTracks = (Int_t) tracks -> GetEntriesFast();
160 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
161 {
162 track = (AliMUONTrack*) tracks -> At(iTrack);
163 hits = track -> GetHitForRecAtHit();
164 detEltOld = 0;
165 //Loop on hit
166 nHits = (Int_t) hits -> GetEntriesFast();
167
168 for ( iHit = 0; iHit < nHits; ++iHit )
169 {
170 hit = (AliMUONHitForRec*) hits -> At(iHit);
171 chamberNbr = hit -> GetChamberNumber();
172 detElt = hit -> GetDetElemId();
173 detEltNew = int(detElt/100);
174 if( chamberNbr >= firstChamber-1 ) {
175 if( chamberNbr < lastChamber ) {
176 if( detEltNew == detEltOld )
177 continue;
178 else {
179 chamber[chamberNbr] = chamber[chamberNbr] + 1;
180 detEltOld = detEltNew;
181 }
182 }
183 }
184 }
185 //End loop on hit
186
187 }
188 //End loop on track
189 muondata.ResetRecTracks();
190 if (event2Check != 0)
191 iEvent = nEvents;
192 trackNb = trackNb + nTracks;
193 }
194 //End loop on event
195//--------------------------------------------------------------------------------------------------------------
196
197 cout<<"\n\n\n";
198 for (Int_t i = firstChamber-1; i < lastChamber; ++i )
199 {
200 printf ("\nChamber %d has responded %d times on %d tracks", i+1, chamber[i], trackNb);
201 if (trackNb == 0)
202 cout<<"\nEfficiency = ? (IS UNKNOWN)\n";
203 else {
204 Double_t eff = chamber[i]*100./trackNb; cout<<"\nEfficiency = "<<eff<<" %\n";
205 }
206 }
207 cout<<"\n\n\n";
208 MUONLoader->UnloadTracks();
209}
210
211
212
213
214
215/*****************************************************************************************************************/
216/*****************************************************************************************************************/
217 /*EFFICIENCY IN FUNCTION OF THETA AND PHI*/
218
219void efficiencyThetaPhi( Int_t event2Check=0, char * filename="galice.root" )
220{
221 cout<<"\nChamber(s) chosen;\nFirst chamber:";
222 cin>>firstChamber;
223 cout<<"Last chamber:";
224 cin>>lastChamber;
225 cout<<"\n\n";
226
227 Int_t eff [54] = {0};
228 Int_t trackN [54] = {0};
229 Int_t chamber;
230 Int_t oldChamber;
231
232//Creating Run Loader and openning file containing Hits
233//--------------------------------------------------------------------------------------------------------------
234 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
235 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
236
237
238//Getting MUONLoader
239//--------------------------------------------------------------------------------------------------------------
240 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
241 MUONLoader -> LoadTracks("READ");
242 MUONLoader -> LoadRecPoints("READ");
243
244
245//--------------------------------------------------------------------------------------------------------------
246//Set mag field; waiting for mag field in CDB
f7a1cc68 247 if (!TGeoGlobalMagField::Instance()->GetField()) {
248 printf("Loading field map...\n");
249 AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
250 TGeoGlobalMagField::Instance()->SetField(field);
251 }
bd022fc7 252
253//--------------------------------------------------------------------------------------------------------------
254//Set Field Map for track extrapolation
f7a1cc68 255 AliMUONTrackExtrap::SetField()
bd022fc7 256
257
258//Creating a MUON data container
259//--------------------------------------------------------------------------------------------------------------
260 AliMUONData muondata(MUONLoader,"MUON","MUON");
261
262
263 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
264
265 //Loop on events
266 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
267 {
268 if ( event2Check!=0 )
269 iEvent=event2Check;
270 printf("\r>>> Event %d ",iEvent);
271 RunLoader->GetEvent(iEvent);
272
273 //Addressing
274 muondata.SetTreeAddress("RT");
275 muondata.GetRecTracks();
276 tracks = muondata.RecTracks();
277
278 //Loop on track
279 nTracks = (Int_t) tracks -> GetEntriesFast();
280 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
281 {
282 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 283 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 284 hits = track ->GetHitForRecAtHit() ;
285 chamber = firstChamber-1;
286 oldChamber = -1;
287 Int_t k = 1;
288
289 //Loop on hits
290 nHits = (Int_t) hits->GetEntriesFast();
291 for ( iHit = 0; iHit<nHits; ++iHit )
292 {
293 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
294 hit = (AliMUONHitForRec* ) hits -> At(iHit);
295 chamberNbr = hit -> GetChamberNumber();
296
297 if ( chamberNbr == oldChamber )
298 continue;
299 else {
300 oldChamber = chamberNbr;
301 if ( chamberNbr > chamber - k ) {
302 if ( chamber < lastChamber ) {
303 if ( chamberNbr == chamber ) {
304 //Positions
305 Double_t traX, traY, traZ;
306 Double_t hitX, hitY, hitZ;
307 Double_t aveX, aveY ;
308
309 //Angle (Phi)
310 Double_t phi = 0.;
311 Double_t theta = 0.;
312 Int_t iPhi = 0 ;
313 Int_t iTheta = 0 ;
314
315 traX = trackParam -> GetNonBendingCoor();
316 traY = trackParam -> GetBendingCoor() ;
317 traZ = trackParam -> GetZ() ;
318
319 hitX = hit -> GetNonBendingCoor();
320 hitY = hit -> GetBendingCoor() ;
321 hitZ = hit -> GetZ() ;
322
323 aveX = 1./2.*(traX + hitX);
324 aveY = 1./2.*(traY + hitY);
325
326 //The calculation of phi:
327 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
328
329 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
330 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
331
8ca62886 332 if ( phi < 0.) phi = 360 - fabs(phi);
bd022fc7 333 iPhi = int( phi/72. );
334 iTheta = int( theta );
335 if( theta > 10 ) iTheta = 9;
336 if( theta < 1 ) iTheta = 1;
337
338 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 1;
339 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
340 chamber = chamber + 1;
341 }
342
343 else {
344 //Positions
345 Double_t chamberZpos;
346 Double_t exXpos, exYpos;
347
348 //Angles
349 Double_t phi = 0.;
350 Double_t theta = 0.;
351 Int_t iPhi = 0 ;
352 Int_t iTheta = 0 ;
353
354 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
355 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
356 exXpos = (Double_t) trackParam->GetNonBendingCoor();
357 exYpos = (Double_t) trackParam->GetBendingCoor();
358
359 //The calculation of phi:
360 phi = 180. * kInvPi * (TMath::ATan2( exYpos, exXpos ));
361
362 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
363 theta = 180. * kInvPi * (TMath::ATan2( sqrt(exXpos*exXpos+exYpos*exYpos), - chamberZpos ));
364
8ca62886 365 if ( phi < 0.) phi = 360. - fabs(phi);
bd022fc7 366 iPhi = int( phi/72. );
367 iTheta = int( theta );
368 if( theta > 10 ) iTheta = 9;
369 if( theta < 1 ) iTheta = 1;
370
371 eff [9*iPhi+iTheta-1] = eff [9*iPhi+iTheta-1] + 0;
372 trackN [9*iPhi+iTheta-1] = trackN [9*iPhi+iTheta-1] + 1;
373 chamber = chamber + 1;
374 iHit = iHit - 1;
375 }
376 }
377 }
378 }
379
380 if ( iHit == nHits-1 ) {
381 if ( chamber == lastChamber )
382 continue;
383 else {
384 oldChamber = -1;
385 k = 5;
386 iHit = iHit-1;
387 }
388 }
389 }
390 //End Loop on hits
391
392 }
393 //End Loop on tracks
394
395 muondata.ResetRecTracks();
396 if ( event2Check!=0 )
397 iEvent=nEvents;
398 }
399 //End Loop on events
400
401
402 TCanvas * c1 = new TCanvas();
403 TGraph2D* effPhiTheta = new TGraph2D();
404 Double_t efficiency = 0;
405
406 if ( firstChamber == lastChamber )
407 {
408 for ( Int_t ph = 0; ph < 5; ++ph )
409 {
410 for ( Int_t th = 1; th < 10; ++th )
411 {
412 Int_t i = 9*ph+th-1;
413 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
414 cout<<"\nThe chamber "<<firstChamber<<" has responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
415
416 Double_t e = eff [i] ;
417 Double_t n = trackN[i] ;
418 Double_t p = ph*72.+36.;
419 Double_t t = th*1. +0.5;
420
421 if ( trackN[i] == 0 ) {
422 efficiency = 0.;
423 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
424 }
425 else {
426 efficiency = 100.*e/n;
427 cout<<"Efficiency = "<<efficiency<<" %\n";
428 }
429
430 effPhiTheta -> SetPoint( i, p, t, efficiency);
431 }
432 }
433 }
434
435 else{
436 for ( Int_t ph = 0; ph < 5; ++ph )
437 {
438 for ( Int_t th = 1; th < 10; ++th )
439 {
440 Int_t i = 9*ph+th-1;
441 cout<<"\nFor Phi = ["<<ph*72<<","<<ph*72+72<<"] and Theta = ["<<179-th<<","<<180-th<<"]:";
442 cout<<"\nThe chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff[i]<<" times on "<<trackN[i]<<" tracks\n";
443
444 Double_t e = eff [i] ;
445 Double_t n = trackN[i] ;
446 Double_t p = ph*72.+36.;
447 Double_t t = th*1. +0.5;
448
449 if ( trackN[i] == 0 ) {
450 efficiency = 0.;
451 cout<<"Efficiency = ? % ( IS UNKNOWN )\n";
452 }
453 else {
454 efficiency = 100.*e/n;
455 cout<<"Efficiency = "<<efficiency<<" %\n";
456 }
457
458 effPhiTheta -> SetPoint( i, p, t, efficiency);
459 }
460 }
461 }
462
463 gStyle->SetPalette(1);
464 effPhiTheta -> Draw("surf1");
465
466 cout<<"\n\n\n";
467 MUONLoader->UnloadTracks();
468 c1->Update();
469
470}
471
472
473
474
475
476
477
478/*****************************************************************************************************************/
479/*****************************************************************************************************************/
480 /*EFFICIENCY IN FUNCTION OF THETA OF INCIDENCE*/
481
482
483void efficiencyThetaI( Int_t event2Check=0, char * filename="galice.root" )
484{
485 cout<<"\nChamber(s) chosen;\nFirst chamber:";
486 cin>>firstChamber;
487 cout<<"Last chamber:";
488 cin>>lastChamber;
489 cout<<"\n\n";
490
491 Int_t eff [12] = {0};
492 Int_t trackN [12] = {0};
493 Int_t chamber;
494 Int_t oldChamber;
495
496
497//Creating Run Loader and openning file containing Hits
498//--------------------------------------------------------------------------------------------------------------
499 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
500 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
501
502
503//Getting MUONLoader
504//--------------------------------------------------------------------------------------------------------------
505 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
506 MUONLoader -> LoadTracks("READ");
507 MUONLoader -> LoadRecPoints("READ");
508
509
510//--------------------------------------------------------------------------------------------------------------
511//Set mag field; waiting for mag field in CDB
f7a1cc68 512 if (!TGeoGlobalMagField::Instance()->GetField()) {
513 printf("Loading field map...\n");
514 AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
515 TGeoGlobalMagField::Instance()->SetField(field);
516 } printf("Loading field map...\n");
517
bd022fc7 518
519//--------------------------------------------------------------------------------------------------------------
520//Set Field Map for track extrapolation
f7a1cc68 521 AliMUONTrackExtrap::SetField();
bd022fc7 522
523
524//Creating a MUON data container
525//--------------------------------------------------------------------------------------------------------------
526 AliMUONData muondata(MUONLoader,"MUON","MUON");
527
528
529 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
530
531 //Loop on events
532 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
533 {
534 if ( event2Check!=0 )
535 iEvent=event2Check;
536 printf("\r>>> Event %d ",iEvent);
537 RunLoader->GetEvent(iEvent);
538
539 //Addressing
540 muondata.SetTreeAddress("RT");
541 muondata.GetRecTracks();
542 tracks = muondata.RecTracks();
543
544 //Loop on track
545 nTracks = (Int_t) tracks -> GetEntriesFast();
546 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
547 {
548 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 549 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 550 hits = track ->GetHitForRecAtHit() ;
551 chamber = firstChamber - 1;
552 oldChamber = -1;
553 Int_t k = 1;
554
555 //Loop on hits
556 nHits = (Int_t) hits -> GetEntriesFast();
557 for ( iHit = 0; iHit < nHits; ++iHit )
558 {
559 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
560 hit = (AliMUONHitForRec*) hits -> At(iHit);
561 chamberNbr = hit -> GetChamberNumber();
562
563 if ( chamberNbr == oldChamber )
564 continue;
565 else {
566 oldChamber = chamberNbr;
567 if ( chamberNbr > chamber - k ) {
568 if ( chamber < lastChamber ) {
569 if ( chamberNbr == chamber ) {
570 //Momentum
571 Double_t Px, Py, Pz, Pr;
572
573 //Angle
574 Double_t theta;
575 Int_t iTheta;
576
577 Px = trackParam -> Px() ;
578 Py = trackParam -> Py() ;
579 Pz = trackParam -> Pz() ;
580 Pr = sqrt( Px*Px + Py*Py );
581
582 //The calculation of theta, the angle of incidence:
583 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
584
585 if ( theta < 79 ) iTheta = 11;
586 else if ( theta < 90 ) iTheta = int( theta - 79.);
587 else iTheta = 11;
588
589 eff [iTheta] = eff [iTheta] + 1;
590 trackN [iTheta] = trackN [iTheta] + 1;
591 chamber = chamber + 1;
592 }
593
594 else {
595 //Positions
596 Double_t chamberZpos;
597
598 //Momentum
599 Double_t Px, Py, Pz, Pr;
600
601 //Angles
602 Double_t theta = 0.;
603 Int_t iTheta = 0 ;
604
605 chamberZpos = AliMUONConstants::DefaultChamberZ(chamber);
606 AliMUONTrackExtrap::ExtrapToZ(trackParam,chamberZpos);
607
608 Px = trackParam -> Px() ;
609 Py = trackParam -> Py() ;
610 Pz = trackParam -> Pz() ;
611 Pr = sqrt( Px*Px + Py*Py );
612
613 //The calculation of thetaI, the angle of incidence:
614 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr));
615
616 if ( theta < 79 ) iTheta = 11;
617 else if ( theta < 90 ) iTheta = int( theta - 79.);
618 else iTheta = 11;
619
620 eff [iTheta] = eff [iTheta] + 0;
621 trackN [iTheta] = trackN [iTheta] + 1;
622 chamber = chamber + 1;
623 iHit = iHit - 1;
624 }
625 }
626 }
627 }
628
629 if ( iHit == nHits-1 ) {
630 if ( chamber == lastChamber )
631 continue;
632 else {
633 oldChamber = -1;
634 k = 5;
635 iHit = iHit-1;
636 }
637 }
638 }
639 //End loop on hits
640
641 }
642 //End Loop on tracks
643
644 muondata.ResetRecTracks();
645 if ( event2Check!=0 )
646 iEvent=nEvents;
647 }
648 //End Loop on events
649
650 Double_t t [11];
651 Double_t efficiency [11];
652 Int_t i = 11;
653
654 Int_t th;
655 TGraph * effTheta = new TGraph ();
656
657 if ( firstChamber == lastChamber ) {
658 for ( th = 0; th < 11; ++th )
659 {
660 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
661 cout<<"The chamber "<<firstChamber<<" has responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
662
663 t [th] = th + 79.5 ;
664 Double_t e = eff [th];
665 Double_t n = trackN [th];
666
667 if ( n == 0. ) {
668 efficiency [th] = 0.;
669 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
670 }
671 else {
672 efficiency [th] = 100.*e/n;
673 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
674 }
675 }
676 }
677
678 else{
679 for ( th = 0; th < 11; ++th )
680 {
681 cout<<"\nFor Theta (angle of incidence) = ["<<th+79<<","<<th+80<<"]:\n";
682 cout<<"The chambers "<<firstChamber<<" to "<<lastChamber<<" have responded "<<eff [th]<<" times on "<<trackN [th]<<" tracks\n";
683
684 t [th] = th + 79.5 ;
685 Double_t e = eff [th];
686 Double_t n = trackN [th];
687
688 if ( n == 0. ) {
689 efficiency [th] = 0.;
690 cout<<"Efficiency = ? % (IS UNKNOWN) \n";
691 }
692 else {
693 efficiency [th] = 100.*e/n;
694 cout<<"Efficiency = "<<efficiency [th]<<" %\n";
695 }
696 }
697 }
698
699 TCanvas * c1 = new TCanvas ();
700 effTheta = new TGraph( i, t, efficiency );
701
702 effTheta -> Draw("ALP");
703
704 MUONLoader->UnloadTracks();
705 c1->Update();
706}
707
708
709
710
711
712
713/*****************************************************************************************************************/
714/*****************************************************************************************************************/
715 /*RESOLUTION*/
716
717void resolution( Int_t event2Check=0, char * filename="galice.root" )
718{
719 cout<<"\nChamber(s) chosen;\nFirst chamber:";
720 cin>>firstChamber;
721 cout<<"Last chamber:";
722 cin>>lastChamber;
723 cout<<"\n\n";
724
725 TH1F * hDeltax;
726 TH1F * hDeltay;
727 TH2 * hDelta3D;
728
729 hDeltax = new TH1F ("hDeltax " , "No Interest", 100, -0.3, 0.3);
730 hDeltay = new TH1F ("hDeltay " , "No Interest", 500, -0.1, 0.1);
731 hDelta3D = new TH2F ("hDelta3D" , "No Interest", 100, -0.3, 0.3, 500, -0.1, 0.1);
732
733
734//Creating Run Loader and openning file containing Hits
735//--------------------------------------------------------------------------------------------------------------
736 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
737 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
738
739
740//Getting MUONLoader
741//--------------------------------------------------------------------------------------------------------------
742 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
743 MUONLoader -> LoadTracks("READ");
744 MUONLoader -> LoadRecPoints("READ");
745
746
747//Creating a MUON data container
748//--------------------------------------------------------------------------------------------------------------
749 AliMUONData muondata(MUONLoader,"MUON","MUON");
750
751
752 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
753
754 //Loop on events
755 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
756 {
757 if (event2Check!=0)
758 iEvent=event2Check;
759 printf("\r>>> Event %d ",iEvent);
760 RunLoader->GetEvent(iEvent);
761
762 //Addressing
763 muondata.SetTreeAddress("RT");
764 muondata.GetRecTracks();
765 tracks = muondata.RecTracks();
766
767 //Loop on track
768 nTracks = (Int_t) tracks -> GetEntriesFast();
769 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
770 {
771 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 772 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 773 hits = track ->GetHitForRecAtHit() ;
774
775 //Loop on hits
776 nHits = (Int_t) hits -> GetEntriesFast();
777 for ( iHit = 0; iHit < nHits; ++iHit )
778 {
779 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
780 hit = (AliMUONHitForRec*) hits -> At(iHit);
781 chamberNbr = hit -> GetChamberNumber();
782 if ( chamberNbr >= firstChamber-1 ) {
783 if ( chamberNbr < lastChamber ) {
784 //Positions
785 Double_t traX, traY;
786 Double_t hitX, hitY;
787
788 //Resolution
789 Double_t deltaX, deltaY;
790
791 traX = trackParam -> GetNonBendingCoor();
792 traY = trackParam -> GetBendingCoor();
793 hitX = hit -> GetNonBendingCoor();
794 hitY = hit -> GetBendingCoor();
795
796 deltaX = traX - hitX;
797 deltaY = traY - hitY;
798
799 hDeltax -> Fill (deltaX);
800 hDeltay -> Fill (deltaY);
801 hDelta3D -> Fill (deltaX, deltaY);
802 }
803 }
804 }
805 //End loop on hits
806 }
807 //End loop on tracks
808 muondata.ResetRecTracks();
809 if ( event2Check!=0 )
810 iEvent=nEvents;
811 }
812 //End loop on events
813
814 char hXtitle[80];
815 char hYtitle[80];
816 char h3title[80];
817
818 if ( firstChamber == lastChamber ) {
819 sprintf(hXtitle,"DeltaX Chamber %d",firstChamber);
820 sprintf(hYtitle,"DeltaY Chamber %d",lastChamber);
821 sprintf(h3title,"DeltaX and Delta Y Chamber %d",lastChamber);
822 }
823 else{
824 sprintf(hXtitle,"DeltaX Chamber %d to %d",firstChamber,lastChamber);
825 sprintf(hYtitle,"DeltaY Chamber %d to %d",firstChamber,lastChamber);
826 sprintf(h3title,"DeltaX and Delta Y Chamber %d to %d",firstChamber,lastChamber);
827 }
828
829
830 hDeltax -> SetTitle (hXtitle);
831 hDeltay -> SetTitle (hYtitle);
832 hDelta3D -> SetTitle (h3title);
833
834 Double_t rmsX = hDeltax -> GetRMS ();
835 Double_t errX = hDeltax -> GetRMSError();
836
837 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
838 hDeltay -> Fit("fY","R","E");
839 Double_t sigY = fY -> GetParameter(2);
840 Double_t errY = fY -> GetParError (2);
841
842 if ( firstChamber == lastChamber ) {
843 cout<<"\nThe resolution (Non Bending plane) of the chamber "<<firstChamber<<" = "<<rmsX<<" cm +-"<<errX;
844 cout<<"\nThe resolution (Bending plane) of the chamber "<<firstChamber<<" = "<<sigY<<" cm +- "<<errY;
845 }
846 else {
847 cout<<"\nThe resolution (Non Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<rmsX<<" cm +- "<<errX;
848 cout<<"\nThe resolution (Bending plane) of the chambers "<<firstChamber<<" to "<<lastChamber<<" = "<<sigY<<" cm +- "<<errY;
849 }
850
851 cout<<"\n\n\n";
852 MUONLoader->UnloadTracks();
853
854}
855
856
857
858
859
860
861
862/*****************************************************************************************************************/
863/*****************************************************************************************************************/
864 /*RESOLUTION IN FUNCTION OF PHI*/
865
866void resolutionPhi( Int_t event2Check=0, char * filename="galice.root" )
867{
868 cout<<"\nChamber(s) chosen;\nFirst chamber:";
869 cin>>firstChamber;
870 cout<<"Last chamber:";
871 cin>>lastChamber;
872 cout<<"\n\n";
873
874 TH1F * hDeltaX[5];
875 TH1F * hDeltaY[5];
876
877 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
878 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
879 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
880 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
881 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
882
883 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
884 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
885 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
886 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
887 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
888
889
890//Creating Run Loader and openning file containing Hits
891//--------------------------------------------------------------------------------------------------------------
892 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
893 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
894
895
896//Getting MUONLoader
897//--------------------------------------------------------------------------------------------------------------
898 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
899 MUONLoader -> LoadTracks("READ");
900 MUONLoader -> LoadRecPoints("READ");
901
902
903//Creating a MUON data container
904//--------------------------------------------------------------------------------------------------------------
905 AliMUONData muondata(MUONLoader,"MUON","MUON");
906
907
908 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
909
910 //Loop on events
911 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
912 {
913 if ( event2Check!=0 )
914 iEvent=event2Check;
915 printf("\r>>> Event %d ",iEvent);
916 RunLoader->GetEvent(iEvent);
917
918 //Addressing
919 muondata.SetTreeAddress("RT");
920 muondata.GetRecTracks();
921 tracks = muondata.RecTracks();
922
923 //Loop on track
924 nTracks = (Int_t) tracks -> GetEntriesFast();
925 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
926 {
927 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 928 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 929 hits = track ->GetHitForRecAtHit() ;
930
931 //Loop on hits
932 nHits = (Int_t) hits -> GetEntriesFast();
933 for ( iHit = 0; iHit < nHits; ++iHit )
934 {
935 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
936 hit = (AliMUONHitForRec* ) hits -> At(iHit);
937 chamberNbr = hit -> GetChamberNumber();
938 if ( chamberNbr >= firstChamber -1 ) {
939 if ( chamberNbr < lastChamber ) {
940 //Positions
941 Double_t traX, traY;
942 Double_t hitX, hitY;
943 Double_t aveY, aveX;
944
945 //Angles
946 Double_t phi;
947 Int_t iPhi;
948
949 //Resolution
950 Double_t deltaX;
951 Double_t deltaY;
952
953 traX = trackParam -> GetNonBendingCoor();
954 traY = trackParam -> GetBendingCoor() ;
955
956 hitX = hit -> GetNonBendingCoor();
957 hitY = hit -> GetBendingCoor() ;
958
959 aveX = 1./2.*(traX + hitX);
960 aveY = 1./2.*(traY + hitY);
961
962 //The calculation of phi:
963 phi = 180. * kInvPi * (TMath::ATan2( aveY, aveX ));
964
8ca62886 965 if ( phi < 0.) phi = 360. - fabs(phi);
bd022fc7 966 iPhi = int( phi/72. );
967
968 deltaX = traX - hitX;
969 deltaY = traY - hitY;
970
971 hDeltaX [iPhi] -> Fill( deltaX );
972 hDeltaY [iPhi] -> Fill( deltaY );
973 }
974 }
975
976 }
977 //End loop on hits
978
979 }
980 //End loop on tracks
981 muondata.ResetRecTracks();
982 if ( event2Check!=0 )
983 iEvent=nEvents;
984
985 }
986 //End loop on events
987
988
989 Int_t iPhi;
990 Int_t iPhiMax = 5;
991 Int_t phiMin, phiMax;
992
993 Float_t phi[5];
994 Float_t sigmaY[5];
995 Float_t errSY [5];
996 Float_t rmsX [5];
997 Float_t errSX [5];
998 Float_t errPh [5];
999
1000 for ( iPhi = 0; iPhi < iPhiMax; ++iPhi )
1001 {
1002 char hXtitle[80];
1003 char hYtitle[80];
1004
1005 phiMin = 72*iPhi ;
1006 phiMax = 72*iPhi + 72;
1007
1008 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1009
1010 if ( firstChamber == lastChamber ) {
1011 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d",phiMin,phiMax,firstChamber);
1012 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d",phiMin,phiMax,lastChamber);
1013 }
1014 else{
1015 sprintf(hXtitle,"DeltaX (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1016 sprintf(hYtitle,"DeltaY (phi=[%d, %d]) Chamber %d to %d",phiMin,phiMax,firstChamber,lastChamber);
1017 }
1018
1019 hDeltaY [iPhi] -> SetTitle(hYtitle);
1020 hDeltaX [iPhi] -> SetTitle(hXtitle);
1021
1022 hDeltaY [iPhi] -> Fit("fY","R","E");
1023 sigmaY [iPhi] = fY -> GetParameter(2) ;
1024 errSY [iPhi] = fY -> GetParError(2) ;
1025
1026 rmsX [iPhi] = hDeltaX [iPhi] -> GetRMS();
1027 errSX [iPhi] = hDeltaX [iPhi] -> GetRMSError();
1028
1029 phi [iPhi] = 72*iPhi + 36 ;
1030 errPh [iPhi] = 36;
1031 }
1032
1033 //--------------------------------------------------------------------------------------------------------------
1034 //For plotting resolution in function of the angle of incidence
1035
1036 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1037 c1-> Divide(1,2);
1038 c1->cd(1);
1039
1040 TGraphErrors * GraphX = new TGraphErrors( iPhiMax, phi, rmsX, errPh, errSX);
1041 GraphX->SetTitle("Resolution en X (Phi)");
1042 GraphX->Draw("ALP");
1043
1044 c1->cd(2);
1045 TGraphErrors * GraphY = new TGraphErrors( iPhiMax, phi, sigmaY, errPh, errSY);
1046 GraphY->SetTitle("Resolution en Y (Phi)");
1047 GraphY->Draw("ALP");
1048
1049 cout<<"\n\n\n";
1050
1051 MUONLoader->UnloadTracks();
1052
1053}
1054
1055
1056
1057
1058
1059
1060
1061/*****************************************************************************************************************/
1062/*****************************************************************************************************************/
1063 /*RESOLUTION IN FUNCTION OF THETA*/
1064
1065void resolutionTheta( Int_t event2Check=0, char * filename="galice.root" )
1066{
1067 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1068 cin>>firstChamber;
1069 cout<<"Last chamber:";
1070 cin>>lastChamber;
1071 cout<<"\n\n";
1072
1073 TH1F * hDeltaX[9];
1074 TH1F * hDeltaY[9];
1075
1076 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1077 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1078 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1079 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1080 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1081 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1082 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1083 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1084 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1085
1086 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1087 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1088 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1089 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1090 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1091 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1092 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1093 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1094 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1095
1096
1097//Creating Run Loader and openning file containing Hits
1098//--------------------------------------------------------------------------------------------------------------
1099 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1100 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1101
1102
1103//Getting MUONLoader
1104//--------------------------------------------------------------------------------------------------------------
1105 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1106 MUONLoader -> LoadTracks("READ");
1107 MUONLoader -> LoadRecPoints("READ");
1108
1109
1110//Creating a MUON data container
1111//--------------------------------------------------------------------------------------------------------------
1112 AliMUONData muondata(MUONLoader,"MUON","MUON");
1113
1114
1115 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1116
1117 //Loop on events
1118 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1119 {
1120 if ( event2Check!=0 )
1121 iEvent=event2Check;
1122 printf("\r>>> Event %d ",iEvent);
1123 RunLoader->GetEvent(iEvent);
1124
1125 //Addressing
1126 muondata.SetTreeAddress("RT");
1127 muondata.GetRecTracks();
1128 tracks = muondata.RecTracks();
1129
1130 //Loop on track
1131 nTracks = (Int_t) tracks -> GetEntriesFast();
1132 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1133 {
1134 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 1135 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 1136 hits = track ->GetHitForRecAtHit() ;
1137
1138 //Loop on hits
1139 nHits = (Int_t) hits -> GetEntriesFast();
1140 for ( iHit = 0; iHit < nHits; ++iHit )
1141 {
1142 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1143 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1144 chamberNbr = hit -> GetChamberNumber();
1145 if ( chamberNbr >= firstChamber -1 ) {
1146 if ( chamberNbr < lastChamber ) {
1147 //Positions
1148 Double_t traX, traY;
1149 Double_t hitX, hitY, hitZ;
1150 Double_t aveY, aveX;
1151
1152 //Angles
1153 Double_t theta;
1154 Int_t iTheta;
1155
1156 //Resolution
1157 Double_t deltaX;
1158 Double_t deltaY;
1159
1160 traX = trackParam -> GetNonBendingCoor();
1161 traY = trackParam -> GetBendingCoor() ;
1162
1163 hitX = hit -> GetNonBendingCoor();
1164 hitY = hit -> GetBendingCoor() ;
1165 hitZ = hit -> GetZ();
1166
1167 aveX = 1./2.*(traX + hitX);
1168 aveY = 1./2.*(traY + hitY);
1169
1170 //The calculation of theta, theta is in fact 180 - Theta(The polar angle):
1171 theta = 180. * kInvPi * (TMath::ATan2( sqrt(aveX*aveX+aveY*aveY), -hitZ ));
1172
1173 iTheta = int( theta );
1174 if( theta > 10 ) iTheta = 9;
1175 if( theta < 1 ) iTheta = 1;
1176
1177 deltaX = traX - hitX;
1178 deltaY = traY - hitY;
1179
1180 hDeltaX [iTheta-1] -> Fill( deltaX );
1181 hDeltaY [iTheta-1] -> Fill( deltaY );
1182 }
1183 }
1184
1185 }
1186 //End loop on hits
1187
1188 }
1189 //End loop on tracks
1190 muondata.ResetRecTracks();
1191 if ( event2Check!=0 )
1192 iEvent=nEvents;
1193
1194 }
1195 //End loop on events
1196
1197
1198 Int_t iTheta;
1199 Int_t iThetaMax = 9;
1200 Int_t thetaMin, thetaMax;
1201
1202 Float_t theta [9];
1203 Float_t sigmaY[9];
1204 Float_t errSY [9];
1205 Float_t rmsX [9];
1206 Float_t errSX [9];
1207 Float_t ErrTh [9];
1208
1209 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1210 {
1211 char hXtitle[80];
1212 char hYtitle[80];
1213
1214 //To find the polar angle
1215 thetaMin = 178 - iTheta;
1216 thetaMax = 179 - iTheta;
1217
1218 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1219
1220 if ( firstChamber == lastChamber ) {
1221 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1222 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1223 }
1224 else{
1225 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1226 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1227 }
1228
1229 hDeltaY [iTheta] -> SetTitle(hYtitle);
1230 hDeltaX [iTheta] -> SetTitle(hXtitle);
1231
1232 hDeltaY [iTheta] -> Fit("fY","R","E");
1233 sigmaY [iTheta] = fY -> GetParameter(2) ;
1234 errSY [iTheta] = fY -> GetParError(2) ;
1235
1236 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1237 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1238
1239 theta [iTheta] = 178.5 - iTheta;
1240 ErrTh [iTheta] = 0.5;
1241 }
1242
1243 //--------------------------------------------------------------------------------------------------------------
1244 //For plotting resolution in function of the angle of incidence
1245
1246 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1247 c1-> Divide(1,2);
1248 c1->cd(1);
1249
1250 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, ErrTh, errSX);
1251 GraphX->SetTitle("Resolution en X (Theta)");
1252 GraphX->Draw("ALP");
1253
1254 c1->cd(2);
1255 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, ErrTh, errSY);
1256 GraphY->SetTitle("Resolution en Y (Theta)");
1257 GraphY->Draw("ALP");
1258
1259 cout<<"\n\n\n";
1260 MUONLoader->UnloadTracks();
1261
1262}
1263
1264
1265
1266
1267
1268/*****************************************************************************************************************/
1269/*****************************************************************************************************************/
1270 /*RESOLUTION IN FUNCTION OF THETA OF INCIDENCE*/
1271
1272void resolutionThetaI( Int_t event2Check=0, char * filename="galice.root" )
1273{
1274 cout<<"\nChamber(s) chosen;\nFirst chamber:";
1275 cin>>firstChamber;
1276 cout<<"Last chamber:";
1277 cin>>lastChamber;
1278 cout<<"\n\n";
1279
1280 TH1F * hDeltaX[11];
1281 TH1F * hDeltaY[11];
1282
1283 hDeltaX[0] = new TH1F ("hDeltaX0" , "No Interest", 100, -0.3, 0.3);
1284 hDeltaX[1] = new TH1F ("hDeltaX1" , "No Interest", 100, -0.3, 0.3);
1285 hDeltaX[2] = new TH1F ("hDeltaX2" , "No Interest", 100, -0.3, 0.3);
1286 hDeltaX[3] = new TH1F ("hDeltaX3" , "No Interest", 100, -0.3, 0.3);
1287 hDeltaX[4] = new TH1F ("hDeltaX4" , "No Interest", 100, -0.3, 0.3);
1288 hDeltaX[5] = new TH1F ("hDeltaX5" , "No Interest", 100, -0.3, 0.3);
1289 hDeltaX[6] = new TH1F ("hDeltaX6" , "No Interest", 100, -0.3, 0.3);
1290 hDeltaX[7] = new TH1F ("hDeltaX7" , "No Interest", 100, -0.3, 0.3);
1291 hDeltaX[8] = new TH1F ("hDeltaX8" , "No Interest", 100, -0.3, 0.3);
1292 hDeltaX[9] = new TH1F ("hDeltaX9" , "No Interest", 100, -0.3, 0.3);
1293 hDeltaX[10]= new TH1F ("hDeltaX10", "No Interest", 100, -0.3, 0.3);
1294
1295 hDeltaY[0] = new TH1F ("hDeltaY0" , "No Interest", 500, -0.1, 0.1);
1296 hDeltaY[1] = new TH1F ("hDeltaY1" , "No Interest", 500, -0.1, 0.1);
1297 hDeltaY[2] = new TH1F ("hDeltaY2" , "No Interest", 500, -0.1, 0.1);
1298 hDeltaY[3] = new TH1F ("hDeltaY3" , "No Interest", 500, -0.1, 0.1);
1299 hDeltaY[4] = new TH1F ("hDeltaY4" , "No Interest", 500, -0.1, 0.1);
1300 hDeltaY[5] = new TH1F ("hDeltaY5" , "No Interest", 500, -0.1, 0.1);
1301 hDeltaY[6] = new TH1F ("hDeltaY6" , "No Interest", 500, -0.1, 0.1);
1302 hDeltaY[7] = new TH1F ("hDeltaY7" , "No Interest", 500, -0.1, 0.1);
1303 hDeltaY[8] = new TH1F ("hDeltaY8" , "No Interest", 500, -0.1, 0.1);
1304 hDeltaY[9] = new TH1F ("hDeltaY9" , "No Interest", 500, -0.1, 0.1);
1305 hDeltaY[10]= new TH1F ("hDeltaY10", "No Interest", 500, -0.1, 0.1);
1306
1307
1308//Creating Run Loader and openning file containing Hits
1309//--------------------------------------------------------------------------------------------------------------
1310 AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
1311 if (RunLoader == 0x0) {printf(">>> Error : Error Opening %s file \n",filename);return;}
1312
1313
1314//Getting MUONLoader
1315//--------------------------------------------------------------------------------------------------------------
1316 AliLoader * MUONLoader = RunLoader -> GetLoader("MUONLoader");
1317 MUONLoader -> LoadTracks("READ");
1318 MUONLoader -> LoadRecPoints("READ");
1319
1320
1321//Creating a MUON data container
1322//--------------------------------------------------------------------------------------------------------------
1323 AliMUONData muondata(MUONLoader,"MUON","MUON");
1324
1325
1326 nEvents = (Int_t) RunLoader -> GetNumberOfEvents();
1327
1328 //Loop on events
1329 for ( iEvent = 0; iEvent < nEvents; ++iEvent )
1330 {
1331 if ( event2Check!=0 )
1332 iEvent=event2Check;
1333 printf("\r>>> Event %d ",iEvent);
1334 RunLoader->GetEvent(iEvent);
1335
1336 //Addressing
1337 muondata.SetTreeAddress("RT");
1338 muondata.GetRecTracks();
1339 tracks = muondata.RecTracks();
1340
1341 //Loop on track
1342 nTracks = (Int_t) tracks -> GetEntriesFast();
1343 for ( iTrack = 0; iTrack < nTracks; ++iTrack )
1344 {
1345 track = (AliMUONTrack*) tracks ->At(iTrack) ;
96ebe67e 1346 trackParams = track ->GetTrackParamAtCluster();
bd022fc7 1347 hits = track ->GetHitForRecAtHit() ;
1348
1349 //Loop on hits
1350 nHits = (Int_t) hits -> GetEntriesFast();
1351 for ( iHit = 0; iHit < nHits; ++iHit )
1352 {
1353 trackParam = (AliMUONTrackParam*) trackParams -> At(iHit);
1354 hit = (AliMUONHitForRec* ) hits -> At(iHit);
1355 chamberNbr = hit -> GetChamberNumber();
1356 if ( chamberNbr >= firstChamber - 1 ) {
1357 if ( chamberNbr < lastChamber ) {
1358 //Momentum
1359 Double_t Px, Py, Pz, Pr;
1360
1361 //Positions
1362 Double_t traX, traY;
1363 Double_t hitX, hitY;
1364
1365 //Resolution
1366 Double_t deltaX;
1367 Double_t deltaY;
1368
1369 //Angle
1370 Double_t theta;
1371 Int_t iTheta;
1372
1373 Px = trackParam -> Px() ;
1374 Py = trackParam -> Py() ;
1375 Pz = trackParam -> Pz() ;
1376 Pr = sqrt( Px*Px + Py*Py );
1377
1378 traX = trackParam -> GetNonBendingCoor();
1379 traY = trackParam -> GetBendingCoor();
1380
1381 hitX = hit -> GetNonBendingCoor();
1382 hitY = hit -> GetBendingCoor();
1383
1384 //The calculation of theta, the angle of incidence:
1385 theta = 180. * kInvPi * (TMath::ATan( - Pz / Pr))
1386;
1387 if ( theta < 79 ) iTheta = 0;
1388 else iTheta = int( theta - 79.);
1389
1390 deltaX = traX - hitX;
1391 deltaY = traY - hitY;
1392
1393 hDeltaX [iTheta] -> Fill (deltaX);
1394 hDeltaY [iTheta] -> Fill (deltaY);
1395 }
1396 }
1397
1398 }
1399 //End loop on hits
1400
1401 }
1402 //End loop on tracks
1403 muondata.ResetRecTracks();
1404 if ( event2Check!=0 )
1405 iEvent=nEvents;
1406
1407 }
1408 //End loop on events
1409
1410
1411 Int_t iTheta;
1412 Int_t iThetaMax = 11;
1413 Int_t thetaMin, thetaMax;
1414
1415 Float_t theta [11];
1416 Float_t sigmaY[11];
1417 Float_t errSY [11];
1418 Float_t rmsX [11];
1419 Float_t errSX [11];
1420 Float_t errTh [11];
1421
1422 for ( iTheta = 0; iTheta < iThetaMax; ++iTheta )
1423 {
1424 char hXtitle[80];
1425 char hYtitle[80];
1426
1427 thetaMin = iTheta + 79;
1428 thetaMax = iTheta + 80;
1429
1430 TF1 *fY = new TF1("fY","gaus",-0.3, 0.3);
1431
1432 if ( firstChamber == lastChamber ) {
1433 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,firstChamber);
1434 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d",thetaMin,thetaMax,lastChamber);
1435 }
1436 else{
1437 sprintf(hXtitle,"DeltaX (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1438 sprintf(hYtitle,"DeltaY (theta=[%d, %d]) Chamber %d to %d",thetaMin,thetaMax,firstChamber,lastChamber);
1439 }
1440
1441 hDeltaY [iTheta] -> SetTitle(hYtitle);
1442 hDeltaX [iTheta] -> SetTitle(hXtitle);
1443
1444 hDeltaY [iTheta] -> Fit("fY","R","E");
1445 sigmaY [iTheta] = fY -> GetParameter(2) ;
1446 errSY [iTheta] = fY -> GetParError(2) ;
1447
1448 rmsX [iTheta] = hDeltaX [iTheta] -> GetRMS();
1449 errSX [iTheta] = hDeltaX [iTheta] -> GetRMSError();
1450
1451 theta [iTheta] = iTheta + 79.5 ;
1452 errTh [iTheta] = 0.5;
1453 }
1454
1455 //--------------------------------------------------------------------------------------------------------------
1456 //For plotting resolution in function of the angle of incidence
1457
1458 TCanvas * c1 = new TCanvas("c1","",200,10,700,500);
1459 c1-> Divide(1,2);
1460 c1->cd(1);
1461
1462 TGraphErrors * GraphX = new TGraphErrors( iThetaMax, theta, rmsX, errTh, errSX);
1463 GraphX->SetTitle("Resolution en X (Theta)");
1464 GraphX->Draw("ALP");
1465
1466 c1->cd(2);
1467 TGraphErrors * GraphY = new TGraphErrors( iThetaMax, theta, sigmaY, errTh, errSY);
1468 GraphY->SetTitle("Resolution en Y (Theta)");
1469 GraphY->Draw("ALP");
1470
1471 cout<<"\n\n\n";
1472 MUONLoader->UnloadTracks();
1473
1474}
1475
1476
1477
1478