class TCandList;
class TCandidate;
class TFitParams;
class PndAnaPidSelector;
class PndAnaPidCombiner;

void run_track_his(TString base="trackcheck", int nevts=10000, Int_t pid= 11, Float_t p = 0.3, Float_t th = 14)
{


	gStyle->SetOptFit(1011);

	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
	FairLogger::GetLogger()->SetLogToFile(kFALSE);




	// *** some variables
	int i=0,j=0, k=0, l=0;
	

        // Input file  
        TString identifier;  
        if(pid  > 0){    
        if(th>=100) identifier = Form("/data_pos_%d_%2.1f_%3.0f",pid,p,th);
        else identifier = Form("/data_pos_%d_%2.1f_%2.0f",pid,p,th);
        }
        else {
           pid= -1*pid;      
           if(th>=100) identifier = Form("/data_neg_%d_%2.1f_%3.0f",pid,p,th);
           else identifier = Form("/data_neg_%d_%2.1f_%2.0f",pid,p,th);  
        }
        base=base+identifier.Data();
        base="$HOME/GSI/pandaroot_trunk/macro/online/data/"+base;
        cout<< base <<endl;
	
        TString inDigiFile  = base+"_dig.root";
        TString inSimFile   = base+"_sim.root";
        TString inRecoFile  = base+"_rec.root";
        TString inPidFile   = base+"_pid.root";
        TString inParFile   = base+"_par.root";
        TString dummyFile   = base+"_dum.root";
        TString OutFile     = base+"_his.root";

	
	
	// *** initialization
	FairRunAna* fRun = new FairRunAna();
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
	fRun->SetInputFile(inSimFile);
	fRun->AddFriend(inPidFile);
	fRun->AddFriend(inRecoFile);

	FairParRootFileIo* parIO = new FairParRootFileIo();
	parIO->open(inParFile);
	rtdb->setFirstInput(parIO);
	rtdb->setOutput(parIO);  
	
	fRun->SetOutputFile(dummyFile);
	fRun->Init(); 


        /////////////////////////////////////////////////////
        TFile *inFile = TFile::Open(inRecoFile,"READ");  
        TTree *tree=(TTree *) inFile->Get("cbmsim") ;

        //TClonesArray* cand_array=new TClonesArray("PndPidCandidate");
        //tree->SetBranchAddress("PidChargedCand", &cand_array);
        //TClonesArray* drc_array=new TClonesArray("PndPidProbability");
        //tree->SetBranchAddress("PidAlgoDrc", &drc_array);
        //TClonesArray* mvd_array=new TClonesArray("PndPidProbability");
        //tree->SetBranchAddress("PidAlgoMvd", &mvd_array);
	TClonesArray* mc_array=new TClonesArray("PndMCTrack");
        tree->SetBranchAddress("MCTrack",&mc_array);
	FairMCEventHeader* event;
        tree->SetBranchAddress("MCEventHeader.", &event);



        PndPidProbability *flux = new PndPidProbability(239+237,114+101,2282+2375,35+42,517+1052);
	
	// *** create an output file for all histograms
	TFile *out = TFile::Open(OutFile,"RECREATE");
	

	
	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis();
	if (nevts==0) nevts= theAnalysis->GetEntries();
	
	// *** TCandLists for the analysis
	//TCandList all, chrg, mctrk, el, eplus, eminus, piplus, piminus, jpsi, psi2s;
	
	// *** Mass selector for the jpsi cands
	TPidMassSelector *jpsiMassSel=new TPidMassSelector("jpsi",3.096,1.0);
	
	// *** the MC truth matcher
	PndMcTruthMatch mcm;
	
	// *** the lorentz vector of the initial psi(2S)
	TLorentzVector ini(0, 0, 6.231552, 7.240065);
	
/*
        PndAnaPidCombiner pidComb("PidCombiner");
        pidComb.SetTcaNames("PidAlgoDrc;PidAlgoMvd;PidAlgoDisc;PidAlgoMdtHardCuts;PidAlgoStt;PidAlgoEmcBayes");
*/

	// *** create some histograms

        TH1F *h0= new TH1F("h0","h0", 5000,  -2500,  2500);;
        TH1F *h10= new TH1F("h10","h10", 20,  0,  20);



	const int n = 5;	
        TH1F *h1[n];
	TH1F *h2[n];
	double binp1=0, binp2=0;
	double bint1=0, bint2=0;
	if(pid==11){
        if(p >= 0.2 && p < 0.4) {binp1=-0.1;  binp2=0.25;  bint1=-2.5; bint2=2.5;}
	if(p >= 0.9 && p < 1.1) {binp1=-0.15; binp2=0.25; bint1=-1.0; bint2=1.0;}
        if(p >= 1.4 && p < 1.6) {binp1=-0.2;  binp2=0.3;  bint1=-0.6; bint2=0.6;}
        if(p >= 1.9 && p < 2.1) {binp1=-0.2;  binp2=0.3;  bint1=-0.5; bint2=0.5;}
        if(p >= 2.9 && p < 3.1) {binp1=-0.2;  binp2=0.3;  bint1=-0.4; bint2=0.4;}
        if(p >= 3.9 && p < 4.1) {binp1=-0.2;  binp2=0.3;  bint1=-0.3; bint2=0.3;}
        if(p >= 4.9 && p < 5.1) {binp1=-0.2;  binp2=0.3;  bint1=-0.2; bint2=0.2;}
        }
	else{	
	if(p >= 0.2 && p < 0.4) {binp1=-0.1;  binp2=0.1;  bint1=-2.5; bint2=2.5;}
	if(p >= 0.9 && p < 1.1) {binp1=-0.15; binp2=0.15; bint1=-1.0; bint2=1.0;}
        if(p >= 1.4 && p < 1.6) {binp1=-0.2;  binp2=0.2;  bint1=-0.6; bint2=0.6;}
        if(p >= 1.9 && p < 2.1) {binp1=-0.2;  binp2=0.2;  bint1=-0.5; bint2=0.5;}
        if(p >= 2.9 && p < 3.1) {binp1=-0.2;  binp2=0.2;  bint1=-0.4; bint2=0.4;}
        if(p >= 3.9 && p < 4.1) {binp1=-0.2;  binp2=0.2;  bint1=-0.3; bint2=0.3;}
        if(p >= 4.9 && p < 5.1) {binp1=-0.2;  binp2=0.2;  bint1=-0.2; bint2=0.2;}	
	}
	
	h1[0]= new TH1F("h1[0]","h1[0]", 100,  binp1,  binp2);
	h1[1]= new TH1F("h1[1]","h1[1]", 100,  binp1,  binp2);
	h1[2]= new TH1F("h1[2]","h1[2]", 100,  binp1,  binp2);
	h1[3]= new TH1F("h1[3]","h1[3]", 100,  binp1,  binp2);
	h1[4]= new TH1F("h1[4]","h1[4]", 100,  binp1,  binp2);

        h2[0]= new TH1F("h2[0]","h2[0]", 100,  bint1,  bint2);
        h2[1]= new TH1F("h2[1]","h2[1]", 100,  bint1,  bint2);
        h2[2]= new TH1F("h2[2]","h2[2]", 100,  bint1,  bint2);
        h2[3]= new TH1F("h2[3]","h2[3]", 100,  bint1,  bint2);
        h2[4]= new TH1F("h2[4]","h2[4]", 100,  bint1,  bint2);







/*

        PndAnaPidSelector positSel("PositSelector");	
        positSel.SetSelection("ElectronVeryLoosePlus");

        PndAnaPidSelector electSel("ElectSelector");	
        electSel.SetSelection("ElectronVeryLooseMinus");

        PndAnaPidSelector pimSel("PimSelector");	
        pimSel.SetSelection("PionVeryLooseMinus");

        PndAnaPidSelector pipSel("PipSelector");	
        pipSel.SetSelection("PionVeryLoosePlus");

*/

        TCandList charged, pipLoose, pimLoose, electLoose, positLoose;
 	TCandList all, chrg, mctrk, eplus, eminus, muplus, muminus, piplus, piminus, kplus, kminus, pplus, pminus;

	// ***
	// the event loop
	// ***		
	
	int counter1 = 0;
	int counter2 = 0;
	int counter3 = 0;
	int counter4 = 0;
	int counter5 = 0;
 
        cout<<"----------------------------------------------------------"<<theAnalysis->GetEntries()<<endl; 
	cout<<"----------------------------------------------------------"<<tree->GetEntriesFast()<<endl;
 		 
	while (theAnalysis->GetEvent() && i++<tree->GetEntriesFast())
	{
		
		tree->GetEntry(i-1);	    
		
		
	
		// *** the MC Truth objects
		theAnalysis->FillList(mctrk,"McTruth");
		h0->Fill(mctrk[0].PdgCode());
		theAnalysis->FillList(charged,"Charged");
		
		
		// *** Select with no PID info ('All'); type and mass are set 		
		theAnalysis->FillList(eplus,   "ElectronAllPlus");
		theAnalysis->FillList(eminus,  "ElectronAllMinus");
		theAnalysis->FillList(muplus,  "MuonAllPlus");
		theAnalysis->FillList(muminus, "MuonAllMinus");
		theAnalysis->FillList(piplus,  "PionAllPlus");
		theAnalysis->FillList(piminus, "PionAllMinus");
		theAnalysis->FillList(kplus,   "KaonAllPlus");
		theAnalysis->FillList(kminus,  "KaonAllMinus");
		theAnalysis->FillList(pplus,   "ProtonAllPlus");
		theAnalysis->FillList(pminus,  "ProtonAllMinus");

		//pidComb.Apply(charged);
		//pipSel.Select(charged, pipLoose);
		//pimSel.Select(charged, pimLoose);
		//electSel.Select(charged, electLoose);
		//positSel.Select(charged, positLoose);

		//cout<<mctrk.GetLength()<<endl;
		//		cout<<mctrk[0].P()<<endl;
		//		cout<<mctrk[0].PdgCode() <<endl;
		//for (j=0;j<mctrk.GetLength();++j) 
		//{
		//cout<<mctrk[j].PdgCode() <<endl;
		//}
		
		int number1=0;
		int number2=0;
		int number3=0;
		int number4=0;
		int number5=0;		 
                ///////////// electron - 
		for (j=0;j<eminus.GetLength();++j) 
		{			 
                      Int_t mcIdx = eminus[j].GetMcIdx();		      
		      if ( mcIdx >=0 && mcIdx < mctrk.GetLength())
		      {		       		
		         TCandidate truth = mctrk[mcIdx];
		         if(truth.PdgCode() == eminus[j].PdgCode())
		         {			   
			    TVector3 gen_p = truth.P3();  
			    TVector3 rec_p = eminus[j].P3();
			    			    
			    double deltap = truth.P()  -  eminus[j].P();
                            double deltat = gen_p.Theta()*TMath::RadToDeg() - rec_p.Theta()*TMath::RadToDeg();

			    h1[0]->Fill(deltap/truth.P());
			    h2[0]->Fill(deltat);

			    number1++;
 
                         } 
		      }
		      
		      
                } 

                ///////////// muon - 
		for (j=0;j<muminus.GetLength();++j) 
		{			 
                      Int_t mcIdx = muminus[j].GetMcIdx();		      
		      if ( mcIdx >=0 && mcIdx < mctrk.GetLength())
		      {		       		
		         TCandidate truth = mctrk[mcIdx];
		         if(truth.PdgCode() == muminus[j].PdgCode())
		         {			    
			    TVector3 gen_p = truth.P3();  
			    TVector3 rec_p = muminus[j].P3();
			    			    
			    double deltap = truth.P()  -  muminus[j].P();
                            double deltat = gen_p.Theta()*TMath::RadToDeg() - rec_p.Theta()*TMath::RadToDeg();

			    h1[1]->Fill(deltap/truth.P());
			    h2[1]->Fill(deltat);
			    number2++;

                         } 
		    }
                }                


                ///////////// pi + 
		for (j=0;j<piplus.GetLength();++j) 
		{			 
                       
                      Int_t mcIdx = piplus[j].GetMcIdx();		      
		      if ( mcIdx >=0 && mcIdx < mctrk.GetLength())
		      {		       		
		         TCandidate truth = mctrk[mcIdx];
		         if(truth.PdgCode() == piplus[j].PdgCode())
		         {
			    TVector3 gen_p = truth.P3();  
			    TVector3 rec_p = piplus[j].P3();
			    
			    double deltap = truth.P()  -  piplus[j].P();
                            double deltat = gen_p.Theta()*TMath::RadToDeg() - rec_p.Theta()*TMath::RadToDeg();

			    h1[2]->Fill(deltap/truth.P());
			    h2[2]->Fill(deltat);
			    
			    number3++;

                         } 
		    }
                }                


                ///////////// K +
		for (j=0;j<kplus.GetLength();++j) 
		{			                        
                      Int_t mcIdx = kplus[j].GetMcIdx();		      
		      if ( mcIdx >=0 && mcIdx < mctrk.GetLength())
		      {		       		
		         TCandidate truth = mctrk[mcIdx];
		         if(truth.PdgCode() == kplus[j].PdgCode())
		         {
			    TVector3 gen_p = truth.P3();  
			    TVector3 rec_p = kplus[j].P3();
			    
			    double deltap = truth.P()  -  kplus[j].P();
                            double deltat = gen_p.Theta()*TMath::RadToDeg() - rec_p.Theta()*TMath::RadToDeg();
			    h1[3]->Fill(deltap/truth.P());
			    h2[3]->Fill(deltat);			    
			    number4++;

                         } 
		    }
                }                


                ///////////// proton
		for (j=0;j<pplus.GetLength();++j) 
		{			   
                      Int_t mcIdx = pplus[j].GetMcIdx();		      
		      if ( mcIdx >=0 && mcIdx < mctrk.GetLength())
		      {		       		
		         TCandidate truth = mctrk[mcIdx];
		         if(truth.PdgCode() == pplus[j].PdgCode())
		         {
			    TVector3 gen_p = truth.P3();  
			    TVector3 rec_p = pplus[j].P3();
			    
			    double deltap = truth.P()  -  pplus[j].P();
                            double deltat = gen_p.Theta()*TMath::RadToDeg() - rec_p.Theta()*TMath::RadToDeg();
			    h1[4]->Fill(deltap/truth.P());
			    h2[4]->Fill(deltat);			    
			    number5++;

                         } 
		    }
                }                

   

                if( number1 > 0 ) {counter1++;   h10->Fill(0,1);}
                if( number2 > 0 ) {counter2++;   h10->Fill(1,1);}
                if( number3 > 0 ) {counter3++;   h10->Fill(2,1);}
                if( number4 > 0 ) {counter4++;   h10->Fill(3,1);}
                if( number5 > 0 ) {counter5++;   h10->Fill(4,1);}



 	}
	
	cout<<"------------------------------------------"<<endl;

       
	
	
	// *** write out all the histos
	out->cd();
	h0->Write();	
	h10->Write();
	for (j=0;j<n;++j)
	{ 
	  h1[j]->Write();
	  h2[j]->Write();


        }

	out->Save();
	
}