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.
//
//  --- Modified on January 2007 ---
//
//      see the comment at the beginning of example1.C      
//
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->GetPamTree(&f);         // << get Pamela trees from file f
    Int_t nevent = T->GetEntries();     
    // ********************
    // load magnetic field
    // ********************
    pam_event->GetTrkLevel2()->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->GetToFLevel2()->GetZTOF(pam_event->GetToFLevel2()->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++){
        //
        pam_event->Clear();

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