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