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