StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
scan_embed_mudst.C
1 //Original Version done by Cristina May 2007
2 //Modification extra histograms by Xianglei
3 //Extra 3D Histograms Dca : ETA: pt and Nfit: Eta: pt are being added
4 
5 #ifndef __CINT__
6 #include "TROOT.h"
7 #include "TSystem.h"
8 #include <iostream.h>
9 #include "TH1.h"
10 #include "TH2.h"
11 #include "TH3.h"
12 #include "TFile.h"
13 #include "TTree.h"
14 #include "TChain.h"
15 #include "TTreeHelper.h"
16 #endif
17 
18 void scan_embed_mudst(
19 
20  const TString inputMuDst = "/star/data06/embed/andrewar/hiroshi/MuDst/*.MuDst.root",
21 
22  const Char_t* production = "P08ic",
23 
24  const Char_t* date = "041909",
25 
26  )
27 {
28  gSystem->Load("StarRoot");
29 
30  // -------------------- A n a l y s i s C u t s ---------------------------------
31 
32  const float VZCut = 30.; // Cut for Z Vertex
33  const float yCut = 1; // Cut for pseudoRapidity. we do not need this cuts for MC tracks and matched tracks
34  const float NFitCut = 25.; //Cut for Number of Fit Points
35  const int NCommCut = 10; // Cut for Number of common hits
36 
37  //------------------Ranges of the histograms ---------------------------------------
38 
39  const int nchNch = 200; const float maxNch = 1000;
40  const int nchPt = 80; const float maxPt = 8.0;
41  const int nchNhits = 50; const float minNhits = 0.; const float maxNhits = 50.;
42  const int nchDca = 100; const float maxDca = 3.;
43  const float maxDedx = 20.;
44  const float NSigmaCut = 2.0;
45  float mass2 = -1.0;
46 
47  //-------------------------MC---------------------------------
48 
49 
50  const Char_t* outputFileName = Form("../MuDst/%s_MuDst_%s.root", production, date);
51  TString oFile(outputFileName);
52 
53 
54  TString filePathR(inputMuDst);
55  cout << endl;
56  cout << endl;
57  cout << endl;
58  cout << " Input MuDst (real data) : " << filePathR.Data() <<endl;
59 
60  TTreeHelper THmu("MuDst");
61  THmu.AddFile(filePathR);
62 
63  TH1D* hMultR=new TH1D("hMultR","Multiplicity",1000,0,1000);
64 
65  TH1D* nChTotalR = new TH1D("nChTotalR", "nCharged Distribution", 1200, 0, maxNch);
66 
67  TH2F* dedxR = new TH2F("dedxR", "de/dx vs P", 400, 0.0, 2., 400, 0., maxDedx);
68  TH1D* hDcax = new TH1D("hDcax", "Distance of Closest Approach x", 1000, 0, maxDca);
69  TH1D* hDcay = new TH1D("hDcay", "Distance of Closest Approach y", 1000, 0, maxDca);
70  TH1D* hDcaz = new TH1D("hDcaz", "Distance of Closest Approach z", 1000, 0, maxDca);
71  TH3D* hDcaR = new TH3D("hDcaR","Nch:Pt:Dca",nchNch,0,maxNch,nchPt,0,maxPt,nchDca,0, maxDca);
72  TH3D* hNfitR = new TH3D("hNfitR","Nch:Pt:Nhits",nchNch,0,maxNch,nchPt,0,maxPt,nchNhits,minNhits,maxNhits);
73 
74 
75  TH1D* hNSigmaPi = new TH1D("hNSigmaPi","Sigma Pion",1000,-5000,5000);
76  TH1D* hNSigmaK = new TH1D("hNSigmaK","Sigma Kaon",1000,-5000,5000);
77  TH1D* hNSigmaP = new TH1D("hNSigmaP","Sigma Proton",1000,-5000,5000);
78  //Real DCA:Eta: pt and Nfit : Eta : pt
79 
80  TH3D* hPtEtaDcaR = new TH3D("hPtEtaDcaR","Pt:Eta:Dca", nchPt,0,maxPt,200,-1.8,1.8, nchDca,0, maxDca);// for Real
81  TH3D* hPtEtaNfitR = new TH3D("hPtEtaNfitR","Pt:Eta:Nfit",nchPt,0,maxPt,200,-1.8,1.8,nchNhits,0,maxNhits);// for Real
82  hPtEtaDcaR->Sumw2();
83  hPtEtaNfitR->Sumw2();
84 
85 
86 
87  const UShort_t *&nNeg= THmu("MuEvent.mRefMultNeg");
88  const UShort_t *&nPos= THmu("MuEvent.mRefMultPos");
89 
90  const Float_t *&VXR = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX1");
91  const Float_t *&VYR = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX2");
92  const Float_t *&VZR = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX3");
93 
94  const Int_t &ntrkR = THmu("GlobalTracks");
95  const UChar_t *&NFitR = THmu("GlobalTracks.mNHitsFit");
96  const Float_t *&dEdxR = THmu("GlobalTracks.mdEdx");
97  const Float_t *&PtR = THmu("GlobalTracks.mPt");
98  const Float_t *&EtaR = THmu("GlobalTracks.mEta");
99 
100  const Float_t *&dcax = THmu("GlobalTracks.mDCAGlobal.mX1");
101  const Float_t *&dcay = THmu("GlobalTracks.mDCAGlobal.mX2");
102  const Float_t *&dcaz = THmu("GlobalTracks.mDCAGlobal.mX3");
103  const Int_t *&vtxid = THmu("GlobalTracks.mVertexIndex");
104 
105  const UChar_t *&NHitsR = THmu("GlobalTracks.mNHits");
106  const Int_t *&NSigmaPion = THmu("GlobalTracks.mNSigmaPion");
107  const Int_t *&NSigmaKaon = THmu("GlobalTracks.mNSigmaKaon");
108  const Int_t *&NSigmaProton = THmu("GlobalTracks.mNSigmaProton");
109 
110 
111  /* if (id ==2 ||id == 8 || id == 9)
112  {const Int_t *&NSigma = THmu("GlobalTracks.mNSigmaPion"); }
113 
114  if (id == 11|| id ==12)
115  {const Int_t *&NSigma = THmu("GlobalTracks.mNSigmaKaon"); }
116 
117  if (id == 14|| id ==15)
118  {const Int_t *&NSigma = THmu("GlobalTracks.mNSigmaProton"); }
119  */
120 
121  int evReal = 0;
122  int evRealGood = 0;
123 
124  while(THmu.Next() )// && ev ==0)
125  {
126  hMultR->Fill(*nPos + *nNeg);
127  evReal++;
128 
129  if(fabs(*VZR)> VZCut) continue;
130 
131  evRealGood++;
132 
133  if( evReal % 100 == 0 ){
134  cout << Form("MuDst: (Good Event)/(Event number) = %5d/%5d, VZ = %1.4f, Number of primary tracks = %10d",
135  evRealGood, evReal, VZR, (Int_t)ntrkR)
136  << endl;
137  }
138 
139  for (Int_t itr=0; itr<ntrkR; itr++) {
140  float dca = sqrt(dcax[itr]*dcax[itr]+dcay[itr]*dcay[itr]+dcaz[itr]*dcaz[itr]);
141 
142  if (fabs(EtaR[itr]) > yCut) continue;
143 
144  float p = PtR[itr]*cosh(EtaR[itr]);
145 
146  hDcax->Fill(dcax[itr]);
147  hDcay->Fill(dcay[itr]);
148  hDcaz->Fill(dcaz[itr]);
149 
150  hNSigmaPi->Fill(NSigmaPion[itr]);
151  hNSigmaK->Fill(NSigmaKaon[itr]);
152  hNSigmaP->Fill(NSigmaProton[itr]);
153 
154  //filling DCA
155  //if(NHitsR[itr]>= NFitCut && fabs(1.*NSigma[itr]/1000.0) < NSigmaCut)
156  if(NFitR[itr]>= NFitCut && vtxid[itr]==0 && fabs(NSigmaPion[itr]/1000.0) < NSigmaCut)
157  {
158  hDcaR->Fill((*nNeg +*nPos), PtR[itr],dca);
159  hPtEtaDcaR->Fill( PtR[itr],EtaR[itr],dca);
160  }
161 
162  //filling Nfit
163  //if (dca < maxDca && fabs(NSigma[itr]/1000.0) < NSigmaCut && NHitsR[itr]>=10)
164  if (dca < maxDca && NFitR[itr]>=10 && vtxid[itr]==0 && fabs(NSigmaPion[itr]/1000.0) < NSigmaCut)
165  {
166  hNfitR ->Fill((*nNeg +*nPos), PtR[itr], NFitR[itr]);
167  hPtEtaNfitR ->Fill(PtR[itr], EtaR[itr], NFitR[itr]);
168  }
169 
170  //to fill de/dx vs P
171  //if(dca< maxDca && NFitR[itr]>= NFitCut && vtxid[itr]==0)
172  if(dca< maxDca && NFitR[itr]>= NFitCut)
173  //if(dca< maxDca && NHitsR[itr]>= NFitCut && fabs(EtaR[itr])< yCut)
174  dedxR->Fill(p,dEdxR[itr]*1e6);
175 
176  }//close for loop
177 
178  }//close while loop
179 
180 
181  //-------Writing-----------
182 
183  cout<< "Creating output file .... "<<oFile; flush(cout);
184 
185  TFile *fout=new TFile(oFile,"recreate");
186  fout->cd();
187 
188 
189  //-----REAL
190 
191  hMultR->Write();
192  nChTotalR->Write();
193  dedxR->Write();
194  hDcaR->Write();
195  hNfitR-> Write();
196 
197  hPtEtaNfitR->Write();
198  hPtEtaDcaR->Write();
199 
200  fout->GetList()->Sort();
201 
202  fout->Close();
203  cout<<" done"<<endl; flush(cout);
204 }
TDataSet::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Definition: TDataSet.cxx:893
Generated by   doxygen 1.8.5