//_____________________________________________________________________________
// example of photon analysis module by Ray Culbertson
// one can do in the interactive session
//_____________________________________________________________________________

/*

  ****** to run: 
  (note naming convention: class_name = "T"+module_name+"Module"), 
  the source code of the class is stored in 

  .X photon_ana.C
//  .X temp.C
  .L TPhotonAnaModule.cc+
  //TStnAna x("stn.root");
  TStnAna x("/cdf/data04/s3/commission/1x8/stntuple/111215_stnmaker.root");
  //TStnAna x("/cdf/data02/s6/exotic/rlc/pythia_dipho_stn.root");
  x.AddModule(new TPhotonAnaModule);
  x.Run();

  Temp_strip->Draw(); Pho_convmiss->Draw(); Pho_convr->Draw();

  //   to display a histogram

  TPhotonAnaModule* m  = (TPhotonAnaModule*) x.GetModule("PhotonAna");
  m->SavePlots(7);


  //TCesAnaModule::Hist_t* c;
  //c = m->GetHist();
  //c->fAdcCounts->Draw();

  ****** to process one event (#11 in the ntuple) and to print the data

  x.ProcessEvent(11)
  TPhotonAnaModule* m = (TPhotonAnaModule*) x.GetModule("PhotonAna");
  TStnClusterBlock*   clu = m->GetClusterBlock();
  clu->Print();

*/

#include <iostream>
#include "TF1.h"
#include "TCanvas.h"
#include "Stntuple/loop/TStnAna.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/obj/TStnJetBlock.hh"
#include "Stntuple/obj/TStnMetBlock.hh"

#include "TPhotonAnaModule.hh"

// ClassImp(TPhotonAnaModule)
//_____________________________________________________________________________
TPhotonAnaModule::TPhotonAnaModule(const char* name, const char* title):
  TStnModule(name,title)
{
}

//_____________________________________________________________________________
TPhotonAnaModule::~TPhotonAnaModule() {
}

//_____________________________________________________________________________
void TPhotonAnaModule::BookHistograms() {
  //char name [200];
  //char title[200];

  //GetListOfHistograms()->Delete();
  //sprintf(name, "%s_ces_wire_adc",GetName());
  //sprintf(title,"%s: Cluster: Ces Wire ADC",GetName());
  //fHist.fCesWireAdc= new TH1F(name, title,500 ,0.001,500);
  //AddHistogram(fHist.fCesWireAdc);

  fHist.fNpho = new TH1F("Pho_npho","N Pho",11,-0.5,10.5);
  fHist.fDet = new TH1F("Pho_det","Det",2,-0.5,1.5);
  fHist.fEt = new TH1F("Pho_et","Et",50,0,100.0);
  fHist.fCo4 = new TH1F("Pho_co4","Co4",50,0.0,50.0);
  fHist.fCo4PJW = new TH1F("Pho_co4PJW","Co4 PJW",50,0.0,50.0);
  fHist.fCesx = new TH1F("Pho_cesx","Cesx",50,-28.0,28.0);
  fHist.fCesz = new TH1F("Pho_cesz","Cesz",50,-280.0,280.0);
  fHist.fCpr5ps = new TH1F("Pho_cpr5ps","Cpr5ps",50,-28.,28.);
  fHist.fCpr5ph = new TH1F("Pho_cpr5ph","Cpr5ph",50,0.,30000.);
  fHist.fDteta = new TH1F("Pho_Dteta","Dteta",50,-4.0,4.0);
  fHist.fCese = new TH1F("Pho_Cese","Ces E",50,0.,150.);
  fHist.fCeseRat = new TH1F("Pho_ceseRat","Cese/CEM",50,0.,3.0);
  fHist.fPhi = new TH1F("Pho_phi","Phi",50,0.,6.29);
  fHist.fHadem = new TH1F("Pho_hadem","Hadem",50,0.,1.0);
  fHist.fLshr = new TH1F("Pho_lshr","Lshr",50,-3.5,3.0);
  fHist.fSumpt4 = new TH1F("Pho_sumpt4","Sumpt4",50,0.,50.);
  fHist.fN3d = new TH1F("Pho_n3d","N3d",11,-0.5,10.5);
  fHist.fChi = new TH1F("Pho_chi","Chi",50,0.,30.);
  fHist.fCescprx = new TH1F("Pho_cescprx","Ces CPR x",50,-28.,28.);
  fHist.fCescprq = new TH1F("Pho_cescprq","Ces CPR q",50,0.,30000.);
  fHist.fPt = new TH1F("Pho_pt","Pt",50,0.,50.);
  //  fHist.f = new TH1F("Pho_det","Det",2,-0.5,1.5);
  fHist.fRho = new TH1F("Pho_rho","Rho",50,0.,3.0);
  fHist.fDipho = new TH1F("Pho_dipho","Dipho",100,0.,100.0);

  // for Steve
  fHist.fCesx2 = new TH1F("Pho_cesx2","Cesx e>2",50,-28.0,28.0);
  fHist.fCesz2 = new TH1F("Pho_cesz2","Cesz e>2",50,-280.0,280.0);
  fHist.fCeseRat2 = new TH1F("Pho_ceseRat2","Cese/CEM isoc",50,0.,3.0);
  fHist.fCeseRat3 = new TH1F("Pho_ceseRat3","Cese/CEM isoc,Et>10",50,0.,3.0);

  fHist.fCesCprDiff = new TH1F("Pho_cescprdiff","Ces-Cpr",50,-28.0,28.0);

  fHist.fPesE   = new TH1F("Pho_pese",  "Pes E",50,0.0,8000.0);
  fHist.fPesX   = new TH1F("Pho_pesx",  "Pes x",50,-180.0,180.0);
  fHist.fPesY   = new TH1F("Pho_pesy",  "Pes y",50,-180.0,180.0);
  fHist.fPesPhi = new TH1F("Pho_pesphi","Pes phi",50,0.0,6.29);
  fHist.fPesR   = new TH1F("Pho_pesxr", "Pes R",50,0.0,180.0);

  fHist.fStrip = new TProfile("Temp_strip","Strip frac",11,-0.5,10.5);

  fHist.fConvMiss  = new TH1F("Pho_convmiss", "Conv miss",50,0.0,10.0);
  fHist.fConvR  = new TH1F("Pho_convr", "Conv R",50,0.0,100.0);
  fHist.fConvMissWS  = new TH1F("Pho_convmissWS", "Conv miss WS",50,0.0,10.0);
  fHist.fConvRWS  = new TH1F("Pho_convrWS", "Conv R WS",50,0.0,100.0);
}


//_____________________________________________________________________________
int TPhotonAnaModule::BeginJob() {

				// register the data block

  RegisterDataBlock("CalDataBlock","TCalDataBlock",&fCalDataBlock);
  RegisterDataBlock("CesDataBlock","TCesDataBlock",&fCesDataBlock);
  RegisterDataBlock("CprDataBlock","TCprDataBlock",&fCprDataBlock);
  //  fJetBlock =     (TStnJetBlock*) 
  //        RegisterDataBlock("JetBlock","TStnJetBlock");

  RegisterDataBlock("MetBlock"    ,"TStnMetBlock",&fMetBlock    );
  RegisterDataBlock("TrackBlock"  ,"TStnTrackBlock",&fTrackBlock  );
  RegisterDataBlock("ClusterBlock","TStnClusterBlock",&fClusterBlock);
  RegisterDataBlock("PhotonBlock" ,"TStnPhotonBlock",&fPhotonBlock );

					// book histograms
  BookHistograms();

  return 0;
}


//_____________________________________________________________________________
int TPhotonAnaModule::BeginRun() {
  return 0;
}

//_____________________________________________________________________________
int TPhotonAnaModule::Event(int ientry) {

  if (fCalDataBlock) fCalDataBlock->GetEntry(ientry);
  if (fCesDataBlock) fCesDataBlock->GetEntry(ientry);
  if (fCprDataBlock) fCprDataBlock->GetEntry(ientry);

  //fJetBlock->GetEntry(ientry);
  fTrackBlock->GetEntry(ientry);
  fClusterBlock->GetEntry(ientry);

  if (! fPhotonBlock) return -1;

  fPhotonBlock-> GetEntry(ientry);
  fMetBlock->GetEntry(ientry);

  int npho = fPhotonBlock->NPhotons();

  //if(npho>1) fPhotonBlock->Print();

  fHist.fNpho->Fill(npho);
  for (int i=0; i<npho; i++) {
    TStnPhoton& pho = *(fPhotonBlock->Photon(i));
    fHist.fDet->Fill(pho.fDetector);
    fHist.fEt->Fill(pho.fEt);
    fHist.fCo4->Fill(pho.fCo4);
    fHist.fCo4PJW->Fill(pho.fCo4PJW);
    fHist.fCesx->Fill(pho.fCesx);
    fHist.fCesz->Fill(pho.fCesz);
    fHist.fCpr5ph->Fill(pho.fCpr5ph);
    fHist.fCpr5ps->Fill(pho.fCpr5ps);
    fHist.fDteta->Fill(pho.fDteta);
    fHist.fCese->Fill(pho.fCese);
    float rat = (pho.fE>0 ? pho.fCese/pho.fE : -1.0);
    fHist.fCeseRat->Fill(rat);
    fHist.fPhi->Fill(pho.fPhi);
    fHist.fHadem->Fill(pho.fHadem);
    fHist.fLshr->Fill(pho.fLshr);
    fHist.fSumpt4->Fill(pho.fSumpt4);
    fHist.fN3d->Fill(pho.fN3d);
    fHist.fChi->Fill(pho.fChi);
    fHist.fCescprx->Fill(pho.fCescprx);
    fHist.fCescprq->Fill(pho.fCescprq);
    fHist.fPt->Fill(pho.fPt);
    //   fHist.fDet->Fill(pho.fEt);

    // for Steve
    if(pho.fCo4PJW<2.0) fHist.fCeseRat2->Fill(rat);
    if(pho.fCo4PJW<2.0 && pho.fEt>10.0) fHist.fCeseRat3->Fill(rat);
    if(pho.fCese>2.0) {
      fHist.fCesx2->Fill(pho.fCesx);
      fHist.fCesz2->Fill(pho.fCesz);
    }
    if(pho.fCese>0.1 && pho.fCescprq>1000.0) 
      fHist.fCesCprDiff->Fill(-pho.fCesx-pho.fCescprx);

    fHist.fPesE->Fill(pho.fPesE);
    if(pho.fPesE>0.0) {
      fHist.fPesX->Fill(pho.fPesX);
      fHist.fPesY->Fill(pho.fPesY);
      float phi = atan2(-pho.fPesY,-pho.fPesX)+M_PI;
      fHist.fPesPhi->Fill(phi);
      float r = sqrt(pho.fPesX*pho.fPesX+pho.fPesY*pho.fPesY);
      fHist.fPesR->Fill(r);
    }
  }

  //RhoAna();
  //PhysAna();
  SetFilter();

return 0;
}


//_____________________________________________________________________________
void TPhotonAnaModule::RhoAna() {

  //std::cout << " begin rho" << std::endl;
  int ntrk = fTrackBlock->NTracks();
  int npho = fPhotonBlock->NPhotons();

  //std::cout << fPhotonBlock << " " << fTrackBlock << std::endl;
  TLorentzVector tvec,pvec;
  for (int i=0; i<npho; i++) {
    TStnPhoton& pho = *(fPhotonBlock->Photon(i));
    //if(pho.fCo4PJW<2.0 && pho.fChi<10.0 && pho.fEt>20.0) {
    if(pho.fCo4PJW<2.0 && pho.fDetector==0  && pho.fEt>10.0) {


      for(int j=0; j<ntrk; j++) {
	TStnTrack& trk = *(fTrackBlock->Track(j));
	if(trk.Pt()>1.5) {
	//	if(trk.NAxSeg()>2 && trk.NStSeg()>2) {
	//std::cout << i << " " << j << std::endl;
	//std::cout << &pho<< " " << &trk << std::endl;
	//std::cout << fTrackBlock->NTracks() << std::endl;

	/*
	  */
	  //	  double phi = trk.fPhi0;
	  //  double pt  = trk.Pt();
	  double z   = trk.Z0();

	  tvec = *trk.Momentum();  // 0.139570
	  //tvec.SetPtEtaPhiM(trk.Pt(),trk.fEveta,trk.fPhi,trk.fE);
	  double ev = pho.fEveta;
	  double cot = (TMath::Exp(ev)-TMath::Exp(-ev))/2.0;
	  double z0p = cot*184.0;
	  z0p = z0p - z;
	  double pr = TMath::Sqrt(z0p*z0p + 184*184);
	  double pmag = TMath::Sqrt(fabs(pho.fE*pho.fE - 0.134976*0.134976));
	  double pt = pmag/TMath::Sqrt(1+cot*cot);
	  double px = pt*TMath::Cos(pho.fPhi);
	  double py = pt*TMath::Sin(pho.fPhi);
	  double pz = TMath::Sqrt(fabs(pmag*pmag-pt*pt));

	  pvec= TLorentzVector(px,py,pz,pho.fE);
	  //pvec.SetPtEtaPhiM(pho.fEt,pho.fEveta,pho.fPhi,pho.fE);

	  //pvec.SetVectM(pho.fMomentum.Vect(),0.134976);
	  float mass = sqrt(fabs(tvec*pvec));


	  //	float mass=0;
	  //std::cout << fHist.fRho << std::endl;

	  fHist.fRho->Fill(mass);

      } // end trk if
      } // end trk loop
    } // end pho if
  } // end pho loop
  //std::cout << " done rho" << std::endl;

}

//_____________________________________________________________________________
void TPhotonAnaModule::PhysAna() {

  int npho = fPhotonBlock->NPhotons();

  //std::cout << fPhotonBlock << " " << fTrackBlock << std::endl;
  TLorentzVector pvec1,pvec2;
  for (int i=0; i<(npho-1); i++) {
    TStnPhoton& pho1 = *(fPhotonBlock->Photon(i));
    if(pho1.fCo4PJW<2.0 && pho1.fPt<1.0 && pho1.fDetector==0) {
      for (int j=i+1; j<npho; j++) {
        TStnPhoton& pho2 = *(fPhotonBlock->Photon(j));
        if(pho2.fCo4PJW<2.0 && pho2.fPt<1.0 && pho2.fDetector==0) {
	  //pvec1.SetVectM(pho1.fMomentum.Vect(),0.0);
	  //pvec2.SetVectM(pho2.fMomentum.Vect(),0.0);
	  pvec1.SetPtEtaPhiM(pho1.fEt,pho1.fEveta,pho1.fPhi,pho1.fE);
	  pvec2.SetPtEtaPhiM(pho2.fEt,pho2.fEveta,pho2.fPhi,pho2.fE);
	  float mass = sqrt(fabs(pvec1*pvec2));
	  //	  std::cout << mass << " " << pvec1.X() << " " << pvec2.X() << std::endl;
	  std::cout << mass << std::endl;
	  std::cout << pvec1.X() << " " << pvec1.Y() << " " << pvec1.Z() << " "<< pvec1.T() << std::endl;
	  std::cout << pho1.fEt << " " << pho1.fEveta << " " <<pho1.fPhi  << " "<< pho1.fE << std::endl;
	  std::cout << pvec2.X() << " " << pvec2.Y() << " " << pvec2.Z() << " "<< pvec2.T() << std::endl;
	  std::cout << pho2.fEt << " " << pho2.fEveta << " " <<pho2.fPhi  << " "<< pho2.fE << std::endl;


	  fHist.fDipho->Fill(mass);

      } // end pho2 if
      } // end pho2 loop
    } // end pho1 if
  } // end pho1 loop

}

void TPhotonAnaModule::SavePlots(int page) {

  if(page==0 || page==1) {
    TCanvas * c1= new TCanvas("Pho_page1","Pho Page 1",100,100,600,600);
    c1->Divide(2,3);
    c1->cd(1); fHist.fNpho->Draw();
    c1->cd(2); fHist.fDet->Draw();
    c1->cd(3); fHist.fEt->Draw();
    c1->cd(4); fHist.fDteta->Draw();
    c1->cd(5); fHist.fPhi->Draw();
    c1->cd(6); fHist.fHadem->Draw();
  }
  if(page==0 || page==2) {
    TCanvas * c2= new TCanvas("Pho_page2","Pho Page 2",150,150,600,600);
    c2->Divide(2,2);
    c2->cd(1); fHist.fCo4->Draw();
    c2->cd(2); fHist.fCo4PJW->Draw();
    c2->cd(3); fHist.fSumpt4->Draw();
    c2->cd(4); fHist.fPt->Draw();
  }
  if(page==0 || page==3) {
    TCanvas * c3= new TCanvas("Pho_page3","Pho Page 3",200,200,600,600);
    c3->Divide(2,3);
    c3->cd(1); fHist.fCese->Draw();
    c3->cd(2); fHist.fCesx->Draw();
    c3->cd(3); fHist.fCesz->Draw();
    c3->cd(4); fHist.fCeseRat->Draw();
    c3->cd(5); fHist.fChi->Draw();
  //  c3->cd(6); fHist.f->Draw();
  }
  if(page==0 || page==4) {
    TCanvas * c4= new TCanvas("Pho_page4","Pho Page 4",250,250,600,600);
    c4->Divide(2,3);
    c4->cd(1); fHist.fCpr5ps->Draw();
    c4->cd(2); fHist.fCpr5ph->Draw();
    c4->cd(3); fHist.fCescprx->Draw();
    c4->cd(4); fHist.fCescprq->Draw();
    c4->cd(5); fHist.fCesCprDiff->Draw();
  // c4->cd(6); fHist.f->Draw();
  }
  if(page==0 || page==5) {
    TCanvas * c5= new TCanvas("Pho_page5","Pho Page 5",300,300,600,600);
    c5->Divide(2,3);
    c5->cd(1); fHist.fRho->Draw();
    c5->cd(2); fHist.fDipho->Draw();
  }
  if(page==0 || page==6) {
    TCanvas * c6= new TCanvas("Pho_page6","Pho Page 6",350,350,600,600);
    c6->Divide(2,4);
    c6->cd(1); fHist.fCesx->Draw();
    c6->cd(2); fHist.fCesz->Draw();
    c6->cd(3); fHist.fCesx2->Draw();
    c6->cd(4); fHist.fCesz2->Draw();
    c6->cd(5); fHist.fCeseRat->Draw();
    c6->cd(6); fHist.fCeseRat2->Draw();
    c6->cd(7); fHist.fCeseRat3->Draw();
  }
  if(page==0 || page==7) {
    TCanvas * c7= new TCanvas("Pho_page7","Pho Page 7",400,400,600,600);
    c7->Divide(2,3);
    c7->cd(1); fHist.fPesE->Draw();
    c7->cd(2); fHist.fPesX->Draw();
    c7->cd(3); fHist.fPesY->Draw();
    c7->cd(4); fHist.fPesPhi->Draw();
    c7->cd(5); fHist.fPesR->Draw();
  }

}


//_____________________________________________________________________________
void TPhotonAnaModule::SetFilter() {

  bool pass = false;
  bool qconv = false;
  bool qconv2 = false;

  // npho is the same as the number of CdfEmobjects...
  int npho = fPhotonBlock->NPhotons();
  if(npho>0) { 
    fHist.fStrip->Fill(1,1);
    pass = true;
  } else {
    fHist.fStrip->Fill(1,0);
  }

  double ri,xi,yi,rj,xj,yj,ff,dd,sep,dt;
  bool rs;
				// coordinates of the conversion point

  double xc, yc, rc, dx, dy ,sqr;

  int ntrk = fTrackBlock->NTracks();
  for(int i=0; i<ntrk-1; i++) {
    TStnTrack& trki = *(fTrackBlock->Track(i));
    //  std::cout << trki.NAxSeg() << " " <<trki.NStSeg() << std::endl;
  for(int j=i+1; j<ntrk; j++) {
    TStnTrack& trkj = *(fTrackBlock->Track(j));
    dt = fabs(trki.Lam0()-trkj.Lam0());
    rs = (trki.Charge()*trkj.Charge()<0.0);
  if(dt<0.1) {
      //   && trki.NAxSeg()>1 && trki.NStSeg()>1 
      // && trkj.NAxSeg()>1 && trkj.NStSeg()>1) {
    ri = 1.0/fabs(2*trki.C0());
    ff = trki.Charge()*ri+trki.D0();
    xi = -ff*sin(trki.Phi0());
    yi =  ff*cos(trki.Phi0());
    rj = 1.0/fabs(2*trkj.C0());
    ff = trkj.Charge()*rj+trkj.D0();
    xj = -ff*sin(trkj.Phi0());
    yj =  ff*cos(trkj.Phi0());
    dd  = sqrt((xj-xi)*(xj-xi) + (yj-yi)*(yj-yi));
    sep = fabs(dd - ri - rj);
				// calculate radius of conversion. Define
				// it at a point where the tracks are parallel
				// to each other in XY, ignore track error
				// matrices for the time being
    dx  = xj-xi;
    dy  = yj-yi;
    sqr = sqrt(dx*dx+dy*dy);
    xc  = xi+dx/sqr*(ri+sep/2);
    yc  = yi+dy/sqr*(ri+sep/2);

    rc = sqrt(xc*xc+yc*yc);

    if(rs) {
      if(rc>28.0) fHist.fConvMiss->Fill(sep);
      if(sep<0.4)  fHist.fConvR->Fill(rc);
    } else {
      if(rc>28.0) fHist.fConvMissWS->Fill(sep);
      if(sep<0.4)  fHist.fConvRWS->Fill(rc);
    }

    if(sep<0.4 && rc>28.0) {
      qconv = true;
      if(trki.Pt()>1.0 && trkj.Pt()>1.0) qconv2 = true;
    }

  }
  }
  }

  if(qconv) { 
    fHist.fStrip->Fill(2,1);
    // pass = true;
  } else {
    fHist.fStrip->Fill(2,0);
  }
  
  if(qconv2) { 
    fHist.fStrip->Fill(3,1);
    pass = true;
  } else {
    fHist.fStrip->Fill(3,0);
  }
  

  if(pass) { 
    fHist.fStrip->Fill(0.,1.);
  } else {
    fHist.fStrip->Fill(0.,0.);
  }
  

  if(pass) SetPassed(1);

  return;
}


//_____________________________________________________________________________
int TPhotonAnaModule::EndJob() {
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}