HarpoAnalyseVertex3Dsim.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseVertex3Dsim.cxx
//
#include "HarpoAnalyseVertex3Dsim.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"
#include "HarpoSimEvent.h"
#include "TpcSimMCEvent.h"
#include "MakeNiceHisto.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TRootEmbeddedCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"
#include "TGLabel.h"
#include "TRandom.h"
#include "TSystem.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

ClassImp(HarpoAnalyseVertex3Dsim);

void   HarpoAnalyseVertex3Dsim::print()
 {
   unsigned int nd; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();

   assert(hdr != NULL);
   hdr->print();
   
   for (nd = 0; nd < gkNDetectors; nd++) {
     //     if (fEvt->isdataExist(nd)) {
       HarpoDccMap *plane = fEvt->GetDccMap(nd);
       if (plane != NULL )plane->print();
   }
 }

void   HarpoAnalyseVertex3Dsim::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  if(nEvents%1000==0)
    Info("process","Event %ld",nEvents);
  //???????????????????????????????
  if (!gHDetSet->isExist(SIMDET)) return;
  if(!fEvt->GetRecoEvent()) return;


  Double_t angle_open_sim = 0;
  //  Double_t angle_omega_sim = 0;
  TVector3 possim = TVector3(0,0,0);
  TVector3 pgammasim = TVector3(0,0,1);
  TVector3 ppolasim = TVector3(0,1,0);
  TVector3 pgamma = TVector3(0,0,1);
  Double_t egamma = -1;
  TVector3 ppola = TVector3(0,1,0);
  TVector3 ppositron, pelectron;
  Double_t epositron= -1.0, eelectron = -1.0;
  HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
  TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
  Int_t nMCtr = mcEvent->GetNtracks();
  //    if(nMCtr==2){
  //Double_t norm[2] = {0,0};
  Bool_t palready = kFALSE;
  for(Int_t i = 0; i<nMCtr; i++){
    TpcSimMCTrack* track = mcEvent->GetTrack(i);
    Int_t pdgcode = track->GetPdgCode();
    if(pdgcode == 11 && palready) pdgcode = -11;
    switch(pdgcode){
    case 22:
      possim = track->GetOrigin(); 
      pgammasim = track->GetP(); 
      ppolasim = track->GetPola(); 
      pgamma = track->GetP(); 
      egamma = track->GetE(); 
      ppola = track->GetPola(); 
      break;
    case 11: 
      palready = kTRUE;
      ppositron = track->GetP(); 
      epositron = track->GetE(); 
      break;
    case -11: 
      pelectron = track->GetP(); 
      eelectron = track->GetE(); 
      break;
    }
  }
  Double_t u1u2 = ppositron.Unit()*pelectron.Unit();
  angle_open_sim=acos(u1u2);
  

  // Info("process","Pgamma = (%g,%g,%g)",pgammasim.Px(),pgammasim.Py(),pgammasim.Pz());
  // Info("process","Ppola = (%g,%g,%g)",ppolasim.Px(),ppolasim.Py(),ppolasim.Pz());

  TVector3 ppola2 = ppolasim.Cross(pgammasim.Unit());
  TVector3 pPlusMC = TVector3(ppositron.Dot(pgammasim.Unit()),ppositron.Dot(ppolasim.Unit()),ppositron.Dot(ppola2.Unit())).Unit();
  TVector3 pMinusMC = TVector3(pelectron.Dot(pgammasim.Unit()),pelectron.Dot(ppolasim.Unit()),pelectron.Dot(ppola2.Unit())).Unit();

  Double_t phiPlusMC = TMath::ATan2(pPlusMC.Pz(),pPlusMC.Py());
  Double_t phiMinusMC = TMath::ATan2(pMinusMC.Pz(),pMinusMC.Py());
  Double_t phipmMC = 0.5*(phiPlusMC+phiMinusMC);
  if(gRandom->Rndm()<0.5) phipmMC += TMath::Pi();
  while(phipmMC>TMath::Pi()) phipmMC  -= 2*TMath::Pi();
  while(phipmMC<-TMath::Pi()) phipmMC += 2*TMath::Pi();
  
  hOpenMC->Fill(angle_open_sim);
  hPhiPmMC->Fill(phipmMC);

  HarpoRecoEvent* reco = fEvt->GetRecoEvent();

  TClonesArray* vArray = reco->GetVertex3DArray();
  //  TClonesArray* v2dArray = reco->GetVertexArray();
  Int_t nV = vArray->GetEntries();
  Int_t nIn = 0, nOut = 0;
  for(Int_t iV = 0; iV<nV; iV++){
    HarpoRecoVertex3D* vertex = (HarpoRecoVertex3D*)vArray->At(iV);

    Double_t e1 = eelectron; // TO DO FIX ME
    Double_t e2 = epositron; // TO DO FIX ME
    Double_t e  = egamma; // TO DO FIX ME

    TVector3 p1 = vertex->GetP1Vertex3D();
    TVector3 p2 = vertex->GetP2Vertex3D();
    // TLorentzVector pL1(p1,e1);
    // TLorentzVector pL2(p2,e2);
    Double_t xV = vertex->GetXvertex3D();
    Double_t yV = vertex->GetYvertex3D();
    Double_t zV = vertex->GetZvertex3D();
    // Int_t iVx = vertex->GetVertexX();
    // Int_t iVy = vertex->GetVertexY();
    TVector3 p = p1.Unit() + p2.Unit();
    //    TLorentzVector3 pL = pL1 + pL2;

    // HarpoRecoVertex* vertexX = (HarpoRecoVertex*)v2dArray->At(iVx);
    // HarpoRecoVertex* vertexY = (HarpoRecoVertex*)v2dArray->At(iVy);

    hDistRecoSim->Fill(xV-possim.x()-150,yV-possim.y()-150,zV-possim.z()-250);

    hPosition->Fill(xV,yV,zV);


    Double_t u1u2 = p1.Unit()*p2.Unit();
    Double_t angle_open=acos(u1u2);
    hOpen->Fill(angle_open);

    //    TVector3 ppola2 = ppola.Cross(pgamma);
    TVector3 pPlus = TVector3(p1.Dot(pgammasim.Unit()),p1.Dot(ppolasim.Unit()),p1.Dot(ppola2.Unit()));
    TVector3 pMinus = TVector3(p2.Dot(pgammasim.Unit()),p2.Dot(ppolasim.Unit()),p2.Dot(ppola2.Unit()));
    Double_t phiPlus = TMath::ATan2(pPlus.Pz(),pPlus.Py());
    Double_t phiMinus = TMath::ATan2(pMinus.Pz(),pMinus.Py());
    Double_t phipm = 0.5*(phiPlus+phiMinus);
    if(gRandom->Rndm()<0.5) phipm += TMath::Pi();
    while(phipm>TMath::Pi()) phipm  -= 2*TMath::Pi();
    while(phipm<-TMath::Pi()) phipm += 2*TMath::Pi(); 
    Double_t deltaphipm = phipm-phipmMC;
    while(deltaphipm>TMath::Pi()) deltaphipm  -= 2*TMath::Pi();
    while(deltaphipm<-TMath::Pi()) deltaphipm += 2*TMath::Pi(); 

    hPhiPm->Fill(phipm);
    Double_t theta = pgammasim.Angle(p1.Unit()+p2.Unit());
    //if(theta>TMath::Pi()/2) theta = TMath::Pi()-theta;
    hTheta->Fill(theta);

    Double_t data[19];
    for(Int_t i = 0; i<3; i++){
      data[0+i] = possim[i];
      data[3+i] = 0;
      data[6+i] = p[i];
      data[9+i] = pgammasim[i];
    }
    data[12] = e;
    data[13] = e1;
    data[14] = e2;
    data[15] = theta;
    data[16] = phipm;
    data[17] = angle_open;
    data[18] = 0;
    if(TMath::Abs(xV-possim.x()-150)<10 &&TMath::Abs(yV-possim.y()-150)<10 &&TMath::Abs(zV-possim.z()-250)<10){
      data[18] = 1;
      hPhiPmCut->Fill(phipm);
      hPositionCut->Fill(xV,yV,zV);
      hDeltaPhiPm->Fill(deltaphipm);
      hDeltaOpen->Fill(angle_open-angle_open_sim);

      hOpenCut->Fill(angle_open);
      if(angle_open>0.05)// + 0.5/egamma)
        hThetaCut->Fill(theta);

      nIn++;
    }else
      nOut++;
    fNtuple->Fill(data);
  }


  hNin->Fill(nIn);
  hNout->Fill(nOut);
 

}

void   HarpoAnalyseVertex3Dsim::Init()
{
  
  fChooseNbins = 0;

  hPhiPmMC = new TH1F("hPhiPmMC","",200,-TMath::Pi(),TMath::Pi());
  hPhiPm = new TH1F("hPhiPm","",200,-TMath::Pi(),TMath::Pi());
  hPhiPmCut = new TH1F("hPhiPmCut","",200,-TMath::Pi(),TMath::Pi());
  hDeltaPhiPm = new TH1F("hDeltaPhiPm","",200,-TMath::Pi(),TMath::Pi());
  hOpenMC = new TH1F("hOpenMC","",200,0,TMath::Pi());
  hOpen = new TH1F("hOpen","",200,0,TMath::Pi());
  hDeltaOpen = new TH1F("hDeltaOpen","",200,-TMath::Pi(),TMath::Pi());
  hOpenCut = new TH1F("hOpenCut","",200,0,TMath::Pi());
  hDistRecoSim = new TH3F("hDistRecoSim","",200,-50,50,200,-50,50,200,-50,50);
  hPosition = new TH3F("hPosition","",100,0,300,100,0,300,128,0,512);
  hPositionCut = new TH3F("hPositionCut","",100,0,300,100,0,300,128,0,512);

  hTheta = new TH1F("hTheta","",200,0,TMath::Pi());
  hThetaCut = new TH1F("hThetaCut","",200,0,TMath::Pi());

  hNin = new TH1F("hNin","",20,0,20);
  hNout = new TH1F("hNout","",20,0,20);

  fNtuple = new TNtupleD("fNtuple","Vertex data","x:y:z:px:py:pz:pBx:pBy:pBz:pxsim:pysim:pzsim:e:e1:e2:theta:phipm:open:cut");
  fNtuple->SetDirectory(0);


 }

void   HarpoAnalyseVertex3Dsim::Save(char * /* mode */)
 {

   OpenHistFile("vertex3dsim");
   
   // Save histograms here
   hDistRecoSim->Write();
   hPhiPm->Write();
   hPhiPmCut->Write();
   hPhiPmMC->Write();
   hDeltaPhiPm->Write();
   hOpen->Write();
   hDeltaOpen->Write();
   hOpenCut->Write();
   hOpenMC->Write();
   hPosition->Write();
   hPositionCut->Write();
   hNin->Write();
   hNout->Write();
   hTheta->Write();
   hThetaCut->Write();

   fNtuple->Write();
 }



void HarpoAnalyseVertex3Dsim::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* /* infobox */)
{
  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
  TCanvas* c = ecTab->GetCanvas();
  c->GetPad(1)->Delete();
  c->GetPad(2)->Delete();
  c->Divide(2);

  for(Int_t plane = 0; plane<2;plane++){
    TVirtualPad* cMap = c->GetPad(plane+1);
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    TH2F* h = new TH2F("h","",512,0,512,288,0,288);
    for(Int_t i=1;i<NADC;i++){ //Time bins
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = m->GetData(j,i);
        h->SetBinContent(i+1,j+1,q);
      }
    } 
    h->SetMinimum(0);
    MakeNiceHisto(h,cMap,"colz");
    h->Delete();
  }

  c->Update();
}


void HarpoAnalyseVertex3Dsim::ConfigFrame(TGMainFrame* fMain, Int_t id)
{
  // in hrecomonitorgui: Creates a popup window for analysis configuration
  // You must define SetConfig() properly

  UInt_t xsize = 200;
  UInt_t ysize2 = 20;
  UInt_t ysize = 10*ysize2;
  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
  main->DontCallClose(); // to avoid double deletions.
  
  // use hierarchical cleaning
  main->SetCleanup(kDeepCleanup);

  TGVerticalFrame* frame = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);

  // Object layout options
  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);

  //________ DO NOT MODIFY ABOVE _____________________________




  // Title of the analysis
  TGLabel* fAnalysisLabel = new TGLabel(frame, "Template analysis");

  // Create a line for choosing value of parameter
  TGHorizontalFrame* nBinsFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* nBinsLabel = new TGLabel(nBinsFrame, "fNbins");
  fChooseNbins    = new TGNumberEntry(nBinsFrame);
  fChooseNbins->SetNumber(fNbins);

  main->AddFrame(frame,fLayout4);
  frame->AddFrame(fAnalysisLabel,fLayout3);
  frame->AddFrame(nBinsFrame,fLayout3);
  nBinsFrame->AddFrame(nBinsLabel,fLayout1);
  nBinsFrame->AddFrame(fChooseNbins,fLayout2);



  //________ DO NOT MODIFY BELOW _____________________________
  // Button to validate configuration
  TGTextButton* setConf = new TGTextButton(frame,"Save Config",id);
  setConf->Associate(fMain);

  frame->AddFrame(setConf,fLayout3);

  main->MapSubwindows();
  main->MapWindow();
  main->Resize();
  return;
}



void   HarpoAnalyseVertex3Dsim::SetConfig()
{
  // Update the configuration according to the values in the popup window

  if(!fChooseNbins) return;

  fNbins = fChooseNbins->GetNumber();
}