EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
analysis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file analysis.C
1 
2 #define _NDF_MAX_ 1000
3 
4 void analysis()
5 {
6  // Load basic libraries;
7  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
8 
9  // Input simulated & reconstructed files;
10  TFile *ff = new TFile("simulation.root");
11  TTree *cbmsim = ff->Get("cbmsim");
12  cbmsim->AddFriend("cbmsim", "reconstruction.root");
13 
14  // Figure out max and most popular ndf value;
15  TClonesArray *rctrk = new TClonesArray("PndPidCandidate");
16  cbmsim->SetBranchAddress("PidChargedCand",&rctrk);
17  unsigned nEvt = cbmsim->GetEntries(), ndfMax = 0;
18  unsigned arr[1000]; memset(arr, 0x00, sizeof(arr));
19  for (unsigned evt = 0; evt<nEvt; evt++) {
20  cbmsim->GetEntry(evt);
21 
22  if (rctrk->GetEntriesFast()) {
23  PndPidCandidate *rctrack = rctrk->At(0);
24 
25  int ndf = rctrack->GetDegreesOfFreedom();
26  if (ndf > ndfMax) ndfMax = ndf;
27  if (ndf < _NDF_MAX_) arr[ndf]++;
28  }
29  }
30  unsigned ndfMostPopular = 0, ndfMostPopularStat = 0;
31  for(unsigned iq=0; iq<_NDF_MAX_; iq++)
32  if (arr[iq] > ndfMostPopularStat) {
33  ndfMostPopular = iq;
34  ndfMostPopularStat = arr[iq];
35  }
36 
37  char *expression = "100.*(PidChargedCand.GetMomentum().Mag()-MCTrack.GetMomentum().Mag())/MCTrack.GetMomentum().Mag()";
38  //char *expression = "10000.*(PidChargedCand.GetMomentum().Theta()-MCTrack.GetMomentum().Theta())";
39  //char *expression = "1000.*(PidChargedCand.GetMomentum().Phi()-MCTrack.GetMomentum().Phi())";
40  char cut[1024];
41  // Allow to lose 3 degrees of freedom compared to the max (assume standard) case;
42  sprintf(cut, "EicIdealGenTrack.fNDF>=%d&&EicIdealGenTrack.fChiSquareCCDF>.001", ndfMostPopular-2);
43 
44  double par[100]; memset(par, 0x00, sizeof(par));
45 
46  TH1D *dp1 = new TH1D("dp1", "dp1", 100, -20.0, 20.0);
47  cbmsim->Project("dp1", expression, cut);
48  TF1 *fq1 = new TF1("fq1", "gaus(0)", -20.0, 20.0);
49  par[0] = 100.0; par[1] = 0.0; par[2] = 1.0;
50  fq1->SetParameters(par);
51  dp1->Fit("fq1","R");
52  fq1->GetParameters(par);
53 
54  double sigma2 = fabs(par[2]), min2 = -5*sigma2, max2 = 5*sigma2;
55 
56  TH1D *dp2 = new TH1D("dp2", "dp2", 100, min2, max2);
57  cbmsim->Project("dp2", expression, cut);
58  TF1 *fq2 = new TF1("fq2", "gaus(0)", min2, max2);
59  fq2->SetParameters(par);
60  dp2->Fit("fq2","R");
61  fq2->GetParameters(par);
62 
63  double sigma3 = fabs(par[2]), min3 = -5*sigma3, max3 = 5*sigma3;
64 
65  TH1D *dp3 = new TH1D("dp3", "dp3", 100, min3, max3);
66  cbmsim->Project("dp3", expression, cut);
67  TF1 *fq3 = new TF1("fq3", "gaus(0)", min3, max3);
68  fq3->SetParameters(par);
69  dp3->Fit("fq3","R");
70 } // analysis()