StRoot: macros/embedding/scan_embed

00001 //Original Version done by Cristina May 2007
00002 //Modification extra histograms by Xianglei
00003 //Extra 3D Histograms Dca : ETA: pt and Nfit: Eta: pt are being added
00004 
00005 #ifndef __CINT__
00006 #include "TROOT.h"
00007 #include "TSystem.h"
00008 #include <iostream.h>
00009 #include "TH1.h"
00010 #include "TH2.h"
00011 #include "TH3.h"
00012 #include "TFile.h"
00013 #include "TTree.h"
00014 #include "TChain.h"
00015 #include "TTreeHelper.h"
00016 #endif
00017 
00018 void scan_embed_mc(
00019                 const Int_t id = 8,
00020                 
00021                 const TString inputMiniMc ="/eliza3/starprod/embedding/P08id/dAu/*/*/*.minimc.root", //test 
00022                 
00023                 const Char_t* particle= "Phi",
00024                 const Char_t* production = "P08id",
00025                 const Int_t eid = 90117,
00026                 const Char_t* date = "040709",
00027         
00028                 )
00029 {
00030   gSystem->Load("StarRoot");
00031   
00032   // -------------------- A n a l y s i s  C u t s  ---------------------------------
00033   
00034   const float VZCut   = 30.; // Cut for Z Vertex 
00035   const float yCut    = 1;  // Cut for pseudoRapidity. we do not need this cuts for MC tracks and matched tracks
00036   const float NFitCut = 25.; //Cut for Number of Fit Points
00037   const int  NCommCut = 10;  // Cut for Number of common hits
00038   
00039   //------------------Ranges of the histograms  ---------------------------------------
00040   
00041   const int  nchNch     = 200;   const float maxNch   = 1000;
00042   const int  nchPt      = 80;    const float maxPt    = 8.0;
00043   const int  nchNhits   = 50;    const float minNhits = 0.;    const float maxNhits = 50.;   
00044   const int  nchDca     = 100;   const float maxDca   = 3.;
00045   const float maxDedx   = 20.;
00046   const float NSigmaCut = 2.0;
00047   float mass2 = -1.0;
00048   
00049   if (id == 2)  { TString tag = "Eplus";  mass2 = 0.511; Char_t *Daugther = "Eplus";}
00050   if (id == 8)  { TString tag = "PiPlus";  mass2 = 0.019; Char_t *Daugther = "PiPlus";}
00051   if (id == 9)  { TString tag = "PiMinus"; mass2 = 0.019; Char_t *Daugther = "PiMius";}
00052   if (id == 11) { TString tag = "KPlus";  mass2 = 0.245; Char_t *Daugther = "KPlus";}
00053   if (id == 12) { TString tag = "KMinus";  mass2 = 0.245; Char_t *Daugther = "KMinus";}
00054   if (id == 14) { TString tag = "Proton";  mass2 = 0.880; Char_t *Daugther = "Proton";}
00055   if (id == 15) { TString tag = "Pbar";    mass2 = 0.880; Char_t *Daugther = "Pbar";}
00056   if (id == 50) { TString tag = "Phi";    mass2 = 1.020; }
00057   
00058   
00059   if( mass2 < 0 ) {
00060     cout << "Negative mass2 for pid = " << id << endl;
00061     return;
00062   }
00063   
00064   //-------------------------MC---------------------------------
00065   
00066   TString filePathMc(inputMiniMc);
00067   
00068   cout << endl;
00069   cout << endl;
00070   cout << "Input MiniMc file : " << filePathMc.Data() <<endl;
00071   cout << "#---------------------------------------------------------------------------" << endl;
00072   cout << " Event selection" << endl;
00073   cout << " |VZ| < " << VZCut << endl;
00074   cout << "#---------------------------------------------------------------------------" << endl;
00075   cout << " Track selection (mMatchedPairs)" << endl;
00076   cout << "  For DCA          : mMatchedPairs.mFitPts >= " << NFitCut << " && mMatchedPairs.mNCommonHit > " << NCommCut << endl;
00077   cout << "  For NFit         : mMatchedPairs.mDcaGl < " << maxDca << " && mMatchedPairs.mFitPts/mMatchedPairs.mNCommonHit < 2.5 && mMatchedPairs.mFitPts >= 10" << endl;
00078   cout << "  For dE/dx        : mMatchedPairs.mDcaGl < " << maxDca << " && mMatchedPairs.mFitPts >= " << NFitCut << " && |mMatchedPairs.mEtaPr| < " << yCut << endl;
00079   cout << "  For pt, eta, phi : |mMatchedPairs.mDcaGl| < " << maxDca << endl;
00080   cout << "#---------------------------------------------------------------------------" << endl;
00081   cout << " Track selection (mGhostPairs)" << endl;
00082   cout << "  For dE/dx        : mGhostPairs.mDcaGl < " << maxDca << " && mGhostPairs.mFitPts >= " << NFitCut << " && |mGhostPairs.mEtaPr| < " << yCut << endl;
00083   cout << "#---------------------------------------------------------------------------" << endl;
00084   cout << " Track selection (mMcTracks)" << endl;
00085   cout << "  For pt, eta, phi : |mMatchedPairs.mEtaPr| > " << yCut << endl;
00086   cout << endl;
00087 
00088   const Char_t* outputFileName = Form("./%s_%s_%s_mc_%s.root", production, particle, Daugther, date);
00089   TString oFile(outputFileName);
00090 
00091   TTreeHelper TH("StMiniMcTree"); 
00092   TH.AddFile(filePathMc);
00093 
00094   //All reconstructed Histograms  
00095 
00096   TH1D* hMult=new TH1D("hMult","Multiplicity",1000,0,maxNch); 
00097 
00098   TH1D* nChTotal = new TH1D("nChTotal", "nCharged Distribution", 1200, 0, maxNch);
00099   TH2F* dedx = new TH2F("dedx", "de/dx vs P", 400, 0.0, 2., 400, 0., maxDedx);
00100   TH3D* hDca = new TH3D("hDca","Nch:Pt:Dca",nchNch,0,maxNch,nchPt,0,maxPt,nchDca,0, maxDca);
00101   TH3D* hNfit = new TH3D("hNfit","Nch:Pt:Nhits",nchNch,0,maxNch,nchPt,0,maxPt,nchNhits,0,maxNhits);
00102 
00103   // Checking some Global Variables
00104 
00105   TH1D* vz = new TH1D("vz","Vertex Z",200,0.-50.,50.0);
00106   TH2D* v_xy =  new TH2D("vxy", "Vertez_X_Y", 200, -3.0, 3., 400, -3.0,  3.);
00107 
00108   TH1D* dvz = new TH1D("dvz","Delta Vertex Z MC - Reco",50,0.-5.,5.0);
00109   TH1D* dvx = new TH1D("dvx","Delta Vertex X MC -  Reco",50,0.-5.,5.0);
00110   TH1D* dvy = new TH1D("dvy","Delta Vertex Y MC - Reco",50,0.-5.,5.0);
00111 
00112   TH1D* hPtM = new TH1D("PtM","Pt_Pr",400,0,20);
00113   TH1D* hEtaM = new TH1D("EtaM","EtaPr",200,-1.8,1.8);
00114   TH1D* hPhiM = new TH1D("PhiM","PhiPr",200,-5.0,5.0);
00115   TH1D* hYM = new TH1D("YM","YPr",100,-1.5,1.5);
00116   TH3D* hPtYPhiM = new TH3D("PtYPhiM","Pt:Y:Phi",nchPt,0,maxPt,100,-1.5,1.5,200,-5.0,5.0);
00117   TH1D* hFSectM = new TH1D("FSectM","First Sector",30,0,30.);
00118   TH1D* hLSectM = new TH1D("LSectM","Last Sector",30,0,30.);
00119 
00120   TH1D* hPtMc =new TH1D("PtMc","Pt_Mc",400,0,20);
00121   TH1D* hEtaMc=new TH1D("EtaMc","Eta_Mc",200,-1.8,1.8);
00122   TH1D* hYMc=new TH1D("YMc","Y_Mc",100,-1.5,1.5);
00123   TH1D* hPhiMc=new TH1D("PhiMc","Phi_Mc",200,-5.0,5.0);
00124   TH3D* hPtYPhiMc = new TH3D("PtYPhiMc","Pt:Y:Phi",nchPt,0,maxPt,100,-1.5,1.5,200,-5.0,5.0);
00125 
00126   hPtMc->Sumw2();
00127   hPtM->Sumw2();
00128 
00129   //Energy loss
00130   TString title = Form("ptM-ptMc :ptM ( |vtx|<30, %.1f< |yMc|<%.1f, nFitM>=25, nComm>=10)", 1., 1.);
00131   TH2D* hPtM_E = new  TH2D ("hPtM_E", title.Data(), 50, 0, 10, 50, -0.05, 0.05);
00132   TString title = Form("YM-YMc :YM ( |vtx|<30, %.1f< |yMc|<%.1f, nFitM>=25, nComm>=10)", 1., 1.);
00133   TH2D* hYM_E = new  TH2D ("hYM_E", title.Data(), 50, -1.2, 1.2, 50, -0.05, 0.05);
00134   TString title = Form("PhiM-PhiMc :PhiM ( |vtx|<30, %.1f< |yMc|<%.1f, nFitM>=25, nComm>=10)", 1., 1.);
00135   TH2D* hPhiM_E = new  TH2D ("hPhiM_E", title.Data(), 50, -3.2, 3.2, 50, -0.05, 0.05);
00136 
00137   //Ghost ==Real 
00138   TH2F* dedxG = new TH2F("dedxG","dE/dx vs p", 400,0., 2., 400, 0.,maxDedx);//for MIPS
00139   TH3D* hDcaG = new TH3D("hDcaG","Nch:Pt:Dca",nchNch,0,maxNch,nchPt,0,maxPt,nchDca,0, maxDca);
00140   TH3D* hNfitG = new TH3D("hNfitG","Nch:Pt:Nhits",nchNch,0,maxNch,nchPt,0,maxPt,nchNhits,0,maxNhits);
00141   TH1D* hPhiG = new TH1D("PhiG","PhiG",200,-5.0,5.0);
00142   TH1D* hFSectG = new TH1D("FSectG","First Sector",30,0,30.);
00143   TH1D* hLSectG = new TH1D("LSectG","Last Sector",30,0,30.);
00144 
00145   //3 D histograms for DCa ETA  and Pt, also NFit Eta and Pt
00146 
00147   TH3D* hPtEtaDcaM  = new TH3D("hPtEtaDcaM","Pt:Eta:DcaM", nchPt,0, maxPt,200,-1.8,1.8, nchDca,0, maxDca);// for MonteCarlo
00148   TH3D* hPtEtaDcaG  = new TH3D("hPtEtaDcaG","Pt:Eta:DcaG", nchPt,0, maxPt,200,-1.8,1.8, nchDca,0, maxDca);// for Global
00149   hPtEtaDcaM->Sumw2();
00150   hPtEtaDcaG->Sumw2();
00151 
00152   TH3D* hPtEtaNfitM  = new TH3D("hPtEtaNfitM","Pt:Eta:NfitM",nchPt,0,maxPt,200,-1.8,1.8,nchNhits,0,maxNhits);// for MonteCarlo
00153   TH3D* hPtEtaNfitG  = new TH3D("hPtEtaNfitG","Pt:Eta:NfitG",nchPt,0,maxPt,200,-1.8,1.8,nchNhits,0,maxNhits);// for Global
00154   hPtEtaNfitM->Sumw2();
00155   hPtEtaNfitG->Sumw2();
00156 
00157   //Variables
00158   
00159   const Int_t  &nch = TH("mNUncorrectedPrimaries"); 
00160   const Float_t &VX  = TH("mVertexX");
00161   const Float_t &VY  = TH("mVertexY");
00162   const Float_t &VZ  = TH("mVertexZ");
00163 
00164   const Float_t &VXmc  = TH("mMcVertexX");
00165   const Float_t &VYmc  = TH("mMcVertexY");
00166   const Float_t &VZmc  = TH("mMcVertexZ");
00167   
00168   //Mc Tracks 
00169   
00170   const Int_t  &ntrk1 = TH ("mMcTracks"); 
00171   
00172   const Float_t *&PtMc  = TH("mMcTracks.mPtMc");
00173   const Float_t *&PzMc  = TH("mMcTracks.mPzMc");
00174   const Float_t *&NHits  = TH("mMcTracks.mNHitMc");
00175   const Float_t *&EtaMc  = TH("mMcTracks.mEtaMc");
00176   const Float_t *&PhiMc  = TH("mMcTracks.mPhiMc");
00177 
00178   // Matched Pairs 
00179   /*
00180   const Int_t  &ntrk  = TH("mMatchedPairs"); //Extension Pr - all possible reconstructed tracks
00181   
00182   const Float_t *&PtM  = TH("mMatchedPairs.mPtPr");
00183   const Float_t *&PzM  = TH("mMatchedPairs.mPzPr");
00184   const Float_t *&EtaM = TH("mMatchedPairs.mEtaPr");
00185   const Float_t *&PhiM = TH("mMatchedPairs.mPhiPr");
00186   const Short_t *&FSectM = TH("mMatchedPairs.mFirstSector");
00187   const Short_t *&LSectM = TH("mMatchedPairs.mLastSector");
00188 
00189 
00190   //Ext  Mc -  reconstructed and matched with embedded tracks
00191 
00192   const Float_t *&PtMmc= TH("mMatchedPairs.mPtMc");
00193   const Float_t *&PzMmc= TH("mMatchedPairs.mPzMc");
00194   const Float_t *&EtaMmc = TH("mMatchedPairs.mEtaMc");
00195   const Float_t *&PhiMmc = TH("mMatchedPairs.mPhiMc");
00196   
00197   const Float_t *&dEdx  = TH("mMatchedPairs.mDedx");
00198   const Float_t *&Dca  = TH("mMatchedPairs.mDcaGl");
00199   const Short_t *&NFit  = TH("mMatchedPairs.mFitPts");
00200   const Short_t *&NComm = TH("mMatchedPairs.mNCommonHit");
00201 
00202   const Short_t *&Gid = TH("mMatchedPairs.mGeantId");
00203 
00204   */
00205   //========================================================
00206   //to scan meson Phi, use MatGlobalPairs instead of Matchd pairs
00207   
00208      const Int_t  &ntrk  = TH("mMatGlobPairs");  
00209 
00210      const Float_t *&PtM  = TH("mMatGlobPairs.mPtGl");
00211      const Float_t *&PzM  = TH("mMatGlobPairs.mPzGl");
00212      const Float_t *&EtaM = TH("mMatGlobPairs.mEtaGl");
00213      const Float_t *&PhiM = TH("mMatGlobPairs.mPhiGl");
00214 
00215      const Short_t *&FSectM = TH("mMatGlobPairs.mFirstSector");
00216      const Short_t *&LSectM = TH("mMatGlobPairs.mLastSector");
00217   
00218 
00219      const Float_t *&PtMmc= TH("mMatGlobPairs.mPtMc");
00220      const Float_t *&PzMmc= TH("mMatGlobPairs.mPzMc");
00221      const Float_t *&EtaMmc = TH("mMatGlobPairs.mEtaMc");
00222      const Float_t *&PhiMmc = TH("mMatGlobPairs.mPhiMc");
00223 
00224      const Float_t *&dEdx = TH("mMatGlobPairs.mDedx");
00225      const Float_t *&Dca  = TH("mMatGlobPairs.mDcaGl");
00226      const Short_t *&NFit  = TH("mMatGlobPairs.mFitPts");
00227      const Short_t *&NComm = TH("mMatGlobPairs.mNCommonHit");
00228 
00229      const Short_t *&Gid = TH("mMatGlobPairs.mGeantId");
00230 
00231   //=======================================================
00232    
00233 
00234   //Ghost Tracks
00235 
00236   const Int_t    &ntrkG  = TH("mGhostPairs");  
00237 
00238   const Float_t  *&PtG  = TH("mGhostPairs.mPtPr");
00239   const Float_t  *&PzG  = TH("mGhostPairs.mPzPr");
00240   const Float_t  *&dEdxG= TH("mGhostPairs.mDedx");
00241   const Float_t  *&DcaG = TH("mGhostPairs.mDcaGl");
00242   const Float_t  *&EtaG = TH("mGhostPairs.mEtaPr");
00243   const Float_t  *&PhiG = TH("mGhostPairs.mPhiPr");
00244   const Short_t  *&NFitG = TH("mGhostPairs.mFitPts");
00245   const Short_t *&FSectG = TH("mGhostPairs.mFirstSector");
00246   const Short_t *&LSectG = TH("mGhostPairs.mLastSector");
00247 
00248   int ev = 0;
00249   int evGood = 0;
00250   while(TH.Next())// && ev ==0) 
00251   {    
00252     //checking vertex position
00253 
00254     vz->Fill(VZ);// 1 vertex per event, first fill vz and then apply the cuts
00255     
00256     vxy->Fill(VX,VY);//this is the 2d ok
00257 
00258     dvx->Fill(VX-VXmc);//delta VX
00259     dvy->Fill(VY-VYmc);
00260     dvz->Fill(VZ-VZmc);//delta VZ
00261     
00262     ev++;
00263     
00264     if(fabs(VZ)> VZCut)  continue;
00265     evGood++;
00266     
00267     // cout << Form("MiniMcTree:    (Good Event)/(Event number)  = %5d/%5d, VZ = %1.4f, Number of embedded tracks = %10d, Matched tracks = %10d", evGood, ev, VZ, ntrk1, ntrk) << endl;
00268     
00269     for (Int_t itr=0; itr<ntrk; itr++) 
00270       { 
00271 
00272         if ( Gid[itr]!=14 && Gid[itr]!=15) continue; //use this line  for Phi meson decaying in KK (14, 15)
00273         if (fabs(EtaM[itr]) > yCut) continue;
00274         
00275         //Defining Momentum
00276         
00277         float p = sqrt(PtM[itr]*PtM[itr]+PzM[itr]*PzM[itr]);
00278         
00279         if (NFit[itr] >= NFitCut  && NComm[itr] > NCommCut)
00280           {
00281            
00282             hDca->Fill(nch, PtM[itr],Dca[itr]);
00283             hPtEtaDcaM ->Fill( PtM[itr], EtaM[itr],  Dca[itr]);
00284             //of course x ->Pt, y->Eta, z->DCA
00285           }
00286         
00287         
00288         //if (Dca[itr] < maxDca && (1.*(NFit[itr])/NComm[itr])< 2.5 && NFit[itr]>=10)
00289         if (Dca[itr] < maxDca && NFit[itr]>=10)
00290           {
00291             hNfit ->Fill(nch, PtM[itr], NFit[itr]); 
00292             hPtEtaNfitM ->Fill( PtM[itr], EtaM[itr], NFit[itr]);
00293             //of course x ->Pt, y->Eta, z->NFIT
00294 
00295           }
00296         
00297         if( Dca[itr] < maxDca  && NFit[itr] >= NFitCut && NComm[itr] > NCommCut)
00298           dedx->Fill(p,dEdx[itr]*1e6);
00299         
00300         // filling histograms for global variables and energy loss
00301         
00302         if (Dca[itr]<0. || Dca[itr]>maxDca) continue;
00303         if(NFit[itr] >= NFitCut  && NComm[itr] > NCommCut )
00304           {
00305             hPtM->Fill(PtM[itr]);
00306             hEtaM->Fill(EtaM[itr]);
00307             hFSectM->Fill(FSectM[itr]);
00308             hLSectM->Fill(LSectM[itr]);
00309             Double_t mass0 = 0.13957;
00310             Double_t P0 = sqrt(PtM[itr]*PtM[itr]+PzM[itr]*PzM[itr]+mass0*mass0);
00311             Double_t P0mc = sqrt(PtMmc[itr]*PtMmc[itr]+PzMmc[itr]*PzMmc[itr]+mass0*mass0);
00312             Double_t rap = 0.5*log((P0+PzM[itr])/(P0-PzM[itr]));
00313             Double_t rapmc = 0.5*log((P0mc+PzMmc[itr])/(P0mc-PzMmc[itr]));
00314             hYM->Fill(rap);
00315             hPhiM->Fill(PhiM[itr]);
00316             hPtM_E->Fill(PtM[itr], PtM[itr]-PtMmc[itr]);
00317             hYM_E->Fill(rap, rap-rapmc);
00318             hPhiM_E->Fill(PhiM[itr], PhiM[itr]-PhiMmc[itr]);
00319             hPtYPhiM->Fill(PtM[itr],rap,PhiM[itr]);
00320           }
00321         
00322       }//closing for loop
00323     
00324     //filling Ghosts tracks
00325     
00326     for (Int_t itrG=0; itrG<ntrkG; itrG++)
00327       { 
00328         if (fabs(EtaG[itrG]) > yCut) continue;
00329         float pg = sqrt(PtG[itrG]*PtG[itrG]+PzG[itrG]*PzG[itrG]);
00330         
00331         if (NFitG[itrG] >= NFitCut)
00332           {
00333             hDcaG      ->Fill(nch, PtG[itrG],DcaG[itrG]);         
00334             hPtEtaDcaG ->Fill(PtG[itr], EtaG[itr], DcaG[itr]);
00335             //of course x-> Pt, y->Eta, z->DCA
00336           }
00337         
00338         if (DcaG[itrG] < maxDca && NFitG[itrG] >= 10)
00339           {
00340             hNfitG     ->Fill(nch, PtG[itrG], NFitG[itrG]); 
00341             hPtEtaNfitG->Fill( PtG[itr], EtaG[itr], NFitG[itr]);
00342             //of course x-> Pt, y->Eta, z->Nfit
00343           }
00344         
00345         if( DcaG[itrG] < maxDca && NFitG[itrG] >= NFitCut){
00346           dedxG->Fill(pg,dEdxG[itrG]*1e6);
00347           hPhiG->Fill(PhiG[itrG]);
00348           hFSectG->Fill(FSectG[itrG]);
00349           hLSectG->Fill(LSectG[itrG]);
00350         }
00351       }
00352     
00353     //filling Mc tracks(embedded)
00354     
00355     for (int itr1=0; itr1 < ntrk1; itr1++)
00356       {
00357         //if (fabs(EtaMc[itr1]) > yCut) continue;
00358         
00359         hPtMc->Fill(PtMc[itr1]);
00360         hEtaMc->Fill(EtaMc[itr1]);
00361         Double_t mass0 = 0.13957;
00362         Double_t P0 = sqrt(PtMc[itr1]*PtMc[itr1]+PzMc[itr1]*PzMc[itr1]+mass0*mass0);
00363         hYMc->Fill(0.5*log((P0+PzMc[itr1])/(P0-PzMc[itr1])));
00364         hPhiMc->Fill(PhiMc[itr1]);
00365         hPtYPhiMc->Fill(PtMc[itr1],0.5*log((P0+PzMc[itr1])/(P0-PzMc[itr1])),PhiMc[itr1]);
00366       }
00367     
00368   }//closing while loop
00369  
00370 
00371   //-------Writing-----------
00372 
00373   cout<< "Creating output file .... "<<oFile; flush(cout);  
00374 
00375   TFile *fout=new TFile(oFile,"recreate"); 
00376   fout->cd();
00377 
00378   // -- MC
00379   hMult->Write();
00380   nChTotal->Write();
00381   dedx->Write();
00382   dedxG->Write();
00383   hNfit->Write();
00384   hNfitG->Write();
00385   hDca-> Write();
00386 
00387   //Ghost
00388   hDcaG-> Write();
00389   hPhiG->Write();
00390   hFSectG->Write();
00391   hLSectG->Write();
00392   
00393   hPtEtaNfitG->Write();
00394   hPtEtaDcaG ->Write();
00395  
00396   //Global Variables
00397   v_xy->Write();
00398   vz->Write();
00399 
00400   dvx->Write();
00401   dvy->Write();
00402   dvz->Write();
00403 
00404   //---Matched
00405   hPtM->Write();
00406   hEtaM->Write();
00407   hFSectM->Write();
00408   hLSectM->Write();
00409   hYM->Write();
00410   hPhiM->Write();
00411   hPtM_E->Write();
00412   hYM_E->Write();
00413   hPhiM_E->Write();
00414   hPtYPhiM->Write();
00415   
00416   hPtEtaDcaM->Write();
00417   hPtEtaNfitM-> Write();
00418   
00419   //---Embedded
00420   hEtaMc->Write();
00421   hYMc->Write();
00422   hPtMc->Write();
00423   hPhiMc->Write();
00424   hPtYPhiMc->Write();
00425 
00426    fout->GetList()->Sort();
00427 
00428   fout->Close();
00429   cout<<" done"<<endl; flush(cout);
00430 }
macros/embedding/scan_embed_mc.C
