ana_check_psi.C

Go to the documentation of this file.

class TCandList;
class TCandidate;
class TFitParams;

int ana_check_psi(int nevts=0)
{
  TString OutFile = "output.root";

  gStyle->SetOptFit(1011);

  TStopwatch timer;
  timer.Start();

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

  TString inPidFile  = "evt_pid.root";
  TString inSimFile = "evt_points.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);

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

  TFile *out = TFile::Open(OutFile,"RECREATE");

  // the PndEventReader takes care about file/event handling
  PndEventReader evr(inPidFile);

  TH1F *h_psi_nocut=new TH1F("h_psi_nocut","m(psi), (no cuts);E, GeV",100,3.2,4.4);
  TH1F *h_mD1_nocuts=new TH1F("h_mD1_nocuts","D+: m(#pi+ #pi+ K-) (no cuts)",100,1.5,2.2);
  TH1F *h_mD2_nocuts=new TH1F("h_mD2_nocuts","D-: m(#pi- #pi- K+) (no cuts)",100,1.5,2.2);
  TH1F *h_mD1_pid=new TH1F("h_mD1_pid","D+: m(#pi+ #pi+ K-) (pid)",100,1.5,2.2);
  TH1F *h_mD2_pid=new TH1F("h_mD2_pid","D-: m(#pi- #pi- K+) (pid)",100,1.5,2.2);

  TH1F *h_psi_pid=new TH1F("h_psi_pid","m(#psi), (MC PID);E, GeV",100,3.2,4.4);
  TH1F *h_psi_pid_vtx=new TH1F("h_psi_pid_vtx","m(#psi), (MC PID);E, GeV",100,3.2,4.4);

  TH1F *h_psi_4c=new TH1F("h_psi_4c","m(#psi), 4C-fit",100,3.2,4.4);
  TH1F *h_mD_4c=new TH1F("h_mD_4c","D: m(K #pi #pi) (4C-fit)",100,1.5,2.2);

  TH1F *h_psi_vtx=new TH1F("h_psi_vtx","m(#psi), Vertex fit",100,3.2,4.4);
  TH1F *h_D1_vtx_bef=new TH1F("h_D1_vtx_bef","#phi: m(K #pi #pi) (Vertex fit)",100,1.5,2.2);
  TH1F *h_D2_vtx_bef=new TH1F("h_D2_vtx_bef","#phi: m(K #pi #pi) (Vertex fit)",100,1.5,2.2);
  TH1F *h_D1_vtx=new TH1F("h_D1_vtx","#phi: m(K #pi #pi) (Vertex fit)",100,1.5,2.2);
  TH1F *h_D2_vtx=new TH1F("h_D2_vtx","#phi: m(K #pi #pi) (Vertex fit)",100,1.5,2.2);

  TH1F *h_psi_Dmass=new TH1F("h_psi_Dmass","m(#psi), (cut on D mass);E, GeV",100,3.2,4.4);
  TH1F *h_mD_final=new TH1F("h_mD_final","D: m(K #pi #pi)",100,1.5,2.2);

  TH1F *nc=new TH1F("nc","n charged",20,0,20);

  TH1F *h_chi2_4c=new TH1F("h_chi2_4c","#chi^{2} 4C-fit;#chi^{2}/N_{df}",100,0,100);
  TH1F *h_chi2b_4c=new TH1F("h_chi2b_4c","#chi^{2} 4C-fit;#chi^{2}/N_{df}",100,0,100);
  TH1F *h_chi2_vtx_D1=new TH1F("h_chi2_vtx_D1","#chi^{2} vertex;#chi^{2}/N_{df}",100,0,100);
  TH1F *h_chi2_vtx_D2=new TH1F("h_chi2_vtx_D2","#chi^{2} vertex;#chi^{2}/N_{df}",100,0,100);
  TH2F *hvpos_D1 = new TH2F("hvpos_D1","(x,y) projection of fitted decay vertex",100,-5,5,100,-5,5);
  TH1F *hvzpos_D1 = new TH1F("hvzpos_D1","z position of fitted decay vertex",100,-10,10);
  TH1F *hvtxresX_D1 = new TH1F("hvtxresX_D1","X resolution of fitted decay vertex",100,-0.1,0.1);
  TH1F *hvtxresY_D1 = new TH1F("hvtxresY_D1","Y resolution of fitted decay vertex",100,-0.1,0.1);
  TH1F *hvtxresZ_D1 = new TH1F("hvtxresZ_D1","Z resolution of fitted decay vertex",100,-0.1,0.1);
  TH2F *hvpos_D2 = new TH2F("hvpos_D2","(x,y) projection of fitted decay vertex",100,-5,5,100,-5,5);
  TH1F *hvzpos_D2 = new TH1F("hvzpos_D2","z position of fitted decay vertex",100,-10,10);
  TH1F *hvtxresX_D2 = new TH1F("hvtxresX_D2","X resolution of fitted decay vertex",100,-0.1,0.1);
  TH1F *hvtxresY_D2 = new TH1F("hvtxresY_D2","Y resolution of fitted decay vertex",100,-0.1,0.1);
  TH1F *hvtxresZ_D2 = new TH1F("hvtxresZ_D2","Z resolution of fitted decay vertex",100,-0.1,0.1);

  TPidMassSelector *psiMassSel=new TPidMassSelector("psi",1.87,0.07);

  TPidPlusSelector *pplusSel=new TPidPlusSelector("pplus");
  TPidMinusSelector *pminusSel=new TPidMinusSelector("pminus");

  // the candidates lists we need
  TCandList p1, p2, k1, k2, D1, D1_pid, D2, D2_pid, psi, psi_pid,  psi_nocut;

  int n_reco=0;
  // Number of events in file and number of reconstructed eta_c to store in root file
  TH1F *n_events=new TH1F("n_events","total number of events",1,0,1);
  TH1F *n_psi=new TH1F("n_psi","number of reconstructed eta_c",1,0,1);

  TLorentzVector ini(0,0,6.578800,7.583333);

  if (nevts==0) nevts=evr.GetEntries();

  n_events->SetBinContent(1,nevts);
  // cout << "nevts " << nevts << "\n";
  int i=0,j=0, k=0, l=0;

  // *************
  // this is the loop through the events ... as simple as this...
  // ****************
  while (evr.GetEvent() && i++<nevts)
    {
      //.cout << "evt: " << i << endl;
      if ((i%100)==0)  cout<<"evt " << i << "\n";
      //if (!((i+1)%100)) cout<<"evt " << i << "\n";
      evr.FillList(p1,"Charged");
      evr.FillList(p2,"Charged");
      evr.FillList(k1,"Charged");
      evr.FillList(k2,"Charged");

      p1.Select(pplusSel);
      p2.Select(pminusSel);
      k1.Select(pplusSel);
      k2.Select(pminusSel);

      int nchrg=p1.GetLength()+p2.GetLength();
      nc->Fill(nchrg);

      for (j=0;j<p1.GetLength();++j) {
        p1[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(211)->Mass());
      }
      for (j=0;j<p2.GetLength();++j) {
        p2[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(211)->Mass());
      }
      for (j=0;j<k1.GetLength();++j) {
        k1[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(321)->Mass());
      }
      for (j=0;j<k2.GetLength();++j) {
        k2[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(321)->Mass());
      }


      D1.Combine(p1,p1, k2);
      D2.Combine(p2,p2, k1);

      for (j=0;j<D1.GetLength();++j) h_mD1_nocuts->Fill(D1[j].M());
      for (j=0;j<D2.GetLength();++j) h_mD2_nocuts->Fill(D2[j].M());

      D1.Select(psiMassSel);
      D2.Select(psiMassSel);
      psi_nocut.Combine(D1,D2);

      for (l=0;l<psi_nocut.GetLength();++l) {
        h_psi_nocut->Fill(psi_nocut[l].M());
      }

      tree->GetEntry(i-1);
      TVector3 mcVertex, mcD1Vertex, mcD2vertex;

      if (((PndMCTrack*)mc_array->At(0))!=0) mcD1Vertex = ((PndMCTrack*)mc_array->At(0))->GetStartVertex();
      else
        cout << "Not found k1!" << endl;
      if (((PndMCTrack*)mc_array->At(3))!=0) mcD2Vertex = ((PndMCTrack*)mc_array->At(3))->GetStartVertex();
      else
        cout << "Not found k2!" << endl;

      evthead->GetVertex(mcVertex);

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

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

      n_removed=0;
      ii=0;
      for (l=0;l<k1.GetLength();++l) {
        ii=l-n_removed;
        if (k1[ii].GetMicroCandidate().GetMcIndex()>-1){
          PndMCTrack *mcTrack = (PndMCTrack*)mc_array->At(k1[ii].GetMicroCandidate().GetMcIndex());
          if (mcTrack!=0)
            {
              if ((mcTrack->GetPdgCode()!=321) || (mcTrack->GetMotherID()!=-1))
                {
                  k1.Remove(k1[ii]);
                  n_removed++;
                }
            }
          else
            {
              std::cout<<"stt h: " << k1[ii].GetMicroCandidate().GetSttHits() << std::endl;
              std::cout<<"Kaon list 2, element "<<l<<" has no assosiated mcTRack"<<std::endl;
            }
        }
      }

      n_removed=0;
      ii=0;
      for (l=0;l<k2.GetLength();++l) {
        ii=l-n_removed;
        if (k2[ii].GetMicroCandidate().GetMcIndex()>-1){
          PndMCTrack *mcTrack = (PndMCTrack*)mc_array->At(k2[ii].GetMicroCandidate().GetMcIndex());
          if (mcTrack!=0)
            {
              if ((mcTrack->GetPdgCode()!=-321) || (mcTrack->GetMotherID()!=-1))
                {
                  k2.Remove(k2[ii]);
                  n_removed++;
                }
            }
          else
            {
              std::cout<<"stt h: " << k2[ii].GetMicroCandidate().GetSttHits() << std::endl;
              std::cout<<"Kaon list 2, element "<<l<<" has no assosiated mcTRack"<<std::endl;
            }
        }
      }

      D1_pid.Combine(p1,p1, k2);
      D2_pid.Combine(p2,p2, k1);

      for (j=0;j<D1_pid.GetLength();++j) h_mD1_pid->Fill(D1_pid[j].M());
      for (j=0;j<D2_pid.GetLength();++j) h_mD2_pid->Fill(D2_pid[j].M());
      D1_pid.Select(psiMassSel);
      D2_pid.Select(psiMassSel);
      psi.Combine(D1_pid,D2_pid);

      for (l=0;l<psi.GetLength();++l) {
        h_psi_pid->Fill(psi[l].M());
      }

      cout << " ";// << endl;
      if (use4cfit)
        {
          int best_i=0;
          double best_chi2=100000;
          TCandidate *ccfit;// = new TCandidate();
          double m_D1, m_D2;
          for (l=0;l<psi.GetLength();++l) {
            Rho4CFitter fitter(psi[l],ini);
            fitter.FitConserveMasses();
            double chi2=fitter.GetChi2();
            if (chi2<best_chi2)
              {
                best_chi2=chi2;
                best_i = l;
                ccfit = (TCandidate*)fitter.FittedCand(psi[l]);
                TCandidate *D1best = fitter.FittedCand(*(psi[l].Daughter(0)));
                TCandidate *D2best = fitter.FittedCand(*(psi[l].Daughter(1)));
                TCandidate *p1best = fitter.FittedCand(*(psi[l].Daughter(0)->Daughter(0)));
                TCandidate *p2best = fitter.FittedCand(*(psi[l].Daughter(0)->Daughter(1)));
                TCandidate *k1best = fitter.FittedCand(*(psi[l].Daughter(0)->Daughter(2)));
                TCandidate *p3best = fitter.FittedCand(*(psi[l].Daughter(1)->Daughter(0)));
                TCandidate *p4best = fitter.FittedCand(*(psi[l].Daughter(1)->Daughter(1)));
                TCandidate *k2best = fitter.FittedCand(*(psi[l].Daughter(1)->Daughter(2)));
                TLorentzVector tlvk1 = p1best->P4();
                TLorentzVector tlvk2 = p2best->P4();
                TLorentzVector tlvk3 = k1best->P4();
                TLorentzVector tlvk4 = p3best->P4();
                TLorentzVector tlvk5 = p4best->P4();
                TLorentzVector tlvk6 = k2best->P4();
                m_D1= (tlvk1+tlvk2+tlvk3).M();
                m_D2= (tlvk4+tlvk5+tlvk6).M();
              }
            h_chi2_4c->Fill(chi2/9); // Ndf=3N-3=9
          }

          if (psi.GetLength()!=0) cout << "evt: " << i << "\tbest: " << best_chi2 << "\tlen " << psi.GetLength() << endl;

          if(/*(best_chi2<270)&&*/(psi.GetLength()!=0))
            {
              h_chi2b_4c->Fill(best_chi2/9); // Ndf=3N-3=9
              h_psi_4c->Fill(ccfit->M());
              h_mD_4c->Fill(m_D1);
              h_mD_4c->Fill(m_D2);
              h_mD_final->Fill(m_D1);
              h_mD_final->Fill(m_D2);
              if (((m_D1>1.87-0.07)&&(m_D1<1.87+0.07))&&((m_D2>1.87-0.07)&&(m_D2<1.87+0.07)))
                {
                  h_psi_Dmass->Fill(psi[best_i].M());
                  if ((psi[best_i].M()>3.7)&&(psi[best_i].M()<3.9))
                    n_reco++;
                }
            }
        }
      else // use vertex fit
        {
          TCandList psi_vtx, D1_vtx, D2_vtx;
          TCandidate *ck1, *ck2, *cp1, *cp2, *cp3, *cp4, *D1tmp, *D2tmp, *D1_tmp, *D2_tmp, *psi_tmp;

          //Combine 4 kaons directly to candidates
          for (j=0;j<psi.GetLength();++j)
            {
              D1tmp=psi[j].Daughter(0);
              D2tmp=psi[j].Daughter(1);

              cp1=D1tmp->Daughter(0);
              cp2=D1tmp->Daughter(1);
              ck1=D1tmp->Daughter(2);

              cp3=D2tmp->Daughter(0);
              cp4=D2tmp->Daughter(1);
              ck2=D2tmp->Daughter(2);
              D1_tmp=cp1->Combine(*cp2,*ck1);
              D2_tmp=cp3->Combine(*cp4,*ck2);
              //psi_tmp=cp1->Combine(*cp2,*ck1,*cp3, *cp4, *ck2);
              D1_vtx.Add(*D1_tmp);
              D2_vtx.Add(*D2_tmp);
              //psi_vtx.Add(*psi_tmp);
            }

          for (l=0;l<psi.GetLength();++l) {
            // h_psi_pid_vtx->Fill(psi_vtx[l].M());
          }

          {
            int best_i=0;
            double best_chi2=1000;
            TCandidate *D1fit_best=0;
            TCandidate *k1fit_best, *p1fit_best, *p2fit_best, *k1nofit, *p1nofit, *p2nofit;
            TVector3 bestPos;
            Float_t D1vtx_mass, pre_m_D1;
            for (j=0;j<D1_vtx.GetLength();++j)
              {
                RhoKinVtxFitter vtxfitter(D1_vtx[j]);        // instantiate a vertex fitter
                vtxfitter.Fit();                          // do the vertex fit

                TCandidate *D1fit=vtxfitter.FittedCand(D1_vtx[j]);  // request the fitted Psi candidate
                TVector3 D1Vtx=D1fit->Pos();                    // and the decay vertex position
                double chi2_vtx_D1=vtxfitter.GlobalChi2();
                h_chi2_vtx_D1->Fill(chi2_vtx_D1/5); // Number degree of freedom 2N-3=5
                // plot mass and vtx x,y projection after fit
                hvpos_D1->Fill(D1Vtx.X(),D1Vtx.Y());
                hvzpos_D1->Fill(D1Vtx.Z());
                if(chi2_vtx_D1<best_chi2)
                {
                  best_chi2=chi2_vtx_D1;
                  best_i=l;
                  D1fit_best=D1fit;
                  p1fit_best=vtxfitter.FittedCand(*(D1fit_best->Daughter(0)));
                  p2fit_best=vtxfitter.FittedCand(*(D1fit_best->Daughter(1)));
                  k1fit_best=vtxfitter.FittedCand(*(D1fit_best->Daughter(2)));
                  p1nofit=D1_vtx[j].Daughter(0);
                  p2nofit=D1_vtx[j].Daughter(1);
                  k1nofit=D1_vtx[j].Daughter(2);
                  TLorentzVector tlvp1 = p1nofit->P4();
                  TLorentzVector tlvp2 = p2nofit->P4();
                  TLorentzVector tlvk1 = k1nofit->P4();
                  pre_m_D1= (tlvp1+tlvp2+tlvk1).M();

                  D1vtx_mass = D1fit_best->M();
                  bestPos = D1fit->Pos();
                }
              }
            if((psi.GetLength()!=0))
              {
                h_D1_vtx->Fill(D1vtx_mass);
                h_D1_vtx_bef->Fill(pre_m_D1);

                hvtxresX_D1->Fill(mcD1Vertex.X()-bestPos.X());
                hvtxresY_D1->Fill(mcD1Vertex.Y()-bestPos.Y());
                hvtxresZ_D1->Fill(mcD1Vertex.Z()-bestPos.Z());

            }
          }


          {
            int best_i=0;
            double best_chi2=1000;
            TCandidate *D2fit_best=0;
            TCandidate *k1fit_best, *p1fit_best, *p2fit_best, *k1nofit, *p1nofit, *p2nofit;
            TVector3 bestPos;
            Float_t D2vtx_mass, pre_m_D2;
            for (j=0;j<D2_vtx.GetLength();++j)
              {
                RhoKinVtxFitter vtxfitter(D2_vtx[j]);        // instantiate a vertex fitter
                vtxfitter.Fit();                          // do the vertex fit

                TCandidate *D2fit=vtxfitter.FittedCand(D2_vtx[j]);  // request the fitted Psi candidate
                TVector3 D2Vtx=D2fit->Pos();                    // and the decay vertex position
                double chi2_vtx_D2=vtxfitter.GlobalChi2();
                h_chi2_vtx_D2->Fill(chi2_vtx_D2/5); // Number degree of freedom 2N-3=5
                // plot mass and vtx x,y projection after fit
                hvpos_D2->Fill(D2Vtx.X(),D2Vtx.Y());
                hvzpos_D2->Fill(D2Vtx.Z());
                if(chi2_vtx_D2<best_chi2)
                {
                  best_chi2=chi2_vtx_D2;
                  best_i=l;
                  D2fit_best=D2fit;
                  p1fit_best=vtxfitter.FittedCand(*(D2fit_best->Daughter(0)));
                  p2fit_best=vtxfitter.FittedCand(*(D2fit_best->Daughter(1)));
                  k1fit_best=vtxfitter.FittedCand(*(D2fit_best->Daughter(2)));
                  p1nofit=D2_vtx[j].Daughter(0);
                  p2nofit=D2_vtx[j].Daughter(1);
                  k1nofit=D2_vtx[j].Daughter(2);
                  TLorentzVector tlvp1 = p1nofit->P4();
                  TLorentzVector tlvp2 = p2nofit->P4();
                  TLorentzVector tlvk1 = k1nofit->P4();
                  pre_m_D2= (tlvp1+tlvp2+tlvk1).M();

                  D2vtx_mass = D2fit_best->M();
                  bestPos = D2fit->Pos();
                }
              }
            if((psi.GetLength()!=0))
              {
                h_D2_vtx->Fill(D2vtx_mass);
                h_D2_vtx_bef->Fill(pre_m_D2);
                hvtxresX_D2->Fill(mcD2Vertex.X()-bestPos.X());
                hvtxresY_D2->Fill(mcD2Vertex.Y()-bestPos.Y());
                hvtxresZ_D2->Fill(mcD2Vertex.Z()-bestPos.Z());

            }
          }
        }

    }
  std::cout<<"Number of reconstructed eta_c = "<<n_reco<<std::endl;
  n_psi->SetBinContent(1,n_reco);


  out->cd();

  h_mD1_nocuts->Write();
  h_mD1_pid->Write();
  h_mD2_nocuts->Write();
  h_mD2_pid->Write();
  h_psi_nocut->Write();
  h_psi_pid->Write();


  n_psi->Write();
  n_events->Write();

  h_psi_Dmass->Write();
  h_psi_vtx->Write();
  h_D1_vtx_bef->Write();
  h_D2_vtx_bef->Write();
  h_D1_vtx->Write();
  h_D2_vtx->Write();

  h_psi_4c->Write();

  //h_mD_vtx->Write();
  h_mD_4c->Write();
  h_mD_final->Write();

  nc->Write();

  h_chi2_4c->Write();
  h_chi2b_4c->Write();
  h_chi2_vtx_D1->Write();
  hvzpos_D1->Write();
  hvpos_D1->Write();
  hvtxresX_D1->Write();
  hvtxresY_D1->Write();
  hvtxresZ_D1->Write();
  hvzpos_D2->Write();
  hvpos_D2->Write();
  hvtxresX_D2->Write();
  hvtxresY_D2->Write();
  hvtxresZ_D2->Write();


  out->Save();

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

  return 0;
}

void removeCombinatoric(TCandList &psi_list)
{
  TCandidate *phi1, *phi2, *k1, *k2, *k3, *k4;
  bool isOverlap;

  int len1=psi_list.GetLength();
  int ii=0;

  for (int i=0; i<len1; ++i)
    {
      phi1=psi_list[ii].Daughter(0);
      phi2=psi_list[ii].Daughter(1);
      k1=phi1->Daughter(0);
      k2=phi1->Daughter(1);
      k3=phi2->Daughter(0);
      k4=phi2->Daughter(1);

      if ((k1->IsCloneOf(*k3))||(k2->IsCloneOf(*k4)))
        {
          psi_list.Remove(psi_list[ii]);
          std::cout<<"Overlap found, i="<<i<<std::endl;
        }
      else
        ii++;
    }

}