run_ana_invariantmass_4pi_stt.C

Go to the documentation of this file.

class TCandList;
class TCandidate;
class TFitParams;

int run_ana_invariantmass_4pi_stt(int nEntries=0)
{
  gStyle->SetOptFit(1011);

  TStopwatch timer;
  timer.Start();

  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");

  TString inPidFile  = "evt_pid_stt.root";
  TString inSimFile = "evt_points_stt.root";
  TFile *inFile = TFile::Open(inSimFile,"READ");
  TTree *tree=(TTree *) inFile->Get("pndsim") ;
  tree->AddFriend("pndsim",inPidFile);
  TClonesArray* mc_array=new TClonesArray("PndMCTrack");
  tree->SetBranchAddress("MCTrack",&mc_array);

  TClonesArray* cand_array=new TClonesArray("PndPidCandidate");
  tree->SetBranchAddress("PidChargedCand",&cand_array);

    FairMCEventHeader* evthead;
   tree->SetBranchAddress("MCEventHeader.", &evthead);

  TFile *out=TFile::Open("invariantmass_4pi_stt.root","RECREATE");

  PndEventReader evr(inPidFile);

  TH1F *nc=new TH1F("nc","Number of Charged Tracks; Charged Tracks",20,-0.5,19.5);
  TH1F *invmassnosel=new TH1F("invmassnosel","2(#pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,0,10);
  TH1F *invmassnocut=new TH1F("invmassnocut","2(#pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,2,4);
  // TH1F *invmassvtx=new TH1F("invmassvtx","2(#pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,2,4);
  TH1F *invmasswithpid=new TH1F("invmasswithpid","#pi^{+}#pi^{-} Invariant mass;Invariant Mass (GeV)",100,0,10);
  TH1F *invmasswithpid_sel=new TH1F("invmasswithpid_sel","#pi^{+}#pi^{-} Invariant mass;Invariant Mass (GeV)",100,2,4);
TH1F *invmass_trackhighmom= new TH1F("invmass_trackhighmom","2(#pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,2,4);
 TH1F *invmasschicut=new TH1F("invmasschicut","#2(pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,2,4);
 TH1F *invmasschicut_best=new TH1F("invmasschicut_best","2(#pi^{+}#pi^{-}) Invariant mass;Invariant Mass (GeV)",100,2,4);

 TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-5,5,100,-5,5);
 TH1F *hvzpos = new TH1F("hvzpos","z position of fitted decay vertex",100,-4,4);

 TH1F *chivtx=new TH1F("chivtx","Chi Square RhoKinVtxFitter; Chi Square / N_{df}",100,0,100);
  TH1F *hvtxresX = new TH1F("hvtxresX","X resolution of fitted decay vertex",100,-0.3,0.3);
  TH1F *hvtxresY = new TH1F("hvtxresY","Y resolution of fitted decay vertex",100,-0.3,0.3);
  TH1F *hvtxresZ = new TH1F("hvtxresZ","Z resolution of fitted decay vertex",100,-0.3,0.3);


   TPidPlusSelector *piplusSel=new TPidPlusSelector("piplus");
  TPidMinusSelector *piminusSel=new TPidMinusSelector("piminus");

  TCandList pp, pm, pipi,pipinosel, pipiwithpid;

  int n_reco=0;

  TLorentzVector ini(0,0,4.0,5.04684);
  int i=0,j=0,k=0,l=0,y=0;

  TLorentzVector beam;
  beam.SetXYZM (0.0,0.0,4.0,0.938272);
  TLorentzVector target;
  target.SetXYZM (0.0,0.0,0.0,0.938272);
  TLorentzVector pbarp=beam+target;



  float energypip, energypim, momxpip, momypip, momzpip;
  TPidMassSelector *pipisel=new TPidMassSelector("pipisel",3.07,0.5); //energy in the center of mass 3.07GeV
  TPidPlusSelector *piplusSel=new TPidPlusSelector("piplus");
  TPidMinusSelector *piminusSel=new TPidMinusSelector("piminus");

    if (nEntries==0) nEntries=evr.GetEntries();
    while (evr.GetEvent() && i++<nEntries){

      if (!((i+1)%100)) cout<<"evt " << i << "\n";

      evr.FillList(pp,"Charged");
    evr.FillList(pm,"Charged");

    pp.Select(piplusSel);
    pm.Select(piminusSel);

 int nchrg=pp.GetLength();
    nc->Fill(nchrg);


    float momentumpplus, momentumpminus, theta1, theta2, phi1, phi2;

   int n_removed=0;
    int ii=0;
    for (Int_t l=0;l<pp.GetLength();l++){
      ii=l-n_removed;
      if((pp[ii].GetMicroCandidate()->GetSttHits())==0){
        pp.Remove(pp[ii]);
        n_removed++;
      }
    }

    int n_removed=0;
    int ii=0;
    for (Int_t l=0;l<pm.GetLength();l++){
      ii=l-n_removed;
      if((pm[ii].GetMicroCandidate()->GetSttHits())==0){
        pm.Remove(pm[ii]);
        n_removed++;
      }
    }



   for (Int_t l=0;l<pp.GetLength();l++){
      pp[l].SetMass(TRho::Instance()->GetPDG()->GetParticle(211)->Mass());
      momentumpplus=pp[l].GetMicroCandidate().GetMomentum().Mag();

    }
    for (Int_t l=0;l<pm.GetLength();l++){
      pm[l].SetMass(TRho::Instance()->GetPDG()->GetParticle(211)->Mass());
      momentumpminus=pm[l].GetMicroCandidate().GetMomentum().Mag();

}


    pipinosel.Combine(pp,pm,pp,pm);

    for (y=0;y<pipinosel.GetLength();++y){
      invmassnosel->Fill(pipinosel[y].M());
    }

    pipi.Combine(pp,pm,pp,pm);
    pipi.Select(pipisel);

    for (y=0;y<pipi.GetLength();++y){
      invmassnocut->Fill(pipi[y].M());
            if (momentumpplus>0.3 && momentumpminus >0.3){
        invmass_trackhighmom->Fill(pipi[y].M());
        }

}


  //MonteCarlo PID
    tree->GetEntry(i-1);

    TVector3 mcVertex;
    evthead->GetVertex(mcVertex);

// MC PID
    // Leave only pions in particle lists
    int n_removed=0;
    int ii=0;
    for (l=0;l<pp.GetLength();++l) {
      ii=l-n_removed;
      if (pp[ii].GetMicroCandidate().GetMcIndex()>-1){
        PndMCTrack *mcTrack = (PndMCTrack*)mc_array->At(pp[ii].GetMicroCandidate().GetMcIndex());
        if (mcTrack!=0)
          {
            if ((mcTrack->GetPdgCode()!=211) || (mcTrack->GetMotherID()!=-1))
              {
                pp.Remove(pp[ii]);
                n_removed++;
              }
          }
        else
          {
            std::cout<<"Pi plus, element "<<l<<" has no assosiated mcTRack"<<std::endl;
          }
      }
    }


    n_removed=0;
    ii=0;
    for (l=0;l<pm.GetLength();++l) {
      ii=l-n_removed;
      if (pm[ii].GetMicroCandidate().GetMcIndex()>-1){
        PndMCTrack *mcTrack = (PndMCTrack*)mc_array->At(pm[ii].GetMicroCandidate().GetMcIndex());
        if (mcTrack!=0)
          {
            if ((mcTrack->GetPdgCode()!=-211) || (mcTrack->GetMotherID()!=-1))
              {
                pm.Remove(pm[ii]);
                n_removed++;
              }
          }
        else
          {
            std::cout<<"Pi Minus,  element "<<l<<" has no assosiated mcTRack"<<std::endl;
          }
      }
    }


    pipiwithpid.Combine(pp,pm,pp,pm);
    for (y=0;y<pipiwithpid.GetLength();++y){
      invmasswithpid->Fill(pipiwithpid[y].M());
    }


    pipiwithpid.Combine(pp,pm,pp,pm);
    pipiwithpid.Select(pipisel); //Mass selector
    for (y=0;y<pipiwithpid.GetLength();++y){
      invmasswithpid_sel->Fill(pipiwithpid[y].M());
    }


      int best_i=0;
    double best_chi2=1000;
    TCandidate *pipifit_best=0;
    Float_t pipivtx_mass;
  for (y=0;y<pipiwithpid.GetLength();++y){

      RhoKinVtxFitter vtxfitter(pipiwithpid[y]);
      vtxfitter.Fit();
      TCandidate *pipifit=vtxfitter.FittedCand(pipiwithpid[y]);

  TVector3 pipiVtx=pipifit->Pos();
      double chi2_vtx=vtxfitter.GlobalChi2();

      hvpos->Fill(pipiVtx.X(),pipiVtx.Y());
      hvzpos->Fill(pipiVtx.Z());

      if(chi2_vtx<best_chi2)
        {
          best_chi2=chi2_vtx;
          best_i=y;
          pipifit_best=pipifit;
          pipivtx_mass=pipifit_best->M();
          bestPos = pipifit->Pos();
        }
      chivtx->Fill(chi2_vtx/5); // Number degree of freedom 2N-3=5; N=number of charged tracks

  }

  if((pipiwithpid.GetLength()!=0))
    {
      invmasschicut_best->Fill(pipivtx_mass);
      n_reco++;
      hvtxresX->Fill(mcVertex.X()-bestPos.X());
      hvtxresY->Fill(mcVertex.Y()-bestPos.Y());
      hvtxresZ->Fill(mcVertex.Z()-bestPos.Z());
    }

    }

  out->cd();
  out->Write();
  out->Save();


  return 0;
}