doc/examples/example2.C

//      
//      Example to get pamela tracks and integrate the trajectory in the apparatus
//
//      The tracking method evaluates, besides the track intersection coordinates at given z-coordinates, 
//      the track length (total and between two given points along the trajectory) and the projected angles.
//
example2(TString file){
        
        gROOT->Reset();
        //
        // create some histograms
        //
        TH1F* s1x = new TH1F("s1x","Track x impact position on S1",100,-30.,30.);
        TH1F* s1y = new TH1F("s1y","Track y impact position on S1",100,-30.,30.);
        TH1F* s2x = new TH1F("s2x","Track x impact position on S2",100,-30.,30.);
        TH1F* s2y = new TH1F("s2y","Track y impact position on S2",100,-30.,30.);
        TH1F* s3x = new TH1F("s3x","Track x impact position on S3",100,-30.,30.);
        TH1F* s3y = new TH1F("s3y","Track y impact position on S3",100,-30.,30.);
        
        TH1F* s1tx = new TH1F("s1tx","Track x projected angle on S1",50,-20.,20.);
        TH1F* s1ty = new TH1F("s1ty","Track y projected angle on S1",50,-20.,20.);
        TH1F* s2tx = new TH1F("s2tx","Track x projected angle on S2",50,-20.,20.);
        TH1F* s2ty = new TH1F("s2ty","Track y projected angle on S2",50,-20.,20.);
        TH1F* s3tx = new TH1F("s3tx","Track x projected angle on S3",50,-20.,20.);
        TH1F* s3ty = new TH1F("s3ty","Track y projected angle on S3",50,-20.,20.);      
        
        TH1F* trl   = new TH1F("trl","Track length",100,0.,100.);
        TH1F* trl12 = new TH1F("trl12","Track length S1-S2",100,0.,100.);
        TH1F* trl23 = new TH1F("trl23","Track length S2-S3",100,0.,100.);
        
        
        //
        PamLevel2*  pam_event  = new PamLevel2();       // << create pamela event
        //
        TFile f(file);
        //
        TTree *T = pam_event->LoadPamTrees(&f);         // << load Pamela trees from file f
        Int_t nevent = T->GetEntries(); 
        // ********************
        // load magnetic field
        // ********************
        pam_event->LoadField("./magnetic-field/");      // 
        // ********************
        // initialize some trajectories
        // ********************
        // << default trajectory is created with 10 points between the upper and lower tracker planes
        Trajectory *tr1 = new Trajectory() ;                                   // << create a default trajectory 
        // << the number of points can be set by the user
        Trajectory *tr2 = new Trajectory(100) ;                                // << create a trajectory with 100 points
        // << and also the z-coordinates can be set by the user.
        // << for example if we want to evaluate the track intersection points in the TOF planes
        // << we can define the following trajectory
        Int_t nz = 6; Float_t zin[6];                                          // << define TOF z-coordinates
        for(Int_t ip=0; ip<nz; ip++) 
               zin[ip] = pam_event->GetZTOF(pam_event->GetToFPlaneID(ip));     // << read ToF plane z-coordinates
        Trajectory *tr = new Trajectory(nz,zin);                               // << create a trajectory in the apparatus
        //
        cout << endl<< " Start loop over events   ";
        for (Int_t i=0; i<nevent;i++){
                //
                T->GetEntry(i);
                //
                if(pam_event->GetNTracks()==1){                 // << select events with only one track
                    //
                    TrkTrack *track = pam_event->GetTrack(0);   // << retrieve the track
                        //
                        // << perform some track selection
                        //
                        if(                                     
                                track->chi2 > 0    && 
                                track->chi2 < 100  &&
                                track->GetNX() >= 4  &&
                                track->GetNY() >= 3  &&
                                true
                                ){      
                                
                            cout << endl<< "***** First trajectory"<< endl;
                            track->DoTrack2(tr1);       // << calculate the first trajectory in magnetic field
//                          tr1->Dump();                // dump the trajectory
                            cout << "Length: "<< tr1->GetLength()<<endl;  // << get the track length 
                            //
                            cout << endl<< "***** Second trajectory"<< endl;
                            track->DoTrack2(tr2);       // << calculate trajectory in magnetic field
                            cout << "Length: "<< tr2->GetLength()<<endl; // << get the track length 
                            //
                            cout << endl<< "***** Third trajectory"<< endl;
                            track->DoTrack2(tr);        // << calculate trajectory in magnetic field                            
//                            tr->Dump();
                            cout << "Length (S11-S32): "<< tr->GetLength()<<endl;
                            // << The length can be evaluated also bewteen two planes set by the user
                            cout << "Length between S11-S21: "<< tr->GetLength(0,2)<< endl;
                            cout << "Length between S21-S31: "<< tr->GetLength(2,4)<< endl;
                            cout << "Length between S11-S31: "<< tr->GetLength(0,4)<< endl;
                            // now fills some histos:
                            // store calculated coordinates
                            s1x->Fill( tr->x[0] );      
                            s1y->Fill( tr->y[0] );
                            s2x->Fill( tr->x[1] );
                            s2y->Fill( tr->y[1] );
                            s3x->Fill( tr->x[2] );
                            s3y->Fill( tr->y[2] );
                            // store calculated projected angles
                            s1tx->Fill( tr->thx[0] );   
                            s1ty->Fill( tr->thy[0] );
                            s2tx->Fill( tr->thx[1] );
                            s2ty->Fill( tr->thy[1] );
                            s3tx->Fill( tr->thx[2] );
                            s3ty->Fill( tr->thy[2] );
                            
                            // store calculated track lengths
                            trl->Fill( tr->GetLength() );
                            trl12->Fill( tr->GetLength(0,1) );
                            trl23->Fill( tr->GetLength(1,2) );
                            };
                };
        };
        cout << endl << " Done "<< endl<<endl;
        //
        // close file
        //
        f.Close();      
        //
        // plot histos
        //
        }  
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			}