#include <CaloProfile.h>

////////////////////////////////////////////////////////////////////////            
/**
 * 1-dimension function describing lateral distribution of the 
 * shower as viewed by calorimeter
 * (projection of 3d function in one direction)
 *  
 *  xi[0] = x or y coordinate relative to shower axis
 *  parmin[0] = rt
 *  parmin[1] = p
 *  parmin[2] = rc
 *
 */
////////////////////////////////////////////////////////////////////////            
Double_t cfradx(Double_t *xi, Double_t *parmin) {

        double fradxmin2,p,rt,rc,es,x,pig,norm,c;
        x=*xi;
        pig = acos(-1.);
        rt=parmin[0];
        p=parmin[1];
        rc=parmin[2];
        norm=parmin[3];
        c=parmin[4];
        x=x-c;
        es=1.5;
        fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es);
        fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es);
        fradxmin2=norm*fradxmin2/(2*pig);     
        //cout<<"x,fradxmin2 "<< x<<" "<<fradxmin2  <<endl;
  return fradxmin2;
}


Double_t cfunc(Double_t *tin,Double_t *par)
{

 Double_t gaa,a,tm,et,value,t,t0;
 t=*tin;
 et=par[0];
 a=par[1];
 tm=par[2];
 t0=par[3];
 gaa=TMath::Gamma(a);
 
 value=et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm);
 return value;
}


//--------------------------------------
/**
 * Default constructor 
 */
CaloLat::CaloLat(){
  Clear();
};

/**
 * Default constructor 
 */
CaloLong::CaloLong(){
  Clear();
};

CaloLat::CaloLat(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  debug = false;
  //
};

CaloLong::CaloLong(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  debug = false;
  //
};

void CaloLat::Clear(){
  //
};

void CaloLong::Clear(){
  //
};

void CaloLat::Print(){
  //
  Process();
  //
  printf("==================== Calorimeter Lateral Profile =======================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
//  printf(" nx [number of X combination]:.. %i\n",nx);
//  printf(" ny [number of Y combination]:.. %i\n",ny);
  printf("========================================================================\n");
  //
};

void CaloLong::Print(){
  //
  Process();
  //
  printf("==================== Calorimeter Lateral Profile =======================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
//  printf(" nx [number of X combination]:.. %i\n",nx);
//  printf(" ny [number of Y combination]:.. %i\n",ny);
  printf("========================================================================\n");
  //
};

void CaloLat::Draw(){
  //
  Process();
  Draw(-1,-1);
};

void CaloLong::Draw(){
  //
  Process();
  Draw(-1);
};


void CaloLat::Draw(Int_t view,Int_t plane){
  //
  Int_t minv = 0;
  Int_t maxv = 0;
  Int_t minp = 0;
  Int_t maxp = 0;
  //
  if ( view == -1 ){
	minv = 0;
	maxv = 2;
  } else {
	minv = view;
	maxv = view+1;
  };

  if ( plane == -1 ){
	minp = 0;
	maxp = 22;
  } else {
	minp = plane;
	maxp = plane+1;
  };
  // 
  Process();
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  for (Int_t v=minv; v<maxv;v++){
     for (Int_t p=minp; p<maxp;p++){
	TString hid = Form("clatv%ip%i",v,p);	
        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
	if ( tc ){
//	 tc->Clear();
	} else {
	 tc = new TCanvas(hid,hid);
	};
	//
	TString thid = Form("hlatv%ip%i",v,p);	
        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
	if ( th ) th->Delete();
//	 th->Clear();
//	 th->Reset();
//	} else {
	 th = new TH1F(thid,thid,96,-0.5,95.5);
//	};
	tc->cd();
	//
	for (Int_t st=0;st<96;st++){
	 th->Fill(st,estrip[v][p][st]);
	};
	th->Draw();
	tc->Modified();
	tc->Update();
     };
  };
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void CaloLong::Draw(Int_t view){
  //
  Int_t minv = 0;
  Int_t maxv = 0;
  //
  if ( view == -1 ){
	maxv = -1;
  } else {
	minv = view;
	maxv = view+1;
  };
  // 
  Process();
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  if ( maxv != -1 ){
  for (Int_t v=minv; v<maxv;v++){
	TString hid = Form("clongv%i",v);	
        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
	if ( tc ){
//	 tc->Clear();
	} else {
	 tc = new TCanvas(hid,hid);
	};
	//
	TString thid = Form("hlongv%i",v);	
        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
	if ( th ) th->Delete();
//	 th->Clear();
//	 th->Reset();
//	} else {
	 th = new TH1F(thid,thid,22,-0.5,21.5);
//	};
	tc->cd();
	//
	for (Int_t st=0;st<22;st++){
	 th->Fill(st,eplane[v][st]);
	};
	th->Draw();
	tc->Modified();
	tc->Update();
  };
  } else {
	//
	TString hid = Form("clongvyvx");	
        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
	if ( tc ){
	} else {
	 tc = new TCanvas(hid,hid);
	};
	//
	TString thid = Form("hlongvyvx");	
        TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
	if ( th ) th->Delete();
        th = new TH1F(thid,thid,44,-0.5,43.5);
	tc->cd();
        Int_t pp=0;
	for (Int_t st=0;st<22;st++){
	  for (Int_t v=1; v>=0;v--){
		//
	        th->Fill(pp,eplane[v][st]);
		//
		pp++;
	  };
	};
	th->Draw();
	tc->Modified();
	tc->Update();
  };
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void CaloLat::Delete(){
  Clear();
  //delete this;
};

void CaloLong::Delete(){
  Clear();
  //delete this;
};

void CaloLat::Process(){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloHough variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  // let's start
  //
  memset(estrip,0, 4224*sizeof(Float_t));
  Float_t mip1 = 0.;
  Int_t view1 = 0;
  Int_t plane1 = 0;
  Int_t strip1 = 0;
  //
  for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
    mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
    estrip[view1][plane1][strip1] = mip1;
  };
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
  //
};

void CaloLong::Process(){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloHough variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  // let's start
  //
  memset(eplane,0, 2*22*sizeof(Float_t));
  //
  for (Int_t v=0; v<2; v++){
   for (Int_t i=0; i<22; i++){
     eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i);
   };
  };
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
  //
};