/[PAMELA software]/PamelaLevel2/doc/examples/example2.C
ViewVC logotype

Annotation of /PamelaLevel2/doc/examples/example2.C

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Fri Jun 16 16:43:55 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
Branch point for: PamelaLevel2
File MIME type: text/plain
Initial revision

1 pam-fi 1.1 //
2     // Example to get pamela tracks and integrate the trajectory in the apparatus
3     //
4     // The tracking method evaluates, besides the track intersection coordinates at given z-coordinates,
5     // the track length (total and between two given points along the trajectory) and the projected angles.
6     //
7     example2(TString file){
8    
9     gROOT->Reset();
10     //
11     // create some histograms
12     //
13     TH1F* s1x = new TH1F("s1x","Track x impact position on S1",100,-30.,30.);
14     TH1F* s1y = new TH1F("s1y","Track y impact position on S1",100,-30.,30.);
15     TH1F* s2x = new TH1F("s2x","Track x impact position on S2",100,-30.,30.);
16     TH1F* s2y = new TH1F("s2y","Track y impact position on S2",100,-30.,30.);
17     TH1F* s3x = new TH1F("s3x","Track x impact position on S3",100,-30.,30.);
18     TH1F* s3y = new TH1F("s3y","Track y impact position on S3",100,-30.,30.);
19    
20     TH1F* s1tx = new TH1F("s1tx","Track x projected angle on S1",50,-20.,20.);
21     TH1F* s1ty = new TH1F("s1ty","Track y projected angle on S1",50,-20.,20.);
22     TH1F* s2tx = new TH1F("s2tx","Track x projected angle on S2",50,-20.,20.);
23     TH1F* s2ty = new TH1F("s2ty","Track y projected angle on S2",50,-20.,20.);
24     TH1F* s3tx = new TH1F("s3tx","Track x projected angle on S3",50,-20.,20.);
25     TH1F* s3ty = new TH1F("s3ty","Track y projected angle on S3",50,-20.,20.);
26    
27     TH1F* trl = new TH1F("trl","Track length",100,0.,100.);
28     TH1F* trl12 = new TH1F("trl12","Track length S1-S2",100,0.,100.);
29     TH1F* trl23 = new TH1F("trl23","Track length S2-S3",100,0.,100.);
30    
31    
32     //
33     PamLevel2* pam_event = new PamLevel2(); // << create pamela event
34     //
35     TFile f(file);
36     //
37     TTree *T = pam_event->LoadPamTrees(&f); // << load Pamela trees from file f
38     Int_t nevent = T->GetEntries();
39     // ********************
40     // load magnetic field
41     // ********************
42     pam_event->LoadField("./magnetic-field/"); //
43     // ********************
44     // initialize some trajectories
45     // ********************
46     // << default trajectory is created with 10 points between the upper and lower tracker planes
47     Trajectory *tr1 = new Trajectory() ; // << create a default trajectory
48     // << the number of points can be set by the user
49     Trajectory *tr2 = new Trajectory(100) ; // << create a trajectory with 100 points
50     // << and also the z-coordinates can be set by the user.
51     // << for example if we want to evaluate the track intersection points in the TOF planes
52     // << we can define the following trajectory
53     Int_t nz = 6; Float_t zin[6]; // << define TOF z-coordinates
54     for(Int_t ip=0; ip<nz; ip++)
55     zin[ip] = pam_event->GetZTOF(pam_event->GetToFPlaneID(ip)); // << read ToF plane z-coordinates
56     Trajectory *tr = new Trajectory(nz,zin); // << create a trajectory in the apparatus
57     //
58     cout << endl<< " Start loop over events ";
59     for (Int_t i=0; i<nevent;i++){
60     //
61     T->GetEntry(i);
62     //
63     if(pam_event->GetNTracks()==1){ // << select events with only one track
64     //
65     TrkTrack *track = pam_event->GetTrack(0); // << retrieve the track
66     //
67     // << perform some track selection
68     //
69     if(
70     track->chi2 > 0 &&
71     track->chi2 < 100 &&
72     track->GetNX() >= 4 &&
73     track->GetNY() >= 3 &&
74     true
75     ){
76    
77     cout << endl<< "***** First trajectory"<< endl;
78     track->DoTrack2(tr1); // << calculate the first trajectory in magnetic field
79     // tr1->Dump(); // dump the trajectory
80     cout << "Length: "<< tr1->GetLength()<<endl; // << get the track length
81     //
82     cout << endl<< "***** Second trajectory"<< endl;
83     track->DoTrack2(tr2); // << calculate trajectory in magnetic field
84     cout << "Length: "<< tr2->GetLength()<<endl; // << get the track length
85     //
86     cout << endl<< "***** Third trajectory"<< endl;
87     track->DoTrack2(tr); // << calculate trajectory in magnetic field
88     // tr->Dump();
89     cout << "Length (S11-S32): "<< tr->GetLength()<<endl;
90     // << The length can be evaluated also bewteen two planes set by the user
91     cout << "Length between S11-S21: "<< tr->GetLength(0,2)<< endl;
92     cout << "Length between S21-S31: "<< tr->GetLength(2,4)<< endl;
93     cout << "Length between S11-S31: "<< tr->GetLength(0,4)<< endl;
94     // now fills some histos:
95     // store calculated coordinates
96     s1x->Fill( tr->x[0] );
97     s1y->Fill( tr->y[0] );
98     s2x->Fill( tr->x[1] );
99     s2y->Fill( tr->y[1] );
100     s3x->Fill( tr->x[2] );
101     s3y->Fill( tr->y[2] );
102     // store calculated projected angles
103     s1tx->Fill( tr->thx[0] );
104     s1ty->Fill( tr->thy[0] );
105     s2tx->Fill( tr->thx[1] );
106     s2ty->Fill( tr->thy[1] );
107     s3tx->Fill( tr->thx[2] );
108     s3ty->Fill( tr->thy[2] );
109    
110     // store calculated track lengths
111     trl->Fill( tr->GetLength() );
112     trl12->Fill( tr->GetLength(0,1) );
113     trl23->Fill( tr->GetLength(1,2) );
114     };
115     };
116     };
117     cout << endl << " Done "<< endl<<endl;
118     //
119     // close file
120     //
121     f.Close();
122     //
123     // plot histos
124     //
125     }

  ViewVC Help
Powered by ViewVC 1.1.23