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