//	
//	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
	//
	}