//---------------This macro shows the result of the bumps after splitting in the shahslyk and the barrel emc.
// on the left you see the 3 largest bumps of an event in the barrel
// on the right the two largest in the shashlyk in c1.
//
//c2 shows the difference of two bumps in shashlyk in space and momentum

gROOT->Reset();
Double_t E = 0.3;
Double_t mp;
Double_t p;
Double_t mpi;
Double_t dmpi=-10;
Double_t EsumP;




void Fehlersuche(TString full, Double_t Momentum)
{
   gROOT->SetStyle("Plain");

        gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
        basiclibs();
        gSystem->Load("libGeoBase");
        gSystem->Load("libParBase");
        gSystem->Load("libBase");
        gSystem->Load("libMCStack");
        gSystem->Load("libField");
        gSystem->Load("libGen");
        gSystem->Load("libPassive");
       
        gSystem->Load("libEmc");
        gSystem->Load("libGeom.so");

     	TFile* f = new TFile(full);
	TTree *t=(TTree *) f->Get("cbmsim") ;

 	t->SetBranchStatus("*",0);
        t->SetBranchStatus("EmcBump.*",1);

	TClonesArray* cluster_array=new TClonesArray("PndEmcBump");
	t->SetBranchAddress("EmcBump",&cluster_array);


        Double_t threshold=0.02;
        mp=0.938271;
        p=Momentum;
	Double_t d=0;
	E=sqrt(2*mp*mp+2*mp*(sqrt(mp*mp+p*p)));
	TLorentzVector pvec(0,0,p,sqrt(p*p+mp*mp)+mp);
        TLorentzVector pvec2, pvec3;

        TLorentzVector gvec(pvec), gvec2(pvec),gvec3(pvec), gvec21(pvec),gvec31(pvec);


        TH1F *h2= new TH1F("h2","Momentum difference between 2. and 3. largest bump in Shashlyk",200,0,E);
        TH1F *h3= new TH1F("h3","Distance between 2. and 3. largest bump in Shashlyk",200,0,60);
        TH1F *h10= new TH1F("h10","Barrel",200,0,E);
        TH1F *h11= new TH1F("h11","Barrel",200,0,E);  
	TH1F *h12= new TH1F("h12","Barrel",200,0,E);
        TH1F *h20= new TH1F("h20","Shashlyk",200,0,E);
	TH1F *h21= new TH1F("h21","Shashlyk",200,0,E);
	TH1F *h22= new TH1F("h22","Shashlyk",200,0,E);
  
       
        Int_t i=0;
        Int_t k=0;
        const int n=80;
        Int_t cnr=0, Zsum;
        Double_t A[n]=0;
        Double_t B[n]=0;
	Double_t C[n]=0;


	Long64_t *idx = 0;
 	if (!(idx = new Long64_t[n])) {return 1;}
	Long64_t *idx1 = 0;
 	if (!(idx1 = new Long64_t[n])) {return 1;}
	Long64_t *idx2 = 0;
 	if (!(idx2 = new Long64_t[n])) {return 1;}
 
        TVector3 v1=0,v2,v3, boostV1;
	PndEmcBump *cluster;


        boostV1 = pvec.BoostVector();
 
        cout<<"Number of events : "<<t->GetEntriesFast()<<endl;
        for (Int_t j=0; j< t->GetEntriesFast(); j++)
        {
                cnr=0;
                t->GetEntry(j);
	
            
                for (i=0; i<cluster_array->GetEntriesFast(); i++)
                {
                        cluster=(PndEmcBump*)cluster_array->At(i);
	
                        if (cluster->energy()<threshold)continue;
                        v1=cluster->where();
			if(v1.Z()>200)
			{	v1.SetMag(cluster->energy());
                        	gvec=TLorentzVector(v1,cluster->energy());
                              	gvec.Boost(-boostV1);
                        	 A[i]=gvec.T();
			}
			else
			{ 	v1.SetMag(cluster->energy());
				gvec=TLorentzVector(v1,cluster->energy());
				gvec.Boost(-boostV1);
				B[i]=gvec.T();
			}
			C[i]=gvec.T();
                        cnr++;
			

                }
		TMath::Sort(n,B,idx2);
                TMath::Sort(n,A,idx);
		TMath::Sort(n,C,idx1);
	
		if(A[idx[0]]+A[idx[1]]!=0 )h20->Fill(A[idx[0]]+A[idx[1]]);
		if(A[idx[0]]!=0  )h21->Fill(A[idx[0]]);
		if(A[idx[1]]!=0 )h22->Fill(A[idx[1]]);
		if((B[idx2[0]]+B[idx2[1]])!=0 )h10->Fill(B[idx2[0]]+B[idx2[1]]);
		if(B[idx2[0]]!=0 )h11->Fill(B[idx2[0]]);
		if(B[idx2[1]]!=0 )h12->Fill(B[idx2[1]]);


		if(cnr>1){
// 				cluster=(PndEmcBump*)cluster_array->At(idx1[0]);
// 				v1=cluster->where();if(v1.Z()>200)
// 				v1.SetMag(cluster->energy());
// 				gvec=TLorentzVector(v1,cluster->energy());
// 				gvec.Boost(-boostV1);
		
				cluster=(PndEmcBump*)cluster_array->At(idx1[1]);
				v1=cluster->where();if(v1.Z()>200)Zsum++;
				
				gvec21=TLorentzVector(v1,cluster->energy());
				gvec21.Boost(-boostV1);

				v1.SetMag(cluster->energy());
				
				gvec2=TLorentzVector(v1,cluster->energy());
				gvec2.Boost(-boostV1);
				if(cnr>2){
					cluster=(PndEmcBump*)cluster_array->At(idx1[2]);
					v1=cluster->where();if(v1.Z()>200)Zsum++;

					gvec31=TLorentzVector(v1,cluster->energy());
					gvec31.Boost(-boostV1);

					v1.SetMag(cluster->energy());

					gvec3=TLorentzVector(v1,cluster->energy());
					gvec3.Boost(-boostV1);
	
				
				d=(gvec2.Vect()-gvec3.Vect()).Mag();d1=(gvec21.Vect()-gvec31.Vect()).Mag();
				if(Zsum=2){h2->Fill(d);h3->Fill(d1);}
				}
		}

 		for(i=0; i < n ; i++){
		 A[i]=0;B[i]=0;C[i]=0;

                }

	}



       TCanvas* c1 = new TCanvas("c1", "", 10, 10, 1250, 450);
	c1->Divide(2);
	c1->cd(1); 

        h10->SetStats(kTRUE);
        h10->GetXaxis()->SetTitle("sqrt(s)");
        h10->GetYaxis()->SetTitle("N");
        h10->SetFillColor(kBlue);
        h11->SetFillColor(kRed);

        h11->Draw();
	h10->Draw("same");
        h12->SetFillColor(kGreen);
        h12->Draw("same");


	c1->cd(2); 
        h20->SetStats(kTRUE);
        h20->GetXaxis()->SetTitle("sqrt(s)");
        h20->GetYaxis()->SetTitle("N");

        h20->SetFillColor(kBlue);
        h20->Draw();
        h21->SetFillColor(kRed);
        h21->Draw("same");
        h22->SetFillColor(kGreen);
        h22->Draw("same");

	TCanvas* c2 = new TCanvas("c2", "", 10, 10, 1250, 900);
	c2->Divide(1,2);
	c2->cd(1);
        h2->SetStats(kTRUE);
        h2->GetXaxis()->SetTitle("GeV");
        h2->GetYaxis()->SetTitle("N");
        h2->Draw();
	c2->cd(2);
        h3->SetStats(kTRUE);
        h3->GetXaxis()->SetTitle("GeV");
        h3->GetYaxis()->SetTitle("N");
        h3->Draw();

}