HarpoAnalyseVertex3D.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseVertex3D.cxx
//
#include "HarpoAnalyseVertex3D.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 "TSystem.h"

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

ClassImp(HarpoAnalyseVertex3D);

void   HarpoAnalyseVertex3D::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   HarpoAnalyseVertex3D::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  //???????????????????????????????
  if(!fEvt->GetRecoEvent()) return;


  Double_t angle_open_sim = 0;
  Double_t angle_omega_sim = 0;
  // Double_t psim[2][3];
  TVector3 possim = TVector3(0,0,0);
  TVector3 pgammasim = TVector3(0,0,1);
  TVector3 ppolasim = TVector3(0,1,0);
  TVector3 pgamma = TVector3(0,0,1);
  TVector3 ppola = TVector3(0,1,0);
  TVector3 ppositron, pelectron;
  Double_t epositron = -1.0, eelectron = -1.0;
  if (gHDetSet->isExist(SIMDET)){
    HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
    Int_t nIonTr = simEvt->GetNtracks();
    for(Int_t i = 0; i<nIonTr; i++){
      TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
      Int_t nPoints = tr->GetNpoints();
      for(Int_t j = 0; j<nPoints; j++){
        //      TpcSimIonisationPoint* point = tr->GetPoint(j);
        //hNeSim->Fill(point->GetNe());
      }
    }
    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(); 
        ppola = track->GetPola(); 
        break;
      case 11: 
        palready = kTRUE;
        ppositron = track->GetP(); 
        epositron = track->GetE(); 
        break;
      case -11: 
        pelectron = track->GetP(); 
        eelectron = track->GetE(); 
        break;
      // default:
      //        psim[i][0] = track->GetPx();
      //        psim[i][1] = track->GetPy();
      //        psim[i][2] = track->GetPz();
      //        norm[i] = TMath::Sqrt(psim[i][0]*psim[i][0]+psim[i][1]*psim[i][1]+psim[i][2]*psim[i][2]);
        
      }
      //hPsim->Fill(track->GetPx(),track->GetPy(),track->GetPz());
    }
    Double_t u1u2 = ppositron.Unit()*pelectron.Unit();
    
    angle_open_sim=acos(u1u2);

    TVector3 ppola2 = ppolasim.Cross(pgammasim.Unit());
    TVector3 pPlus = TVector3(ppositron.Dot(pgammasim.Unit()),ppositron.Dot(ppolasim.Unit()),ppositron.Dot(ppola2.Unit())).Unit();
    TVector3 pMinus = TVector3(pelectron.Dot(pgammasim.Unit()),pelectron.Dot(ppolasim.Unit()),pelectron.Dot(ppola2.Unit())).Unit();
    angle_omega_sim = atan2((pPlus.Px()*pMinus.Py() - pMinus.Px()*pPlus.Py()),
                            (pPlus.Px()*pMinus.Pz() - pMinus.Px()*pPlus.Pz()));

    hOpenSim->Fill(angle_open_sim);
    hOmegaSim->Fill(angle_omega_sim);
    //}
  }





  HarpoRecoEvent* reco = fEvt->GetRecoEvent();

  const Int_t kNdata = 41;
  Double_t data[kNdata];
  for(Int_t i = 0; i<kNdata; i++) data[i] = -1000;
  for(Int_t i = 0; i<3; i++){
    data[10+i] = ppositron[i];
    data[13+i] = pelectron[i];
  }
  data[32] = epositron;
  data[33] = eelectron;

  TVector3 psumsim = ppositron + pelectron;
  data[23] = psumsim.Unit().x();
  data[24] = psumsim.Unit().y();
  data[25] = psumsim.Unit().z();

  Double_t xStart[2] = {-1,-1};
  Double_t zStart[2] = {1000,1000};
  TClonesArray* clArray = reco->GetClustersArray();
  Int_t nCl = clArray->GetEntries();
  for(Int_t iCl = 0; iCl<nCl; iCl++){
    HarpoRecoClusters* cl = (HarpoRecoClusters*)clArray->At(iCl);
    Int_t plane = cl->GetPlane();
    if(cl->GetZ()<zStart[plane]){
      xStart[plane] = cl->GetX();
      zStart[plane] = cl->GetZ();
    }
  }

  if(gHDetSet->isExist(PMM2)){
    Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
    data[34] = anaEvt->GetTriggerType();
  }

  Double_t fTstart = 90;
  Double_t fTend = 430;
  Double_t fThr = 1e4;
  for(Int_t plane = 0; plane<2; plane++){
    if (!gHDetSet->isExist(plane)) continue;
    // Process RAW data
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    if ( m == NULL ) continue;
    //    Int_t i2 = 0;
    Double_t qTotIn = 0;
    Double_t qTotBefore = 0;
    Double_t qTotAfter = 0;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      for(Int_t i=0;i<NADC;i++){ //Time bins
        Double_t q = m->GetData(j,i);
        if(q<0) continue;
        if(i< fTstart            ) qTotBefore += q;
        if(              i >fTend) qTotAfter += q;
        if(i>=fTstart && i<=fTend) qTotIn += q;
      }
    }
    if(qTotBefore>fThr || qTotAfter>fThr)
      data[34] = -2;
    data[35] = qTotBefore;
    data[36] = qTotAfter;
  }






  data[0] = xStart[XPLANE];
  data[1] = xStart[YPLANE];
  data[2] = 0.5*(zStart[XPLANE]+zStart[XPLANE]);
  data[26] = xStart[XPLANE];
  data[27] = xStart[YPLANE];
  data[28] = 0.5*(zStart[XPLANE]+zStart[XPLANE]);
  data[29] = possim.x();
  data[30] = possim.y();
  data[31] = possim.z();
  data[5] = angle_omega_sim;
  data[6] = angle_open_sim;
  data[22] = -1;
  fNtuple->Fill(data);
  

  TClonesArray* vArray = reco->GetVertex3DArray();
  Int_t nV = vArray->GetEntries();
  for(Int_t iV = 0; iV<nV; iV++){
    HarpoRecoVertex3D* vertex = (HarpoRecoVertex3D*)vArray->At(iV);

    TVector3 p1 = vertex->GetP1Vertex3D();
    TVector3 p2 = vertex->GetP2Vertex3D();

    Double_t u1u2 = p1.Unit()*p2.Unit();
    Double_t angle_open=acos(u1u2);

    TVector3 ppola2 = ppola.Cross(pgamma);
    TVector3 pPlus = TVector3(p1.Dot(pgamma.Unit()),p1.Dot(ppola.Unit()),p1.Dot(ppola2.Unit()));
    TVector3 pMinus = TVector3(p2.Dot(pgamma.Unit()),p2.Dot(ppola.Unit()),p2.Dot(ppola2.Unit()));
    Double_t angle_omega = atan2((pPlus.Px()*pMinus.Py() - pMinus.Px()*pPlus.Py()),
                                 (pPlus.Px()*pMinus.Pz() - pMinus.Px()*pPlus.Pz()));
    
    data[0] = vertex->GetXvertex3D();
    data[1] = vertex->GetYvertex3D();
    data[2] = vertex->GetZvertex3D();
    data[3] = angle_omega;
    data[4] = angle_open;
    TVector3 psum = (p1 + p2).Unit();
    data[7] = psum.x();
    data[8] = psum.y();
    data[9] = psum.z();

    Double_t d1[2] = {0,0};
    Double_t d2[2] = {0,0};
    d1[0] = (ppositron - p1).Mag();
    d2[0] = (ppositron - p2).Mag();
    d1[1] = (pelectron - p1).Mag();
    d2[1] = (pelectron - p2).Mag();

    Int_t j = 0;
    if(d1[0]+d2[1]<d1[1]+d2[0]){
      j = 1;
      data[32] = eelectron;
      data[33] = epositron;
    }
    for(Int_t i = 0; i<3; i++){
      if(j){
        data[10+i] = pelectron[i];
        data[13+i] = ppositron[i];
      }else{
        data[10+i] = ppositron[i];
        data[13+i] = pelectron[i];
      }
      data[16+i] = p1[i];
      data[19+i] = p2[i];
    }
    data[22] = iV;






    Double_t fWidthX = 50, fWidthZ = 50;
    Double_t q1[2] = {0,0}, q2[2] = {0,0};
    for(Int_t plane = 0; plane<2; plane++){
      HarpoDccMap* m = fEvt->GetDccMap(plane);
      for(Int_t i=0;i<NADC;i++){ //Time bins
        for(Int_t j=0;j<NCHAN;j++){ // Channels
          Double_t q = m->GetData(j,i);
          if(q<=0) continue;
          if(TMath::Abs(j-data[plane]-p1[plane]/p1[2]*(i-data[2]))>fWidthX  && TMath::Abs(j-data[plane]-p2[plane]/p2[2]*(i-data[2]))>fWidthX ) continue;
          if(TMath::Abs(i-data[2]-p1[2]/p1[plane]*(j-data[plane]))>fWidthZ && TMath::Abs(i-data[2]-p2[2]/p2[plane]*(j-data[plane]))>fWidthZ) continue;
          Double_t r = 1; // r is the fraction of area in the sector containing the corner i,j
          if((i-data[2])*psum[plane]-(j-data[plane])*psum[2]>0) r = 1. - r;
          q1[plane] += q*r;
          q2[plane] += q*(1.-r);
        }
      }
    }
    

    data[37] = q1[0];
    data[38] = q1[1];
    data[39] = q2[0];
    data[40] = q2[1];








    fNtuple->Fill(data);

  }
 

}

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

  hDist = new TH1F("hDist","",500,0,500);
  hP = new TH2F("hP","",100,-1,1,100,-1,1);

  hOmega = new TH1F("hOmega","",200,-3.5,3.5);
  hOpen = new TH1F("hOpen","",200,0,3.5);
  hDeltaOmega = new TH1F("hDeltaOmega","",200,-3.5,3.5);
  hDeltaOpen = new TH1F("hDeltaOpen","",200,-3.5,3.5);
  hOmegaSim = new TH1F("hOmegaSim","",200,-3.5,3.5);
  hOpenSim = new TH1F("hOpenSim","",200,0,3.5);
  hDiffOpen = new TH1F("hDiffOpen","",200,-0.5,0.5);


  gBeam = new TGraph2D();

  fNtuple = new TNtupleD("fNtuple","Vertex data","x0:y0:z0:omega:alpha:omegasim:alphasim:px:py:pz:px1sim:py1sim:pz1sim:px2sim:py2sim:pz2sim:px1:py1:pz1:px2:py2:pz2:iV:pxsim:pysim:pzsim:xS:yS:zS:xsim:ysim:zsim:e1:e2:trig:qBefore:qAfter:qx1:qy1:qx2:qy2");
  fNtuple->SetDirectory(0);


 }

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

   OpenHistFile("vertex3d");

   // Save histograms here
   hDist->Write();
   gBeam->Write("gBeam");
   hP->Write();
   hOmega->Write();
   hOpen->Write();
   hDeltaOmega->Write();
   hDeltaOpen->Write();
   hOmegaSim->Write();
   hOpenSim->Write();
   hDiffOpen->Write();

   fNtuple->Write();
 }



void HarpoAnalyseVertex3D::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 HarpoAnalyseVertex3D::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   HarpoAnalyseVertex3D::SetConfig()
{
  // Update the configuration according to the values in the popup window

  if(!fChooseNbins) return;

  fNbins = fChooseNbins->GetNumber();
}