StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Scanning_Embed.C
1 // $Id: Scanning_Embed.C,v 1.1 2007/04/27 15:00:30 cristina Exp $
3 // owner: Cristina
4 //
5 // What it does: scans files from MonteCarlo Embeddded Data and Ral ,fills the 3D
6 // histogram (nch,pt,dca), 3D histogram (nch,pt,nfit) and 2D histogram of dEdx vs P
7 // also fills the histograms to check some global variables such as z vertex, x-y vertex
8 // and Transverese momentum, Eta, phi and creates a root file as the output
9 //
10 // When running this macro, just include tag(id) for the particle that you are interested,
11 // check the cuts that you want to apply and the ranges for the histograms that you need.
12 // This macro uses Tree helper : "StMiniMcTre" and MuDst
14 
15 #ifndef __CINT__
16 #include "TROOT.h"
17 #include "TSystem.h"
18 #include <iostream.h>
19 #include "TH1.h"
20 #include "TH2.h"
21 #include "TH3.h"
22 #include "TFile.h"
23 #include "TTree.h"
24 #include "TChain.h"
25 #include "TTreeHelper.h"
26 #endif
27 
28 void scan_embed_mc(Int_t id=9){
29 
30 // -------------------- A n a l y s i s C u t s ---------------------------------
31 
32  float VZCut = 25; // Cut for Z Vertex
33  float yCut = 0.5; // Cut for Rapidity
34  float NFitCut = 25; //Cut for Number of Fit Points
35  int NCommCut= 10; // Cut for Number of common hits
36 
37  //------------------Ranges of the histograms ---------------------------------------
38 
39  int nchNch = 200; float maxNch=1000;
40 
41  int nchPt = 80; float maxPt = 15.0;
42 
43  int nchNhits= 50; float minNhits = 0.; float maxNhits = 50.;
44 
45  int nchDca = 100; float maxDca = 3.;
46 
47  float maxDedx = 20.;
48 
49  float mass2;
50 
51  if (id == 8) { TString tag = "PiPlus"; mass2 = 0.019;}
52  if (id == 9) { TString tag = "PiMinus"; mass2 = 0.019;}
53  if (id == 11) { TString tag = "KPlus"; mass2 = 0.245;}
54  if (id == 12) { TString tag = "KMinus"; mass2 = 0.245;}
55  if (id == 14) { TString tag = "Proton"; mass2 = 0.880;}
56  if (id == 15) { TString tag = "Pbar"; mass2 = 0.880;}
57  if (id == 50) { TString tag = "Phi"; mass2 = 1.020;}
58  if (id == 2) { TString tag = "Eplus"; mass2 = 0.511;}
59 
60 //location of files that you are scanning---------------------------------------
61 
62  TString filePathMC;
63 
64  filePath="~/Phi_101/Mini_Mc/*root";
65 
66  cout <<filePath<<endl;
67 
68 //oFile is where the output will be saved---------------------------------------------
69 
70  TString oFile;
71 
72  oFile="~/outputs/Phi/out_scan_embed_mc_"+tag+"1.root";
73 
74 //------------------------------------------
75 
76  TTreeHelper TH("StMiniMcTree");
77  TH.AddFile(filePath);
78 
79 //Definition of Histograms
80 
81 
82  TH1D* hMult=new TH1D("hMult","Multiplicity",1000,0,maxNch);
83 
84  TH1D* nChTotal = new TH1D("nChTotal", "nCharged Distribution", 1200, 0, maxNch);
85  TH2F* dedx = new TH2F("dedx", "de/dx vs P", 500, 0.0, maxPt, 400, 0., maxDedx);
86  TH3D* hDca = new TH3D("hDca3","Nch:Pt:Dca",nchNch,0,maxNch,nchPt,0,maxPt,nchDca,0, maxDca);
87  TH3D* hNfit = new TH3D("hNfit3","Nch:Pt:Nhits",nchNch,0,maxNch,nchPt,0,maxPt,nchNhits,0,maxNhits);
88 
89 //Defining the histograms to Check some Global Variables
90 
91  TH1D* vz = new TH1D("vz","Vertex Z",200,0.-50.,50.0);
92  TH2D* v_xy = new TH2D("vxy", "Vertez_X_Y", 200, -3.0, 3., 400, -3.0, 3.);
93 
94  TH1D* dvz = new TH1D("dvz","Delta Vertex Z MC - Reco",50,0.-5.,5.0);
95  TH1D* dvx = new TH1D("dvx","Delat Vertex X MC - Reco",50,0.-5.,5.0);
96  TH1D* dvy = new TH1D("dvy","Delta Vertex Y Mc - Reco",50,0.-5.,5.0);
97 
98  //Histograms for Embedded Tracks
99 
100  TH1D* hPtMc =new TH1D("PtMc","Pt_Mc",400,0,20);
101  TH1D* hEtaMc =new TH1D("EtaMc","Eta_Mc",200,-1.8,1.8);
102  TH1D* hPhiMc =new TH1D("PhiMc","Phi_Mc",200,-5.0,5.0);
103 
104 
105  //Histograms for Matched tracks
106 
107  TH1D* hPtPrMatched = new TH1D("PtPrMatched","Pt_Pr",400,0,20);
108  TH1D* hEtaPrMatched = new TH1D("EtaPrMatched","EtaPr",200,-1.8,1.8);
109  TH1D* hPhiPrMatched = new TH1D("PhiPrMatched","PhiPr",200,-5.0,5.0);
110 
111  //Histogram for Energy Loss
112 
113  TH2D* hPtM_E_Pr =new TH2D("hPtM_E_Pr",title,100,0,10,200,-0.1,0.1);// for Primary Tracks
114  TH2D* hPtM_E_Gl =new TH2D("hPtM_E_Gl",title,100,0,10,200,-0.1,0.1);// for Global Tracks
115 
116 // Definition of variables-----------------------------------------
117 
118  const Int_t &nch = TH("mNUncorrectedPrimaries");
119  const Float_t &VX = TH("mVertexX");
120  const Float_t &VY = TH("mVertexY");
121  const Float_t &VZ = TH("mVertexZ");
122 
123  const Float_t &VXmc = TH("mMcVertexX");
124  const Float_t &VYmc = TH("mMcVertexY");
125  const Float_t &VZmc = TH("mMcVertexZ");
126 
127 
128  // Branch : McTracks // branch is related with the embedded tracks
129 
130  const Int_t &ntrk1 = TH ("mMcTracks");
131 
132  const Float_t *&PtMc = TH("mMcTracks.mPtMc");
133  const Float_t *&PzMc = TH("mMcTracks.mPzMc");
134  const Float_t *&NHits = TH("mMcTracks.mNHitMc");
135  const Float_t *&EtaMc = TH("mMcTracks.mEtaMc");
136  const Float_t *&PhiMc = TH("mMcTracks.mPhiMc");
137 
138  // Branch : Matched Pairs //
139 
140  //Extension Pr ( Pr stands for Primary and it's related with all possible reconstructed tracks)
141  //Ext Mc ( Mc stands for MonteCarlo and it's related with reconstructed and matched with embedded tracks)
142 
143  const Int_t &ntrk = TH("mMatchedPairs");
144 
145  const Float_t *&PtPrMatched = TH("mMatchedPairs.mPtPr");
146  const Float_t *&PzPrMatched = TH("mMatchedPairs.mPzPr");
147  const Float_t *&EtaPrMatched = TH("mMatchedPairs.mEtaPr");
148  const Float_t *&PhiPrMatched= TH("mMatchedPairs.mPhiPr");
149 
150  const Float_t *&PtMcMatched = TH("mMatchedPairs.mPtMc");
151  const Float_t *&PzMcMatched = TH("mMatchedPairs.mPzMc");
152  const Float_t *&EtaMcMatched = TH("mMatchedPairs.mEtaMc");
153  const Float_t *&PhiMcMatched= TH("mMatchedPairs.mPhiMc");
154 
155  const Float_t *&PtGlMatched= TH("mMatchedPairs.mPtGl");
156  const Float_t *&PzGlMatched= TH("mMatchedPairs.mPzGl");
157  const Float_t *&EtaGlMatched = TH("mMatchedPairs.mEtaMGl");
158  const Float_t *&PhiGlMatched = TH("mMatchedPairs.mPhiMGl");
159 
160 
161  const Float_t *&dEdx = TH("mMatchedPairs.mDedx");
162  const Float_t *&Dca = TH("mMatchedPairs.mDcaGl");
163  const Short_t *&NFit = TH("mMatchedPairs.mFitPts");
164  const Short_t *&NComm = TH("mMatchedPairs.mNCommonHit");
165 
166  const Short_t *&GidMatched = TH("mMatchedPairs.mGeantId");
167 
168  //========================================================
169  // Branch : MatGlobalPairs
170  /*
171  const Int_t &ntrk = TH("mMatGlobPairs");
172 
173  const Float_t *&PtGlGlobal = TH("mMatGlobPairs.mPtGl");
174  const Float_t *&PzGlGlobal = TH("mMatGlobPairs.mPzGl");
175  const Float_t *&EtaGlGlobal = TH("mMatGlobPairs.mEtaGl");
176  const Float_t *&PhiGlGlobal = TH("mMatGlobPairs.mPhiGl");
177 
178  const Float_t *&PtPrGlobal= TH("mMatGlobPairs.mPtPr");
179  const Float_t *&PzPrGlobal= TH("mMatGlobPairs.mPzPr );
180  const Float_t *&EtaPrGlobal = TH("mMatGlobPairs.mEtaPr");
181  const Float_t *&PhiPrGlobal = TH("mMatGlobPairs.mPhiPr"
182 
183  const Float_t *&dEdx = TH("mMatGlobPairs.mDedx");
184  const Float_t *&Dca = TH("mMatGlobPairs.mDcaGl");
185  const Short_t *&NFit = TH("mMatGlobPairs.mFitPts");
186  const Short_t *&NComm = TH("mMatGlobPairs.mNCommonHit");
187 
188  const Short_t *&GidGlobal = TH("mMatGlobPairs.mGeantId");
189  */
190  //=======================================================
191  /*
192  // Branch : Ghost Pairs (Real Tracks with no matches)
193 
194  const Int_t &ntrk = TH("mGhostPairs");
195 
196  const Float_t *&PtGlGlobal = TH("mGhostPairs.mPtGl");
197  const Float_t *&PzGlGlobal = TH("mGhostPairs.mPzGl");
198  const Float_t *&EtaGlGlobal = TH("mGhostPairs.mEtaGl");
199  const Float_t *&PhiGlGlobal = TH("mGhostPairs.mPhiGl");
200 
201  const Float_t *&PtPrGlobal= TH("mGhostPairs.mPtPr");
202  const Float_t *&PzPrGlobal= TH("mGhostPairs.mPzPr" );
203  const Float_t *&EtaPrGlobal = TH("mGhostPairs.mEtaPr");
204  const Float_t *&PhiPrGlobal = TH("MmGhostPairs.mPhiPr"
205 
206  const Float_t *&dEdx = TH("mGhostPairs.mDedx");
207  const Float_t *&Dca = TH("mGhostPairs.mDcaGl");
208  const Short_t *&NFit = TH("mGhostPairs.mFitPts");
209  const Short_t *&NComm = TH("mGhostPairs.mNCommonHit");
210 
211  const Short_t *&GidGhost = TH("mGhostPairs.mGeantId");
212  */
213 //-----------------------------------------------------------------
214 
215 int ev = 0;
216 while(TH.Next())// && ev ==0)
217 {
218  //ev++;
219  vz->Fill(VZ);// 1 vertex per event
220  // first fill vz and then apply the cuts
221 
222  vxy->Fill(VX,VY);
223 
224  dvz->Fill(VZ-VZmc);
225  dvx->Fill(VX-VXmc);
226  dvy->Fill(VY-VYmc);
227 
228  if(fabs(VZ)> VZCut) continue;
229  ev++;
230 
231 for (Int_t itr=0; itr<ntrk; itr++)
232  {
233  //to fill the histograms related with all possible reconstructed tracks
234  hMult ->Fill();
235 
236  if ( GidMatched[itr]!=11 && GidMatched[itr]!=12) continue; //use this line just for Phi
237 
238  if (fabs(EtaPrMatched[itr]) > yCut) continue;
239 
240  //Defining Momentum
241 
242  float p = sqrt(PtPrMatched[itr]*PtPrMatched[itr]+PzPrMatched[itr]*PzPrMatched[itr]);
243 
244 
245  if (NFit[itr] >= NFitCut && NComm[itr] > NCommCut)
246  hDca->Fill(nch, PtPrMatched[itr],Dca[itr]);
247 
248  if (Dca[itr] < maxDca && (1.*(NFit[itr])/NComm[itr]) < 2.5)
249  hNfit ->Fill(nch, PtPrMatched[itr], NFit[itr]);
250 
251  if( Dca[itr] < maxDca && NFit[itr] >= NFitCut && fabs(EtaPrMatched[itr])< yCut)
252  dedx->Fill(p,dEdx[itr]*1e6);
253 
254  //to fill the histograms related with the global variables
255 
256  if(NFit[itr] >= NFitCut && NComm[itr] > NCommCut && Dca[itr] < maxDca )
257  {
258  v_xy->Fill(VX[itr], VY[itr]);
259  hPtPrMatched->Fill(PtPrMatched[itr]);
260  hEtaPrMatched->Fill(EtaPrMatched[itr]);
261  hPhiPrMatched->Fill(PhiPrMatched[itr]);
262  }
263  }//close for
264 
265 // To Fill histograms of embedded tracks
266 
267  for (int itr1=0; itr1 < ntrk1; itr1++)
268  {
269  if (fabs(EtaM[itr]) > yCut) continue;
270 
271  hPtMc->Fill(PtMc[itr1]);
272  hEtaMc->Fill(EtaMc[itr]);
273  hPhiMc->Fill(PhiMc[itr]);
274  }
275 
276 }//close while
277 
278 
279 //--------------------------SCANNING REAL DATA------------------------------
280 
281 //Using TTree Helper MuDst
282 
283 TTreeHelper THmu("MuDst");
284 
285 THmu.AddFile(filePathR);
286 
287 
288 TH1D* nChTotalR = new TH1D("nChTotalR", "nCharged Distribution", 1200, 0, 1200);
289 TH2F* dedxR = new TH2F("dedxR", "de/dx vs P", 500, 0.0, 4., 400, 0., 20);//start at 0
290 TH1D* hDcax=new TH1D("hDcax", "Distance of Closest Approach x", 1000, 0, 5);
291 TH1D* hDcay=new TH1D("hDcay", "Distance of Closest Approach y", 1000, 0, 5);
292 TH1D* hDcaz=new TH1D("hDcaz", "Distance of Closest Approach z", 1000, 0, 5);
293 TH3D* hDcaR=new TH3D("hDcaR","Nch:Pt:Dca",nchNch,0,maxNch,nchPt,0,maxPt,nchDca,0, maxDca);
294 TH3D* hNfitR = new TH3D("hNfitR","Nch:Pt:Nhits",nchNch,0,maxNch,nchPt,0,maxPt,nchNhits,minNhits,maxNhits);
295 TH1D* hNSigmaPi=new TH1D("hNSigmaPi","Sigma Pion",1000,0,1000);
296 TH1D* hNSigmaK=new TH1D("hNSigmaK","Sigma Kaon",1000,0,1000);
297 TH1D* hNSigmaP=new TH1D("hNSigmaP","Sigma Proton",1000,0,1000);
298 
299 
300  const UShort_t *&nNeg= THmu("MuEvent.mRefMultNeg");
301  const UShort_t *&nPos= THmu("MuEvent.mRefMultPos");
302 
303  const Float_t *&VX1 = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX1");
304  const Float_t *&VY1 = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX2");
305  const Float_t *&VZ1 = THmu("MuEvent.mEventSummary.mPrimaryVertexPos.mX3");
306 
307  const Int_t &ntrk1 = THmu("PrimaryTracks");
308 
309  const Float_t *&dEdxR = THmu("PrimaryTracks.mdEdx");
310  const Float_t *&Pt1 = THmu("PrimaryTracks.mPt");
311  const Float_t *&Eta1 = THmu("PrimaryTracks.mEta");
312 
313  const Float_t *&dcax1 = THmu("PrimaryTracks.mDCAGlobal.mX1");
314  const Float_t *&dcay1 = THmu("PrimaryTracks.mDCAGlobal.mX2");
315  const Float_t *&dcaz1 = THmu("PrimaryTracks.mDCAGlobal.mX3");
316 
317  const UChar_t *&NHits1 = THmu("PrimaryTracks.mNHits");
318  const Int_t *&NSigmaPion = THmu("PrimaryTracks.mNSigmaPion");
319  const Int_t *&NSigmaKaon = THmu("PrimaryTracks.mNSigmaKaon");
320  const Int_t *&NSigmaProton = THmu("PrimaryTracks.mNSigmaProton");
321 
322 
323  if (id == 8 || id == 9)
324  {const Int_t *&NSigma = THmu("PrimaryTracks.mNSigmaPion"); }
325 
326  if (id == 11|| id ==12)
327  {const Int_t *&NSigma = THmu("PrimaryTracks.mNSigmaKaon"); }
328 
329  if (id == 14|| id ==15)
330  {const Int_t *&NSigma = THmu("PrimaryTracks.mNSigmaProton"); }
331 
332 //const Int_t &nchR = *nNeg + *nPos;
333 //THmu.Print();
334 //return;
335 
336 int ev = 0;
337 while(THmu.Next() )// && ev ==0)
338 {
339 hMult->Fill(*nPos + *nNeg);
340 
341 for (Int_t itr=0; itr<ntrk1; itr++)
342  {
343  ev++;
344 
345 //cout<<*VZ1<<endl;
346 
347  if(fabs(*VZ1)>25.) continue;//cut in vertex
348 
349  float dca = sqrt(dcax1[itr]*dcax1[itr]+dcay1[itr]*dcay1[itr]+dcaz1[itr]*dcaz1[itr]);
350  float p_r = Pt1[itr]*cosh(Eta1[itr]);
351 
352  hDcax->Fill(dcax1[itr]);
353  hDcay->Fill(dcay1[itr]);
354  hDcaz->Fill(dcaz1[itr]);
355 
356  hNSigmaPi->Fill(NSigmaPion[itr]);
357  hNSigmaK->Fill(NSigmaKaon[itr]);
358  hNSigmaP->Fill(NSigmaProton[itr]);
359 
360  //to fill DCA histogram
361 
362  if (fabs(Eta1[itr]) >0.5) continue;
363 
364  if(NHits1[itr]>= 25)
365  hDcaR->Fill((*nNeg +*nPos), Pt1[itr],dca);
366 
367 
368  //to fill Nfit points histogram
369 
370  if (dca < 3)
371  hNfitR ->Fill(*nNeg +*nPos, Pt1[itr], NHits1[itr]);
372 
373 
374  //to fill de/dx vs P
375 
376  if(dca< 3. && NHits1[itr]>= 25 && fabs(Eta1[itr])< 0.5 && NSigma[itr]/1000 < 2.)
377  dedxR->Fill(p_r,dEdxR[itr]*1e6);
378 
379 }
380 }//close MuDst
381 
382 
383 cout<< "Creating output file .... "<<oFile; flush(cout);
384 
385  TFile *fout=new TFile(oFile,"recreate");
386  fout->cd();
387 
388  //writting all the histograms
389 
390  nChTotal->Write();
391  dedx->Write("dEdx");
392  hNfit->Write();
393  hDca-> Write();
394 
395  //Global Variables
396 
397  v_xy->Write();
398  vz->Write();
399 
400  dvx->Write();
401  dvy->Write();
402  dvz->Write();
403 
404  hPtPrMatched->Write();
405  hEtaPrMatched->Write();
406  hPhiPrMatched->Write();
407 
408  hEtaMc->Write();
409  hPtMc->Write();
410  hPhiMc->Write();
411 
412  fout->Close();
413  cout<<" done"<<endl; flush(cout);
414 }
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