PWGDQ/dielectron/TakuAlberica/pair/ana_sgl.C

   1 #define ana_sgl_cxx
   2 #include "ana_sgl.h"
   3 #include <TH2.h>
   4 #include <TStyle.h>
   5 #include <TCanvas.h>
   6 #include <TMath.h>
   7 #include <vector>
   8 #include <deque>
   9 #include <cstdlib>
  10 #include <TRandom.h>
  11
  12 using namespace std;
  13
  14 using std::vector;
  15
  16 class etrk;
  17
  18 //Buffer for event mixing
  19 static const int NBUF=100; //depth of buffer
  20 static const int NMix=10; //# of events mixed (for +-)
  21 //static const int NMix=2; //# of events mixed (for +-)
  22
  23
  24 static const int NZBIN=10;
  25 static const int NCENT=10;
  26 int d_ibuf[NZBIN][NCENT];
  27 vector<etrk> d_vep[NBUF][NZBIN][NCENT];
  28 vector<etrk> d_vem[NBUF][NZBIN][NCENT];
  29
  30 static const unsigned int MAXPOOL=150;
  31 //static const unsigned int MAXPOOL=50;
  32 static const int MAX_TRY=3;
  33 deque<etrk> d_poolp[NZBIN][NCENT];
  34 deque<etrk> d_poolm[NZBIN][NCENT];
  35
  36 void reshuffle_buffer(vector<etrk> &ve,
  37                       deque<etrk> &pool){
  38   //If there is not enough electron in the pool, give up
  39   unsigned int ne = ve.size();
  40   unsigned int poolsize = pool.size();
  41   if(poolsize < ne) {
  42     cout <<" pool size="<<poolsize<<" ne"<<ne<<endl;
  43     return;
  44   }
  45   for(unsigned int ie=0; ie < ne; ie++) {
  46     int j = rand()%poolsize;
  47     ve[ie] = pool[j];
  48   }
  49 }
  50
  51 void ana_sgl::Loop()
  52 {
  53 //   In a ROOT session, you can do:
  54 //      Root > .L ana_sgl.C
  55 //      Root > ana_sgl t
  56 //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  57 //      Root > t.Show();       // Show values of entry 12
  58 //      Root > t.Show(16);     // Read and show values of entry 16
  59 //      Root > t.Loop();       // Loop on all entries
  60 //
  61
  62 //     This is the loop skeleton where:
  63 //    jentry is the global entry number in the chain
  64 //    ientry is the entry number in the current Tree
  65 //  Note that the argument to GetEntry must be:
  66 //    jentry for TChain::GetEntry
  67 //    ientry for TTree::GetEntry and TBranch::GetEntry
  68 //
  69 //       To read only selected branches, Insert statements like:
  70 // METHOD1:
  71 //    fChain->SetBranchStatus("*",0);  // disable all branches
  72 //    fChain->SetBranchStatus("branchname",1);  // activate branchname
  73 // METHOD2: replace line
  74 //    fChain->GetEntry(jentry);       //read all branches
  75 //by  b_branchname->GetEntry(ientry); //read only this branch
  76    if (fChain == 0) return;
  77
  78    Long64_t nentries = fChain->GetEntriesFast();
  79
  80    Long64_t nbytes = 0, nb = 0;
  81    for (Long64_t jentry=0; jentry<nentries;jentry++) {
  82       Long64_t ientry = LoadTree(jentry);
  83       if (ientry < 0) break;
  84       nb = fChain->GetEntry(jentry);   nbytes += nb;
  85       ana_event(ientry, jentry);
  86       // if (Cut(ientry) < 0) continue;
  87    }
  88 }
  89
  90 //___________________________________________
  91   void ana_sgl::ana_init(char *outname){
  92   fout  = new TFile(outname,"recreate");
  93   char name[100];
  94
  95   hdedx_pt = new TH2D("hdedx_pt","dedx vs. pT",100,0, 10, 100, 0, 200);
  96   hdedx_tof_elec_pt = new TH2D("hdedx_tof_elec_pt","dedx vs. pT with tof veto",100,0, 10, 100, 0, 200);
  97   hdedx_tof_all_pt = new TH2D("hdedx_tof_all_pt","dedx vs. pT with tof veto",100,0, 10, 100, 0, 200);
  98   hdedx_tof_elec_emc_pt = new TH2D("hdedx_tof_elec_emc_pt","dedx vs. pT with tof veto and EMC",100,0, 10, 100, 0, 200);
  99   hdedx_tof_all_emc_pt = new TH2D("hdedx_tof_all_emc_pt","dedx vs. pT with tof veto and EMC",100,0, 10, 100, 0, 200);
 100   hdedx_emc_pt = new TH2D("hdedx_emc_pt","dedx vs. pT with EMC",100,0, 10, 100, 0, 200);
 101
 102
 103   hbetatof_pt = new TH2D("hbetatof_pt","beta vs. pT",100,0, 10, 100, 0, 1);
 104   hbetatof_tof_elec_pt = new TH2D("hbetatof_tof_elec_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1);
 105   hbetatof_tof_all_pt = new TH2D("hbetatof_tof_all_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1);
 106   hbetatof_tof_elec_emc_pt = new TH2D("hbetatof_tof_elec_emc_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1);
 107   hbetatof_tof_all_emc_pt = new TH2D("hbetatof_tof_all_emc_pt","beta vs. pT with tof cut and emc E/p",100,0, 10, 100, 0, 1);
 108   hbetatof_emc_pt = new TH2D("hbetatof_emc_pt","dedx vs. pT with EMC",100,0, 10, 100, 0, 200);
 109
 110   hCentrality = new TH1F("hCentrality","hCentrality",100,0,100);
 111   hV0AC = new TH2F("hV0AC","V0AC",5000,0,15000,5000,0,15000);
 112   hV0AC_Ntrk = new TH2F("hV0AC_Ntrk","V0AC vs. Ntrk",5000,0,30000,5000,0,15000);
 113   hV0AC_NaccTrcklts  = new TH2F("hV0AC_NaccTrcklts","V0AC vs. Ntrklets",5000,0,30000,2500,0,5000);
 114
 115
 116
 117   fBinWidth = 0.05;
 118
 119   nHistos = (int)(10/fBinWidth);
 120
 121   for(int i=0;i<nHistos;i++){
 122     sprintf(name,"hdedx_p%d",i);
 123     hdedx[i] = new TH1D(name, name, 200,0,200);
 124
 125     sprintf(name,"hdedx_tof_elec_p%d",i);
 126     hdedx_tof_elec[i] = new TH1D(name, name, 200,0,200);
 127
 128     sprintf(name,"hdedx_tof_all_p%d",i);
 129     hdedx_tof_all[i] = new TH1D(name, name, 200,0,200);
 130
 131
 132     sprintf(name,"hdedx_tof_elec_emc_p%d",i);
 133     hdedx_tof_elec_emc[i] = new TH1D(name, name, 200,0,200);
 134
 135     sprintf(name,"hdedx_tof_all_emc_p%d",i);
 136     hdedx_tof_all_emc[i] = new TH1D(name, name, 200,0,200);
 137
 138     sprintf(name,"hdedx_emc_p%d",i);
 139     hdedx_emc[i] = new TH1D(name, name, 200,0,200);
 140
 141   }
 142
 143
 144   int nbinx=400;
 145   float max_x=20;
 146   float min_x=0.2;
 147   float binw = (TMath::Log(max_x)-TMath::Log(min_x))/nbinx;
 148   double xbin[401];
 149   for(int ii=0;ii<nbinx+1;ii++){
 150     xbin[ii] = TMath::Exp(TMath::Log(min_x) + 0.5*binw+binw*ii);
 151   }
 152
 153
 154   //// my histo
 155   fEventStat = new TH1D("EventStat ","EventStat ",7,0,7);
 156   fEvent = new TH1D("Event","Number of Events",   60,0,60);
 157   fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 158   fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt","dE/dX normalized to electron vs. pT of TPC",
 159                                  nbinx, xbin, 2000, -10, 10);
 160   fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2);
 161   fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10);
 162
 163   fdEdXvsPtWithCut = new TH2D("dEdXvsPtWithCuts","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200);
 164
 165
 166   fPtElec[0] = new TH1D("hPtElec0","elctron pt",250,0,5);
 167   fPtElec[1] = new TH1D("hPtElec1","elctron pt",250,0,5);
 168
 169
 170   fNelc_pos = new TH2D("hNelc_pos","# of electrons and positrons",40, -0.5, 39.5, 40, -0.5, 39.5);
 171   fNelc_all = new TH1D("hNelc_all","# of electrons and positrons",50, -0.5, 49.5);
 172   fNelc_all_pT =  new TH2D("hNelc_all_pT","# of electrons and positrons in each pT bin",10, 0, 5, 50, -0.5, 49.5);
 173
 174   for(int ic=0;ic<10;ic++){
 175     sprintf(name,"hNelc_pos_cent%d",ic);
 176     fNelc_pos_cent[ic] = new TH2D(name, name, 40, -0.5, 39.5, 40, -0.5, 39.5);
 177     sprintf(name,"hNelc_all_cent%d",ic);
 178     fNelc_all_cent[ic] = new TH1D(name, name, 50, -0.5, 49.5);
 179   }
 180
 181
 182   //// tree branch
 183   d_tree = new TTree("d_tree","single tree");
 184   d_tree->Branch("evt", &d_evt, "evt/I");
 185   d_tree->Branch("cent", &d_cent, "cent/D");
 186   d_tree->Branch("ntrk", &d_ntrk, "ntrk/D");
 187   d_tree->Branch("xvprim", &d_xvprim, "xvprim/D");
 188   d_tree->Branch("yvprim", &d_yvprim, "yvprim/D");
 189   d_tree->Branch("zvprim", &d_zvprim, "zvprim/D");
 190   d_tree->Branch("nacctrklets", &d_nacctrklets, "nacctrklets/D");
 191   d_tree->Branch("xres", &d_xres, "xres/D");
 192   d_tree->Branch("yres", &d_yres, "yres/D");
 193   d_tree->Branch("zres", &d_zres, "zres/D");
 194   d_tree->Branch("nelc", &d_nelc, "nelc/I");
 195   d_tree->Branch("px", d_px, "px[nelc]/D");
 196   d_tree->Branch("py", d_py, "py[nelc]/D");
 197   d_tree->Branch("pz", d_pz, "pz[nelc]/D");
 198   d_tree->Branch("p", d_p, "p[nelc]/D");
 199   d_tree->Branch("pt", d_pt, "pt[nelc]/D");
 200   d_tree->Branch("xv", d_xv, "xv[nelc]/D");
 201   d_tree->Branch("yv", d_yv, "yv[nelc]/D");
 202   d_tree->Branch("zv", d_zv, "zv[nelc]/D");
 203   d_tree->Branch("phi", d_phi, "phi[nelc]/D");
 204   d_tree->Branch("theta", d_theta, "theta[nelc]/D");
 205   d_tree->Branch("eta", d_eta, "eta[nelc]/D");
 206   d_tree->Branch("c", d_c, "c[nelc]/D");
 207   d_tree->Branch("nclusITS", d_nclusITS, "nclusITS[nelc]/D");
 208   d_tree->Branch("nclusTPC", d_nclusTPC, "nclusTPC[nelc]/D");
 209   d_tree->Branch("nclusTPCiter", d_nclusTPCiter, "nclusTPCiter[nelc]/D");
 210   d_tree->Branch("nfclusTPC", d_nfclusTPC, "nfclusTPC[nelc]/D");
 211   d_tree->Branch("nfclusTPCr", d_nfclusTPCr, "nfclusTPCr[nelc]/D");
 212   d_tree->Branch("nfclusTPCrFrac", d_nfclusTPCrFrac, "nfclusTPCrFrac[nelc]/D");
 213   d_tree->Branch("TPCsignalN", d_TPCsignalN, "TPCsignalN[nelc]/D");
 214   d_tree->Branch("TPCsignalNfrac", d_TPCsignalNfrac, "TPCsignalNfrac[nelc]/D");
 215   d_tree->Branch("TPCchi2cl", d_TPCchi2cl, "TPCchi2cl[nelc]/D");
 216   d_tree->Branch("trkstat", d_trkstat, "trkstat[nelc]/D");
 217   d_tree->Branch("nclsTRD", d_nclsTRD, "nclsTRD[nelc]/D");
 218   d_tree->Branch("TRDntracklets", d_TRDntracklets, "TRDntracklets[nelc]/D");
 219   d_tree->Branch("TRDpidquality", d_TRDpidquality, "TRDpidquality[nelc]/D");
 220   d_tree->Branch("TRDprobEle", d_TRDprobEle, "TRDprobEle[nelc]/D");
 221   d_tree->Branch("TRDprobPio", d_TRDprobPio, "TRDprobPio[nelc]/D");
 222   d_tree->Branch("impactXY", d_impactXY, "impactXY[nelc]/D");
 223   d_tree->Branch("impactZ", d_impactZ, "impactZ[nelc]/D");
 224   d_tree->Branch("tracklength", d_tracklength, "tracklength[nelc]/D");
 225   d_tree->Branch("ITSsignal", d_ITSsignal, "ITSsignal[nelc]/D");
 226   d_tree->Branch("ITSnsigmaEle", d_ITSnsigmaEle, "ITSnsigmaEle[nelc]/D");
 227   d_tree->Branch("ITSnsigmaPio", d_ITSnsigmaPio, "ITSnsigmaPio[nelc]/D");
 228   d_tree->Branch("ITSnsigmaMuo", d_ITSnsigmaMuo, "ITSnsigmaMuo[nelc]/D");
 229   d_tree->Branch("ITSnsigmaKao", d_ITSnsigmaKao, "ITSnsigmaKao[nelc]/D");
 230   d_tree->Branch("ITSnsigmaPro", d_ITSnsigmaPro, "ITSnsigmaPro[nelc]/D");
 231   d_tree->Branch("PIn", d_PIn, "PIn[nelc]/D");
 232   d_tree->Branch("TPCsignal", d_TPCsignal, "TPCsignal[nelc]/D");
 233   d_tree->Branch("TOFsignal", d_TOFsignal, "TOFsignal[nelc]/D");
 234   d_tree->Branch("TOFbeta", d_TOFbeta, "TOFbeta[nelc]/D");
 235   d_tree->Branch("TPCnSigmaEle", d_TPCnSigmaEle, "TPCnSigmaEle[nelc]/D");
 236   d_tree->Branch("TPCnSigmaPio", d_TPCnSigmaPio, "TPCnSigmaPio[nelc]/D");
 237   d_tree->Branch("TPCnSigmaMuo", d_TPCnSigmaMuo, "TPCnSigmaMuo[nelc]/D");
 238   d_tree->Branch("TPCnSigmaKao", d_TPCnSigmaKao, "TPCnSigmaKao[nelc]/D");
 239   d_tree->Branch("TPCnSigmaPro", d_TPCnSigmaPro, "TPCnSigmaPro[nelc]/D");
 240   d_tree->Branch("TOFnSigmaEle", d_TOFnSigmaEle, "TOFnSigmaEle[nelc]/D");
 241   d_tree->Branch("TOFnSigmaPio", d_TOFnSigmaPio, "TOFnSigmaPio[nelc]/D");
 242   d_tree->Branch("TOFnSigmaMuo", d_TOFnSigmaMuo, "TOFnSigmaMuo[nelc]/D");
 243   d_tree->Branch("TOFnSigmaKao", d_TOFnSigmaKao, "TOFnSigmaKao[nelc]/D");
 244   d_tree->Branch("TOFnSigmaPro", d_TOFnSigmaPro, "TOFnSigmaPro[nelc]/D");
 245   d_tree->Branch("E", d_E, "E[nelc]/D");
 246   d_tree->Branch("dphi", d_phi, "dphi[nelc]/D");
 247   d_tree->Branch("deta", d_eta, "deta[nelc]/D");
 248
 249   d_tree->Branch("chi2ndf", d_chi2ndf, "chi2ndf[nelc]/D");
 250
 251   cout<<"aaaaaaaaaaa"<<endl;
 252   d_ntpair = new TTree("ntpair","pair");
 253   d_ntpair->Branch("run", &d_run, "run/D");
 254   d_ntpair->Branch("event", &d_event, "event/I");
 255   d_ntpair->Branch("centrality",&d_centrality,"centrality/D");
 256   d_ntpair->Branch("prim_xv",&d_prim_xv,"prim_xv/D");
 257   d_ntpair->Branch("prim_yv",&d_prim_yv,"prim_yv/D");
 258   d_ntpair->Branch("prim_zv",&d_prim_zv,"prim_zv/D");
 259   d_ntpair->Branch("mass",&d_mass,"mass/D");
 260   d_ntpair->Branch("pxpair",&d_pxpair,"pxpair/D");
 261   d_ntpair->Branch("pypair",&d_pypair,"pypair/D");
 262   d_ntpair->Branch("pzpair",&d_pzpair,"pzpair/D");
 263   d_ntpair->Branch("ptpair",&d_ptpair,"ptpair/D");
 264   d_ntpair->Branch("epair",&d_epair,"epair/D");
 265   d_ntpair->Branch("eta",&d_etapair,"eta/D");
 266   d_ntpair->Branch("phi",&d_phipair,"phi/D");
 267   d_ntpair->Branch("cos",&d_cos,"cos/D");
 268   d_ntpair->Branch("phiv",&d_phiv,"phiv/D");
 269   d_ntpair->Branch("psi",&d_psi,"psi/D");
 270   d_ntpair->Branch("pairtype",&d_pairtype,"pairtype/I");
 271   d_ntpair->Branch("cent1",&d_cent1,"cent1/D");
 272   d_ntpair->Branch("xv1",&d_xv1,"xv1/D");
 273   d_ntpair->Branch("yv1",&d_yv1,"yv1/D");
 274   d_ntpair->Branch("zv1",&d_zv1,"zv1/D");
 275   d_ntpair->Branch("px1",&d_px1,"px1/D");
 276   d_ntpair->Branch("py1",&d_py1,"py1/D");
 277   d_ntpair->Branch("pz1",&d_pz1,"pz1/D");
 278   d_ntpair->Branch("pt1",&d_pt1,"pt1/D");
 279   d_ntpair->Branch("eta1",&d_eta1,"eta1/D");
 280   d_ntpair->Branch("phi1",&d_phi1,"phi1/D");
 281   d_ntpair->Branch("theta1",&d_theta1,"theta1/D");
 282   d_ntpair->Branch("tpc1",&d_tpc1,"tpc1/D");
 283   d_ntpair->Branch("ntpc_ele1",&d_ntpc_ele1,"ntpc_ele1/D");
 284   d_ntpair->Branch("ntpc_pio1",&d_ntpc_pio1,"ntpc_pio1/D");
 285   d_ntpair->Branch("ntpc_kao1",&d_ntpc_kao1,"ntpc_kao1/D");
 286   d_ntpair->Branch("ntpc_pro1",&d_ntpc_pro1,"ntpc_pro1/D");
 287   d_ntpair->Branch("beta1",&d_beta1,"beta1/D");
 288   d_ntpair->Branch("ntof_ele1",&d_ntof_ele1,"ntof_ele1/D");
 289   d_ntpair->Branch("ntof_pio1",&d_ntof_pio1,"ntof_pio1/D");
 290   d_ntpair->Branch("ntof_kao1",&d_ntof_kao1,"ntof_kao1/D");
 291   d_ntpair->Branch("ntof_pro1",&d_ntof_pro1,"ntof_pro1/D");
 292   d_ntpair->Branch("its1",&d_its1,"its1/D");
 293   d_ntpair->Branch("nits1",&d_nits1,"nits1/D");
 294   d_ntpair->Branch("ntpc1",&d_ntpc1,"ntpc1/D");
 295   d_ntpair->Branch("e1",&d_e1,"e1/D");
 296   d_ntpair->Branch("dphi1",&d_dphi1,"dphi1/D");
 297   d_ntpair->Branch("deta1",&d_deta1,"deta1/D");
 298   d_ntpair->Branch("dcaxy1",&d_dcaxy1,"dcaxy1/D");
 299   d_ntpair->Branch("dcaz1",&d_dcaz1,"dcaz1/D");
 300   d_ntpair->Branch("conv1",&d_conv1,"conv1/I");
 301
 302
 303
 304
 305   d_ntpair->Branch("cent2",&d_cent2,"cent2/D");
 306   d_ntpair->Branch("xv2",&d_xv2,"xv2/D");
 307   d_ntpair->Branch("yv2",&d_yv2,"yv2/D");
 308   d_ntpair->Branch("zv2",&d_zv2,"zv2/D");
 309   d_ntpair->Branch("px2",&d_px2,"px2/D");
 310   d_ntpair->Branch("py2",&d_py2,"py2/D");
 311   d_ntpair->Branch("pz2",&d_pz2,"pz2/D");
 312   d_ntpair->Branch("pt2",&d_pt2,"pt2/D");
 313   d_ntpair->Branch("eta2",&d_eta2,"eta2/D");
 314   d_ntpair->Branch("phi2",&d_phi2,"phi2/D");
 315   d_ntpair->Branch("theta2",&d_theta2,"theta2/D");
 316   d_ntpair->Branch("tpc2",&d_tpc2,"tpc2/D");
 317   d_ntpair->Branch("ntpc_ele2",&d_ntpc_ele2,"ntpc_ele2/D");
 318   d_ntpair->Branch("ntpc_pio2",&d_ntpc_pio2,"ntpc_pio2/D");
 319   d_ntpair->Branch("ntpc_kao2",&d_ntpc_kao2,"ntpc_kao2/D");
 320   d_ntpair->Branch("ntpc_pro2",&d_ntpc_pro2,"ntpc_pro2/D");
 321   d_ntpair->Branch("beta2",&d_beta2,"beta2/D");
 322   d_ntpair->Branch("ntof_ele2",&d_ntof_ele2,"ntof_ele2/D");
 323   d_ntpair->Branch("ntof_pio2",&d_ntof_pio2,"ntof_pio2/D");
 324   d_ntpair->Branch("ntof_kao2",&d_ntof_kao2,"ntof_kao2/D");
 325   d_ntpair->Branch("ntof_pro2",&d_ntof_pro2,"ntof_pro2/D");
 326   d_ntpair->Branch("its2",&d_its2,"its2/D");
 327   d_ntpair->Branch("nits2",&d_nits2,"nits2/D");
 328   d_ntpair->Branch("ntpc2",&d_ntpc2,"ntpc2/D");
 329   d_ntpair->Branch("e2",&d_e2,"e2/D");
 330   d_ntpair->Branch("dphi2",&d_dphi2,"dphi2/D");
 331   d_ntpair->Branch("deta2",&d_deta2,"deta2/D");
 332   d_ntpair->Branch("dcaxy2",&d_dcaxy2,"dcaxy2/D");
 333   d_ntpair->Branch("dcaz2",&d_dcaz2,"dcaz2/D");
 334   d_ntpair->Branch("conv2",&d_conv2,"conv2/I");
 335
 336   cout<<"aaaaaaaaaaa"<<endl;
 337   for(int i=0;i<7;i++){
 338     for(int j=0;j<11;j++){
 339       sprintf(name,"hmasspt_cent%d_pair%d",j,i);
 340       hmasspt[i][j] = new TH2D(name, name, 500, 0, 5, 500, 0, 5);
 341       sprintf(name,"hmasspt_weight_cent%d_pair%d",j,i);
 342       hmasspt_weight[i][j] = new TH2D(name, name, 500, 0, 5, 500, 0, 5);
 343     }
 344   }
 345
 346   vep.clear();
 347   vem.clear();
 348   vep_tmp.clear();
 349   vem_tmp.clear();
 350
 351   d_evt = 0;
 352   d_event = 0;
 353
 354   simflag = false;
 355   d_conv_flag = false;
 356
 357   for(int i=0;i<10;i++){
 358     nelec_pos[i] = 0;
 359   }
 360
 361   magnetic_field_mm = true;
 362
 363   d_flag_tof_cut = false;
 364   d_flag_emc_cut =false;
 365   d_flag_phiv = false;
 366   d_tpc_dedx_low = 75;
 367   d_tpc_dedx_high = 90;
 368   d_tof_low = -3;
 369   d_tof_high = 3;
 370   d_emc_low = 0.7;
 371   d_emc_high = 1.3;
 372   d_phiv_cut = 0.6;
 373
 374   d_flag_kaon_veto = false;
 375   d_flag_proton_veto = false;
 376   d_dedx_kaon_veto_low = -2;
 377   d_dedx_kaon_veto_high = 2;
 378   d_dedx_proton_veto_low = -2;
 379   d_dedx_proton_veto_high = 2;
 380
 381
 382   cout<<"ana_init end:"<<endl;
 383 }
 384
 385 //___________________________________________
 386 void ana_sgl::ana_end(void){
 387
 388   fout->cd();
 389   /*
 390   hdedx_pt->Write();
 391   hdedx_tof_elec_pt->Write();
 392   hdedx_tof_all_pt->Write();
 393   hdedx_tof_elec_emc_pt->Write();
 394   hdedx_tof_all_emc_pt->Write();
 395   hdedx_emc_pt->Write();
 396
 397   hbetatof_pt->Write();
 398   hbetatof_tof_elec_pt->Write();
 399   hbetatof_tof_all_pt->Write();
 400   hbetatof_tof_elec_emc_pt->Write();
 401   hbetatof_tof_all_emc_pt->Write();
 402   hbetatof_emc_pt->Write();
 403
 404   for(int i=0;i<nHistos;i++){
 405     hdedx[i]->Write();
 406     hdedx_tof_elec[i]->Write();
 407     hdedx_tof_all[i]->Write();
 408     hdedx_tof_elec_emc[i]->Write();
 409     hdedx_tof_all_emc[i]->Write();
 410     hdedx_emc[i]->Write();
 411   }
 412
 413   */
 414
 415   if(simflag==false){
 416     fEventStat->Write();
 417     fEvent->Write();
 418     fdEdXvsPt->Write();
 419     fdEdXnSigmaElecvsPt->Write();
 420     fTOFbetavsPt->Write();
 421     fTOFnSigmaElecvsPt->Write();
 422     fdEdXvsPtWithCut->Write();
 423     fPtElec[0]->Write();
 424     fPtElec[1]->Write();
 425   }
 426   fNelc_pos->Write();
 427   fNelc_all->Write();
 428   for(int i=0;i<10;i++){
 429     fNelc_pos_cent[i]->Write();
 430     fNelc_all_cent[i]->Write();
 431   }
 432   fNelc_all_pT->Write();
 433   //d_tree->Write();
 434   hCentrality->Write();
 435   hV0AC->Write();
 436   hV0AC_Ntrk->Write();
 437   hV0AC_NaccTrcklts->Write();
 438   for(int j=0;j<11;j++){
 439     for(int i=0;i<7;i++){
 440       hmasspt[i][j]->Write();
 441       hmasspt_weight[i][j]->Write();
 442     }
 443   }
 444
 445   //d_ntpair->Write();
 446   fout->Close();
 447 }
 448
 449 //____________________________________________
 450 void ana_sgl::loop_a_file(char *file){
 451
 452   TFile *treefile = TFile::Open(file);
 453   TDirectory *d = (TDirectory*)treefile->Get("PWG3_dielectron");
 454   if(d==0){
 455     cout<<" PWG3_dielectron is not found "<<endl;
 456     return ;
 457   }
 458
 459   //TTree *tree = (TTree*)d->Get("t");
 460   TTree *tree = (TTree*)d->Get("tree_MultiDie_CENT1");
 461   if(tree == 0) {
 462     cout <<"tree is not found in "<<file<<endl;
 463     treefile->Close();
 464     return;
 465   }
 466   cout << file <<" is opened"<<endl;
 467   if(simflag==false){
 468     //add_histograms(treefile);
 469   }
 470   Init(tree);
 471   Loop();
 472   tree->Clear();
 473   d->Clear();
 474   tree->Delete();
 475   d->Delete();
 476   delete d;
 477   treefile->Close();
 478   delete treefile;
 479
 480
 481
 482   cout <<"one file processed"<<endl;
 483 }
 484
 485 //____________________________________________
 486 void ana_sgl::ana_event(int ientry, int jentry){
 487
 488   //select trigger class:
 489   if(sel_trigger==2){
 490     if(kTriggerCent<100){
 491       return ;
 492     }
 493   }else if(sel_trigger==1){
 494     if(kTriggerCent%100<10){
 495       return ;
 496     }
 497   }else if(sel_trigger==0){
 498     if(kTriggerCent%10!=1){
 499       return ;
 500     }
 501   }
 502
 503   //  if(fkRunNumber>138350){
 504   /*
 505   if(kMag>0){
 506     magnetic_field_mm = false;
 507   }else{
 508     magnetic_field_mm = true;
 509   }
 510   */
 511   if(fkRunNumber>169591){
 512     magnetic_field_mm = false;
 513   }else{
 514     magnetic_field_mm = true;
 515   }
 516
 517   if(ientry%1000==0){
 518     cout<<" event processing "<<ientry<<" / "<<d_evt<<" / "<<d_event<<" : trigger "<<kTriggerCent<<endl;
 519   }
 520
 521   d_nelc = 0;
 522   for(int i=0;i<10;i++){
 523     nelec_pos[i] = 0;
 524   }
 525
 526   hCentrality->Fill(fkCentrality);
 527   hV0AC->Fill(fkV0C, fkV0A);
 528   hV0AC_Ntrk->Fill(fkNTrk, fkV0C+fkV0A);
 529   hV0AC_NaccTrcklts->Fill(fkNaccTrcklts, fkV0C+fkV0A);
 530   fill_to_tree_variables();
 531
 532   for(int i=0;i<fkNPar;i++){
 533
 534     if(simflag==false){
 535
 536       //fdEdXvsPtWithCut->Fill(kPIn[i], kTPCsignal[i]);
 537
 538       if(GlobalTrackcut(i)==false){
 539         continue;
 540       }
 541       fill_histograms(i);
 542     }
 543
 544
 545     fill_to_tree_track_variables(i);
 546
 547     //if(PairTrackcut(i)==false){
 548     //continue;
 549     //}
 550
 551     if(fkCentrality>0 && fkCentrality<10){
 552       int iptbin = (int)(kP[i]/0.5);
 553       if(iptbin>=10){
 554         iptbin=9;
 555       }
 556       nelec_pos[iptbin]++;
 557     }
 558
 559     if(kCharge[i]>0){
 560       etrk e = etrk(
 561                     fkCentrality, fkXvPrim, fkYvPrim, fkZvPrim,
 562                     kXv[i], kYv[i], kZv[i],
 563                     kPx[i], kPy[i], kPz[i], kPt[i],
 564                     kEta[i], kPhi[i], kTheta[i],
 565                     kTPCsignal[i], kTOFbeta[i],
 566                     kE[i], kDeltaPhi[i], kDeltaEta[i],
 567                     kTPCnSigmaEle[i], kTPCnSigmaPio[i], kTPCnSigmaKao[i], kTPCnSigmaPro[i],
 568                     kTOFnSigmaEle[i], kTOFnSigmaPio[i], kTOFnSigmaKao[i], kTOFnSigmaPro[i],
 569                     kITSsignal[i], kNclsITS[i], kNclsTPC[i], kLegDistXY[i], kLegDist[i]
 570                   );
 571       vep_tmp.push_back(e);
 572     }else{
 573       etrk e = etrk(
 574                     fkCentrality, fkXvPrim, fkYvPrim, fkZvPrim,
 575                     kXv[i], kYv[i], kZv[i],
 576                     kPx[i], kPy[i], kPz[i], kPt[i],
 577                     kEta[i], kPhi[i], kTheta[i],
 578                     kTPCsignal[i], kTOFbeta[i],
 579                     kE[i], kDeltaPhi[i], kDeltaEta[i],
 580                     kTPCnSigmaEle[i], kTPCnSigmaPio[i], kTPCnSigmaKao[i], kTPCnSigmaPro[i],
 581                     kTOFnSigmaEle[i], kTOFnSigmaPio[i], kTOFnSigmaKao[i], kTOFnSigmaPro[i],
 582                     kITSsignal[i], kNclsITS[i], kNclsTPC[i], kLegDistXY[i], kLegDist[i]
 583                     );
 584       vem_tmp.push_back(e);
 585     }
 586   }
 587
 588   //////// fill to the tree //////////////////
 589   //d_tree->Fill();
 590
 591
 592
 593   if(d_conv_flag==true){
 594     check_conversion_pairs(vep_tmp, vem_tmp);
 595   }
 596
 597   check_ghost_pairs(vep_tmp);
 598   check_ghost_pairs(vem_tmp);
 599   randomize_pool(vep_tmp, vem_tmp);
 600   if(d_conv_flag==false){
 601     check_conversion_pairs(vep, vem);
 602   }
 603
 604   fNelc_pos->Fill(vep.size(), vem.size());
 605   fNelc_all->Fill(vep.size()+vem.size());
 606
 607   if(fkCentrality>0 && fkCentrality<10){
 608     for(int i=0;i<10;i++){
 609       fNelc_all_pT->Fill(0.25+0.5*i, nelec_pos[i]);
 610     }
 611   }
 612
 613
 614   calc_pair(vep, vem);
 615
 616   int icent = (int)(fkCentrality/10.0);
 617   int izbin = (int)((fkZvPrim+10)/2.0);
 618   if(icent<0) icent=0;
 619   if(icent>=NCENT) icent=NCENT-1;
 620   if(izbin<0) izbin=0;
 621   if(izbin>=NZBIN) izbin=NZBIN-1;
 622
 623   fNelc_pos_cent[icent]->Fill(vep.size(), vem.size());
 624   fNelc_all_cent[icent]->Fill(vep.size()+vem.size());
 625
 626
 627
 628   d_vep[d_ibuf[izbin][icent]][izbin][icent].clear();
 629   vector<etrk>::iterator iep;
 630   for(iep = vep.begin();iep != vep.end();++iep) {
 631     d_vep[d_ibuf[izbin][icent]][izbin][icent].push_back(*iep);
 632     d_poolp[izbin][icent].push_back(*iep);
 633     if(d_poolp[izbin][icent].size()>MAXPOOL) {
 634       d_poolp[izbin][icent].pop_front();
 635     }
 636   }
 637   d_vem[d_ibuf[izbin][icent]][izbin][icent].clear();
 638   vector<etrk>::iterator iem;
 639   for(iem = vem.begin();iem != vem.end();++iem) {
 640     d_vem[d_ibuf[izbin][icent]][izbin][icent].push_back(*iem);
 641     d_poolm[izbin][icent].push_back(*iem);
 642     if(d_poolm[izbin][icent].size()>MAXPOOL) {
 643       d_poolm[izbin][icent].pop_front();
 644     }
 645   }
 646   // Update the buffer pointer
 647   d_ibuf[izbin][icent]++;
 648   if(d_ibuf[izbin][icent]>= NBUF) d_ibuf[izbin][icent]=0;
 649
 650   ///////////////////////////////////////////
 651   vem.clear();
 652   vep.clear();
 653   vem_tmp.clear();
 654   vep_tmp.clear();
 655   d_evt ++;
 656   d_event ++;
 657
 658 }
 659
 660
 661 //____________________________________________
 662 bool ana_sgl::kTOFcut(int itrk){
 663
 664   if(kTOFsignal[itrk]!=9999 && kTOFbeta[itrk]>0.3 &&
 665      (kTOFnSigmaEle[itrk]<3 && kTOFnSigmaEle[itrk]>-3)
 666      /*
 667      (kTOFnSigmaKao>3 || kTOFnSigmaKao<-3) &&
 668      (kTOFnSigmaPro>3 || kTOFnSigmaPro<-3)
 669      */
 670
 671      ){
 672
 673     return true;
 674   }else{
 675     return false;
 676   }
 677 }
 678
 679 //____________________________________________
 680 bool ana_sgl::GlobalTrackcut(int itrk){
 681
 682
 683   if((
 684      (d_flag_kaon_veto == true && (kTPCnSigmaKao[itrk]<d_dedx_kaon_veto_low || kTPCnSigmaKao[itrk]>d_dedx_kaon_veto_high)) ||
 685      (d_flag_kaon_veto == false)
 686      )
 687      && (
 688          (d_flag_proton_veto == true && (kTPCnSigmaPro[itrk]<d_dedx_proton_veto_low || kTPCnSigmaPro[itrk]>d_dedx_proton_veto_high)) ||
 689          (d_flag_proton_veto == false)
 690          )
 691      ){
 692     fdEdXvsPtWithCut->Fill(kPIn[itrk], kTPCsignal[itrk]);
 693   }
 694
 695   if( kNclsTPC[itrk]>120
 696       && kTPCsignal[itrk]>d_tpc_dedx_low && kTPCsignal[itrk]<d_tpc_dedx_high
 697       && (
 698           (d_flag_tof_cut==true && kTOFnSigmaEle[itrk]>d_tof_low && kTOFnSigmaEle[itrk]<d_tof_high) ||
 699           (d_flag_tof_cut==false)
 700           )
 701       && (
 702           (d_flag_pt_cut == true && kPt[itrk]>d_pt_cut_low && kPt[itrk]<d_pt_cut_high) ||
 703           (d_flag_pt_cut == false)
 704           )
 705       && (
 706           (d_flag_kaon_veto == true && (kTPCnSigmaKao[itrk]<d_dedx_kaon_veto_low || kTPCnSigmaKao[itrk]>d_dedx_kaon_veto_high)) ||
 707           (d_flag_kaon_veto == false)
 708           )
 709       && (
 710           (d_flag_proton_veto == true && (kTPCnSigmaPro[itrk]<d_dedx_proton_veto_low || kTPCnSigmaPro[itrk]>d_dedx_proton_veto_high)) ||
 711           (d_flag_proton_veto == false)
 712           )
 713       ){
 714     return true;
 715   }else{
 716     return false;
 717   }
 718 }
 719
 720 /*
 721 //____________________________________________
 722 bool ana_sgl::PairTrackcut(int itrk){
 723
 724   if( (kTOFnSigmaEle[itrk]<3 && kTOFnSigmaEle[itrk]>-3) &&
 725       (kTOFnSigmaKao[itrk]>3||kTOFnSigmaKao[itrk]<-3) &&
 726       (kTOFnSigmaPro[itrk]>3||kTOFnSigmaPro[itrk]<-3) //&&
 727       //(kTPCnSigmaEle[itrk]-1.65226*exp(-kPIn[itrk]*kPIn[itrk]*1.60890)+0.838)<1 &&
 728       //(kTPCnSigmaEle[itrk]-1.65226*exp(-kPIn[itrk]*kPIn[itrk]*1.60890)+0.838)>-1
 729       ){
 730     return true;
 731   }else{
 732     return false;
 733   }
 734
 735
 736 //  if(
 737 //     !(kTPCsignal>75-210*(kPIn-0.5) &&
 738 //       kTPCsignal<100-190*(kPIn-0.5)) &&
 739 //     !(kTPCsignal>130-120*(kPIn-0.5) &&
 740 //       kTPCsignal<135-80*(kPIn-0.5)) &&
 741 //     (kTPCnSigmaEle-1.65226*exp(-kPIn*kPIn*1.60890)+0.838)<1 &&
 742 //      (kTPCnSigmaEle-1.65226*exp(-kPIn*kPIn*1.60890)+0.838)>-1
 743 //     ){
 744 //    return true;
 745 //  }else{
 746 //   return false;
 747 //  }
 748
 749
 750 }
 751 */
 752
 753 //____________________________________________
 754 void ana_sgl::fill_histograms(int itrk){
 755   /*
 756   hdedx_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 757   hbetatof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]);
 758
 759   if(kTOFcut(itrk)==true){
 760     hdedx_tof_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 761     hbetatof_tof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]);
 762   }
 763
 764   int iptbin = (int)((kPIn[itrk])/fBinWidth);
 765   if(iptbin>=0 && iptbin<nHistos){
 766     hdedx[iptbin]->Fill(kTPCsignal[itrk]);
 767     if(kTOFcut(itrk)==true)hdedx_tof[iptbin]->Fill(kTPCsignal[itrk]);
 768   }
 769   */
 770
 771   if(kCharge[itrk]>0){
 772     fPtElec[1]->Fill(kPIn[itrk]);
 773   }else{
 774     fPtElec[0]->Fill(kPIn[itrk]);
 775   }
 776
 777   hdedx_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 778   hbetatof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]);
 779
 780   if(kTOFnSigmaEle[itrk]>-3 && kTOFnSigmaEle[itrk]<3){
 781     hdedx_tof_elec_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 782     hbetatof_tof_elec_pt->Fill(kPIn[itrk], kTOFbeta[itrk]);
 783
 784     if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 785       hdedx_tof_elec_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 786       hbetatof_emc_pt->Fill(kPIn[itrk], kTOFbeta[itrk]);
 787     }
 788
 789     if( (kTOFnSigmaKao[itrk]>3 || kTOFnSigmaKao[itrk]<-3) &&
 790         (kTOFnSigmaPro[itrk]>3 || kTOFnSigmaPro[itrk]<-3)
 791         ){
 792       hdedx_tof_all_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 793       hbetatof_tof_all_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 794       if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 795         hdedx_tof_all_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 796         hbetatof_tof_all_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 797       }
 798     }
 799   }
 800
 801   if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 802     hdedx_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 803     hbetatof_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]);
 804   }
 805
 806
 807
 808
 809   int iptbin = (int)((kPIn[itrk])/fBinWidth);
 810   if(iptbin>=0 && iptbin<nHistos){
 811     hdedx[iptbin]->Fill(kTPCsignal[itrk]);
 812     //if(kTOFcut(itrk)==true)hdedx_tof[iptbin]->Fill(kTPCsignal[itrk]);
 813
 814     if(kTOFnSigmaEle[itrk]>-3 && kTOFnSigmaEle[itrk]<3){
 815       hdedx_tof_elec[iptbin]->Fill(kTPCsignal[itrk]);
 816
 817       if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 818         hdedx_tof_elec_emc[iptbin]->Fill(kTPCsignal[itrk]);
 819       }
 820       if( (kTOFnSigmaKao[itrk]>3 || kTOFnSigmaKao[itrk]<-3) &&
 821           (kTOFnSigmaPro[itrk]>3 || kTOFnSigmaPro[itrk]<-3)
 822           ){
 823         hdedx_tof_all[iptbin]->Fill(kTPCsignal[itrk]);
 824         if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 825           hdedx_tof_all_emc[iptbin]->Fill(kTPCsignal[itrk]);
 826         }
 827       }
 828     }
 829     if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){
 830       hdedx_emc[iptbin]->Fill(kTPCsignal[itrk]);
 831     }
 832   }
 833
 834 }
 835
 836 //____________________________________________
 837 void ana_sgl::add_histograms(TFile *fin){
 838
 839   TDirectory *d = (TDirectory*)fin->Get("PWG3_dielectron");
 840   if(d==0){
 841     cout<<" PWG3_dielectron is not found "<<endl;
 842     return ;
 843   }
 844   fEventStat->Add((TH1D*)d->Get("hEventStat_MultiDie_CENT1"));
 845   TList *list = (TList*)d->Get("jpsi_QA_CENT1");
 846   TList *listQA = (TList*)list->FindObject("QAElectron");
 847   fEvent->Add((TH1D*)listQA->FindObject("Event"));
 848   fdEdXvsPt->Add((TH1D*)listQA->FindObject("dEdXvsPt"));
 849   fdEdXnSigmaElecvsPt->Add((TH1D*)listQA->FindObject("fdEdXnSigmaElecvsPt"));
 850   fTOFbetavsPt->Add((TH1D*)listQA->FindObject("fTOFbetavsPt"));
 851   fTOFnSigmaElecvsPt->Add((TH1D*)listQA->FindObject("fTOFnSigmaElecvsPt"));
 852
 853   cout<<" ana_sgl::add_histograms "<<fin->GetName()<<" done "<<endl;
 854
 855 }
 856
 857 //_____________________________________________
 858 void ana_sgl::fill_to_tree_variables(void){
 859
 860   d_cent =  fkCentrality;
 861   d_ntrk= fkNTrk;
 862   d_xvprim= fkXvPrim;
 863   d_yvprim= fkYvPrim;
 864   d_zvprim= fkZvPrim;
 865   d_nacctrklets= fkNaccTrcklts;
 866   d_xres= fkXRes;
 867   d_yres= fkYRes;
 868   d_zres= fkZRes;
 869
 870 }
 871 //_____________________________________________
 872 void ana_sgl::fill_to_tree_track_variables(int itrk){
 873
 874
 875   d_px[d_nelc]= kPx[itrk];
 876   d_py[d_nelc]= kPy[itrk];
 877   d_pz[d_nelc]= kPz[itrk];
 878   d_p[d_nelc]= kP[itrk];
 879   d_pt[d_nelc]= kPt[itrk];
 880   d_xv[d_nelc]= kXv[itrk];
 881   d_yv[d_nelc]= kYv[itrk];
 882   d_zv[d_nelc]= kZv[itrk];
 883   d_phi[d_nelc]= kPhi[itrk];
 884   d_theta[d_nelc]= kTheta[itrk];
 885   d_eta[d_nelc]= kEta[itrk];
 886   d_c[d_nelc]= kCharge[itrk];
 887   d_nclusITS[d_nelc]= kNclsITS[itrk];
 888   d_nclusTPC[d_nelc]= kNclsTPC[itrk];
 889   d_nclusTPCiter[d_nelc]= kNclsTPCiter1[itrk];
 890   d_nfclusTPC[d_nelc]= kNFclsTPC[itrk];
 891   d_nfclusTPCr[d_nelc]= kNFclsTPCr[itrk];
 892   d_nfclusTPCrFrac[d_nelc]= kNFclsTPCrFrac[itrk];
 893   d_TPCsignalN[d_nelc]= kTPCsignalN[itrk];
 894   d_TPCsignalNfrac[d_nelc]= kTPCsignalNfrac[itrk];
 895   d_TPCchi2cl[d_nelc]= kTPCchi2Cl[itrk];
 896   d_trkstat[d_nelc]= kTrackStatus[itrk];
 897   d_nclsTRD[d_nelc]= kNclsTRD[itrk];
 898   d_TRDntracklets[d_nelc]= kTRDntracklets[itrk];
 899   d_TRDpidquality[d_nelc]= kTRDpidQuality[itrk];
 900   d_TRDprobEle[d_nelc]= kTRDprobEle[itrk];
 901   d_TRDprobPio[d_nelc]= kTRDprobPio[itrk];
 902   d_impactXY[d_nelc]= kImpactParXY[itrk];
 903   d_impactZ[d_nelc]= kImpactParZ[itrk];
 904   d_tracklength[d_nelc]= kTrackLength[itrk];
 905   d_ITSsignal[d_nelc]= kITSsignal[itrk];
 906   d_ITSnsigmaEle[d_nelc]= kITSnSigmaEle[itrk];
 907   d_ITSnsigmaPio[d_nelc]= kITSnSigmaPio[itrk];
 908   d_ITSnsigmaMuo[d_nelc]= kITSnSigmaMuo[itrk];
 909   d_ITSnsigmaKao[d_nelc]= kITSnSigmaKao[itrk];
 910   d_ITSnsigmaPro[d_nelc]= kITSnSigmaPro[itrk];
 911   d_PIn[d_nelc]= kPIn[itrk];
 912   d_TPCsignal[d_nelc]= kTPCsignal[itrk];
 913   d_TOFsignal[d_nelc]= kTOFsignal[itrk];
 914   d_TOFbeta[d_nelc]= kTOFbeta[itrk];
 915   d_TPCnSigmaEle[d_nelc]= kTPCnSigmaEle[itrk];
 916   d_TPCnSigmaPio[d_nelc]= kTPCnSigmaPio[itrk];
 917   d_TPCnSigmaMuo[d_nelc]= kTPCnSigmaMuo[itrk];
 918   d_TPCnSigmaKao[d_nelc]= kTPCnSigmaKao[itrk];
 919   d_TPCnSigmaPro[d_nelc]= kTPCnSigmaPro[itrk];
 920   d_TOFnSigmaEle[d_nelc]= kTOFnSigmaEle[itrk];
 921   d_TOFnSigmaPio[d_nelc]= kTOFnSigmaPio[itrk];
 922   d_TOFnSigmaMuo[d_nelc]= kTOFnSigmaMuo[itrk];
 923   d_TOFnSigmaKao[d_nelc]= kTOFnSigmaKao[itrk];
 924   d_TOFnSigmaPro[d_nelc]= kTOFnSigmaPro[itrk];
 925
 926   d_chi2ndf[d_nelc]= kChi2NDF[itrk];
 927   d_E[d_nelc] = kE[itrk];
 928   d_dphi[d_nelc] = kDeltaPhi[itrk];
 929   d_deta[d_nelc] = kDeltaEta[itrk];
 930
 931   d_nelc++;
 932
 933 }
 934
 935 //____________________________________________________________________
 936 void ana_sgl::randomize_pool(vector<etrk> e1, vector<etrk> e2){
 937
 938   int size1 = e1.size();
 939   int used_index[1000];
 940   for(int i=0;i<1000;i++){
 941     used_index[i] = -1;
 942   }
 943   for(int i=0;i<size1;i++){
 944     used_index[i] = 0;
 945   }
 946
 947   for(unsigned int i=0;i<size1;i++){
 948     int j = (int)(gRandom->Uniform(0,size1));
 949     while(used_index[j]==1){
 950       j = (int)(gRandom->Uniform(0,size1));
 951     }
 952     if( (e1[j].ghost_flag==1) &&
 953         (
 954          (d_conv_flag==true && e1[j].conv_flag==1) ||
 955          (d_conv_flag==false)
 956          )
 957         ){
 958       vep.push_back(e1[j]);
 959     }
 960     used_index[j] = 1;
 961   }
 962
 963
 964   int size2 = e2.size();
 965   for(int i=0;i<1000;i++){
 966     used_index[i] = -1;
 967   }
 968   for(int i=0;i<size2;i++){
 969     used_index[i] = 0;
 970   }
 971
 972   for(unsigned int i=0;i<size2;i++){
 973     int j = (int)(gRandom->Uniform(0,size2));
 974     while(used_index[j]==1){
 975       j = (int)(gRandom->Uniform(0,size2));
 976     }
 977     if( (e2[j].ghost_flag==1) &&
 978       (
 979        (d_conv_flag==true && e2[j].conv_flag==1) ||
 980        (d_conv_flag==false)
 981        )
 982         ){
 983       vem.push_back(e2[j]);
 984     }
 985     used_index[j] = 1;
 986   }
 987 }
 988
 989 //____________________________________________________________________
 990 void ana_sgl::calc_pair(vector<etrk> vep, vector<etrk> vem){
 991
 992   //vector<etrk>::iterator iep, iem;
 993
 994   vector<etrk>::iterator iep;
 995   vector<etrk>::iterator iem;
 996
 997   //cout<<vep.size()<<" "<<vem.size()<<endl;
 998
 999   ///unlike pairs
1000   for(iep=vep.begin(); iep!=vep.end(); ++iep){
1001     for(iem=vem.begin(); iem!=vem.end(); ++iem){
1002       if(PairTrackcut(iep, iem)==false) continue;
1003       fill_pair(iep, iem, 0);
1004     }
1005   }
1006   for(iep=vep.begin(); iep!=vep.end(); ++iep){
1007     vector<etrk>::iterator iep2=iep;
1008     ++iep2;
1009     for(iem=iep2; iem!=vep.end(); ++iem){
1010       if(PairTrackcut(iep, iem)==false) continue;
1011       fill_pair(iep, iem, 1);
1012     }
1013   }
1014
1015   for(iep=vem.begin(); iep!=vem.end(); ++iep){
1016     vector<etrk>::iterator iep2=iep;
1017     ++iep2;
1018     for(iem=iep2; iem!=vem.end(); ++iem){
1019       if(PairTrackcut(iep, iem)==false) continue;
1020       fill_pair(iep, iem, 2);
1021     }
1022   }
1023
1024   int icent = (int)(fkCentrality/10.0);
1025   int izbin = (int)((fkZvPrim+10)/2.0);
1026   if(icent<0) icent=0;
1027   if(icent>=NCENT) icent=NCENT-1;
1028   if(izbin<0) izbin=0;
1029   if(izbin>=NZBIN) izbin=NZBIN-1;
1030
1031
1032   int nmixed;
1033   if(vep.size()>0) {
1034     //
1035     // Now mixed event for +- pairs
1036     //
1037     nmixed = 0;
1038     for(int ibuf=0;(nmixed<NMix);ibuf++) {
1039       int ntry = 0;
1040       while(ntry<MAX_TRY) {
1041         reshuffle_buffer(d_vem[ibuf][izbin][icent],d_poolm[izbin][icent]);
1042         ntry++;
1043       }
1044       for(iep=vep.begin(); iep!=vep.end(); ++iep){
1045         for(iem=d_vem[ibuf][izbin][icent].begin();
1046             iem!=d_vem[ibuf][izbin][icent].end(); ++iem){
1047           if(PairTrackcut(iep, iem)==false) continue;
1048           fill_pair(iep,iem,3);
1049         }
1050       }
1051       ++nmixed;
1052     }//for(ibuf)
1053   }
1054
1055   if(vem.size()>0) {
1056     //
1057     // Now mixed event for +- pairs
1058     //
1059     nmixed = 0;
1060     for(int ibuf=0;(nmixed<NMix);ibuf++) {
1061       int ntry = 0;
1062       while(ntry<MAX_TRY) {
1063         reshuffle_buffer(d_vep[ibuf][izbin][icent],d_poolp[izbin][icent]);
1064         ntry++;
1065       }
1066       for(iem=vem.begin(); iem!=vem.end(); ++iem){
1067         for(iep=d_vep[ibuf][izbin][icent].begin();
1068             iep!=d_vep[ibuf][izbin][icent].end(); ++iep){
1069           if(PairTrackcut(iep, iem)==false) continue;
1070           fill_pair(iep,iem,4);
1071         }
1072       }
1073       ++nmixed;
1074     }//for(ibuf)
1075   }
1076
1077
1078   if(vep.size()>0) {
1079     //
1080     // Now mixed event for ++ pairs
1081     //
1082     nmixed = 0;
1083     for(int ibuf=0;(nmixed<NMix);ibuf++) {
1084       int ntry = 0;
1085       while(ntry<MAX_TRY) {
1086         reshuffle_buffer(d_vep[ibuf][izbin][icent],d_poolp[izbin][icent]);
1087         ntry++;
1088       }
1089       for(iep=vep.begin(); iep!=vep.end(); ++iep){
1090         for(iem=d_vep[ibuf][izbin][icent].begin();
1091             iem!=d_vep[ibuf][izbin][icent].end(); ++iem){
1092           if(PairTrackcut(iep, iem)==false) continue;
1093           fill_pair(iep,iem,5);
1094         }
1095       }
1096       ++nmixed;
1097     }//for(ibuf)
1098   }
1099
1100   if(vem.size()>0) {
1101     //
1102     // Now mixed event for +- pairs
1103     //
1104     nmixed = 0;
1105     for(int ibuf=0;(nmixed<NMix);ibuf++) {
1106       int ntry = 0;
1107       while(ntry<MAX_TRY) {
1108         reshuffle_buffer(d_vem[ibuf][izbin][icent],d_poolm[izbin][icent]);
1109         ntry++;
1110       }
1111       for(iem=vem.begin(); iem!=vem.end(); ++iem){
1112         for(iep=d_vem[ibuf][izbin][icent].begin();
1113             iep!=d_vem[ibuf][izbin][icent].end(); ++iep){
1114           if(PairTrackcut(iep, iem)==false) continue;
1115           fill_pair(iep,iem,6);
1116         }
1117       }
1118       ++nmixed;
1119     }//for(ibuf)
1120   }
1121 }
1122
1123
1124 //___________________________________________________________________________
1125 //void ana_sgl::fill_pair(etrk* iep, etrk* iem, int type){
1126 //void ana_sgl::fill_pair(etrk iep, etrk iem, int type){
1127 void ana_sgl::fill_pair(vector<etrk>::iterator iep, vector<etrk>::iterator iem, int type){
1128
1129
1130   d_pairtype = type;
1131
1132   calc_vars(iep, iem, d_mass, d_phiv, d_pxpair, d_pypair, d_pzpair,
1133             d_ptpair, d_epair, d_phipair, d_etapair, d_cos, d_psi);
1134
1135   if(type==0||type==1||type==2||type==3||type==5){
1136     d_centrality = iep->cent;
1137     d_prim_xv = iep->pxv;
1138     d_prim_yv = iep->pyv;
1139     d_prim_zv = iep->pzv;
1140   }else if(type==4 || type==6){
1141     d_centrality = iem->cent;
1142     d_prim_xv = iem->pxv;
1143     d_prim_yv = iem->pyv;
1144     d_prim_zv = iem->pzv;
1145
1146   }
1147
1148
1149   //cout<<iep->px<<" "<<iem->px<<" "<<d_centrality<<endl;
1150
1151   d_cent1 = iep->cent;
1152   d_xv1 = iep->xv;
1153   d_yv1 = iep->yv;
1154   d_zv1 = iep->zv;
1155   d_px1 = iep->px;
1156   d_py1 = iep->py;
1157   d_pz1 = iep->pz;
1158   d_pt1 = iep->pt;
1159   d_eta1 = iep->eta;
1160   d_phi1 = iep->phi;
1161   d_theta1 = iep->theta;
1162   d_tpc1 = iep->tpc;
1163   d_ntpc_ele1 = iep->ntpc_ele;
1164   d_ntpc_pio1 = iep->ntpc_pio;
1165   d_ntpc_kao1 = iep->ntpc_kao;
1166   d_ntpc_pro1 = iep->ntpc_pro;
1167   d_ntof_ele1 = iep->ntof_ele;
1168   d_ntof_pio1 = iep->ntof_pio;
1169   d_ntof_kao1 = iep->ntof_kao;
1170   d_ntof_pro1 = iep->ntof_pro;
1171   d_its1 = iep->its;
1172   d_nits1 = iep->nits;
1173   d_ntpc1 = iep->ntpc;
1174   d_e1 = iep->e;
1175   d_dphi1 = iep->dphi;
1176   d_deta1 = iep->deta;
1177   d_dcaxy1 = iep->dcaxy;
1178   d_dcaz1 = iep->dcaz;
1179   d_conv1 = iep->conv_flag;
1180
1181   d_cent2 = iem->cent;
1182   d_xv2 = iem->xv;
1183   d_yv2 = iem->yv;
1184   d_zv2 = iem->zv;
1185   d_px2 = iem->px;
1186   d_py2 = iem->py;
1187   d_pz2 = iem->pz;
1188   d_pt2 = iem->pt;
1189   d_eta2 = iem->eta;
1190   d_phi2 = iem->phi;
1191   d_theta2 = iem->theta;
1192   d_tpc2 = iem->tpc;
1193   d_ntpc_ele2 = iem->ntpc_ele;
1194   d_ntpc_pio2 = iem->ntpc_pio;
1195   d_ntpc_kao2 = iem->ntpc_kao;
1196   d_ntpc_pro2 = iem->ntpc_pro;
1197   d_ntof_ele2 = iem->ntof_ele;
1198   d_ntof_pio2 = iem->ntof_pio;
1199   d_ntof_kao2 = iem->ntof_kao;
1200   d_ntof_pro2 = iem->ntof_pro;
1201   d_its2 = iem->its;
1202   d_nits2 = iem->nits;
1203   d_ntpc2 = iem->ntpc;
1204   d_e2 = iem->e;
1205   d_dphi2 = iem->dphi;
1206   d_deta2 = iem->deta;
1207   d_dcaxy2 = iem->dcaxy;
1208   d_dcaz2 = iem->dcaz;
1209   d_conv2 = iem->conv_flag;
1210
1211   d_run = fkRunNumber;
1212   //  d_event = d_evt;
1213
1214
1215   int icent = (int)(fkCentrality/10.0);
1216   if(icent<0) icent=0;
1217   if(icent>=NCENT) icent=NCENT-1;
1218   if(pair_cut()==true){
1219     hmasspt[type][icent]->Fill(d_mass, d_ptpair);
1220     hmasspt[type][10]->Fill(d_mass, d_ptpair);
1221   }
1222
1223   ////////////// pt weighting for mixed event
1224   float weight = 1;
1225   float k = 10000;
1226   if(type==0 || type==1 || type==2){
1227     weight = 1;
1228   }else if(type==3 || type==5 || type==6){
1229     k = 18.7*exp(-2.4*d_pt2)+1.24;
1230     weight = 1+1/k;
1231   }else if(type==4){
1232     k = 18.7*exp(-2.4*d_pt1)+1.24;
1233     weight = 1+1/k;
1234   }
1235   if(pair_cut()==true){
1236     hmasspt_weight[type][icent]->Fill(d_mass, d_ptpair, weight);
1237     hmasspt_weight[type][10]->Fill(d_mass, d_ptpair, weight);
1238   }
1239
1240
1241   //d_ntpair->Fill();
1242
1243 }
1244
1245 //____________________________________________
1246 bool ana_sgl::PairTrackcut(vector<etrk>::iterator e1, vector<etrk>::iterator e2){
1247   /*
1248   double p1 = sqrt(pow(e1.px,2)+
1249                    pow(e1.py,2)+
1250                    pow(e1.pz,2));
1251
1252   double p2 = sqrt(pow(e2.px,2)+
1253                    pow(e2.py,2)+
1254                    pow(e2.pz,2));
1255
1256   if(
1257      (e1.e/p1>0.7 && e1.e/p1<1.3) ||
1258      (e2.e/p1>0.7 && e2.e/p2<1.3)
1259      ){
1260     return true;
1261   }else{
1262     return false;
1263   }
1264 o  */
1265   return true;
1266
1267 }
1268
1269 //____________________________________________
1270 void ana_sgl::check_conversion_pairs(vector<etrk> &e1, vector<etrk> &e2){
1271   vector<etrk>::iterator iep;
1272   vector<etrk>::iterator iem;
1273   bool reject = false;
1274   if(e1.size()>0 && e2.size()>0){
1275     for(iep = e1.begin(); iep != e1.end(); ++iep){
1276       reject = false;
1277       for(iem = e2.begin(); iem != e2.end(); ++iem){
1278         double mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi;
1279         calc_vars(iep, iem, mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi);
1280         if(magnetic_field_mm==true){ //pp
1281           if(phiv<0.6 && iep->phi-iem->phi<0){ // this depends on the magntic field
1282             reject = true;
1283             iem->conv_flag = 0;
1284           }
1285         }else{
1286           if(phiv>acos(-1.0)-0.6 && iep->phi-iem->phi>0){ // this depends on the magntic field
1287             reject = true;
1288             iem->conv_flag = 0;
1289           }
1290         }
1291       }
1292       if(reject==true) iep->conv_flag=0;
1293     }
1294   }
1295 }
1296
1297 //____________________________________________
1298 void ana_sgl::check_ghost_pairs(vector<etrk> &e1){
1299   vector<etrk>::iterator iep;
1300   vector<etrk>::iterator iem;
1301   bool reject = false;
1302   if(e1.size()>1){
1303     for(iep = e1.begin(); iep != e1.end(); ++iep){
1304       reject = false;
1305       vector<etrk>::iterator iep2=iep;
1306       ++iep2;
1307       for(iem = iep2; iem != e1.end(); ++iem){
1308         double mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi;
1309         calc_vars(iep, iem, mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi);
1310         if(mass<0.01){
1311           reject = true;
1312           iem->ghost_flag = 0;
1313         }
1314       }
1315       if(reject==true) iep->ghost_flag=0;
1316     }
1317   }
1318 }
1319
1320 //____________________________________________
1321 //void ana_sgl::calc_vars(etrk iep, etrk iem, double &mass, double &phiv, double &px, double &py, double&pz,
1322 void ana_sgl::calc_vars(vector<etrk>::iterator iep, vector<etrk>::iterator iem, double &mass, double &phiv, double &px, double &py, double&pz,
1323                         double &pt, double &e, double &phi, double &eta, double &cos, double &psi){
1324
1325   px = iep->px+iem->px;
1326   py = iep->py+iem->py;
1327   pz = iep->pz+iem->pz;
1328   pt = sqrt(px*px+py*py);
1329   double d_ppair = sqrt(pt*pt+pz*pz);
1330   static const double me=0.0005109989;
1331   e = sqrt(me*me+iep->px*iep->px+iep->py*iep->py+iep->pz*iep->pz)
1332     + sqrt(me*me+iem->px*iem->px+iem->py*iem->py+iem->pz*iem->pz);
1333
1334   mass =  e*e-px*px-py*py-pz*pz;
1335   if(mass<0){
1336     mass = mass;
1337   }else{
1338     mass = sqrt(mass);
1339   }
1340
1341
1342
1343   phi = atan2(py, px);
1344   eta = -0.5*TMath::Log((d_ppair+pz)/(d_ppair-pz));
1345   double p1 = sqrt(pow(iep->px,2)+pow(iep->py,2)+pow(iep->pz,2));
1346   double p2 = sqrt(pow(iem->px,2)+pow(iem->py,2)+pow(iem->pz,2));
1347   cos = acos((iep->px*iem->px+iep->py*iem->py+iep->pz*iem->pz)/(p1*p2));
1348
1349
1350   double dtheta = iep->theta-iem->theta;
1351   psi = asin(dtheta/cos);
1352
1353
1354   //unit vector of (pep+pem)
1355   float pl = d_ppair;
1356   float ux = px/pl;
1357   float uy = py/pl;
1358   float uz = pz/pl;
1359   float ax = uy/sqrt(ux*ux+uy*uy);
1360   float ay = -ux/sqrt(ux*ux+uy*uy);
1361
1362   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
1363   //definition.
1364   float ptep = iep->px*ax + iep->py*ay;
1365   float ptem = iem->px*ax + iem->py*ay;
1366
1367   float pxep = iep->px;
1368   float pyep = iep->py;
1369   float pzep = iep->pz;
1370   float pxem = iem->px;
1371   float pyem = iem->py;
1372   float pzem = iem->pz;
1373
1374
1375   //vector product of pep X pem
1376   float vpx = pyep*pzem - pzep*pyem;
1377   float vpy = pzep*pxem - pxep*pzem;
1378   float vpz = pxep*pyem - pyep*pxem;
1379   float vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
1380   float thev = acos(vpz/vp);
1381
1382   //unit vector of pep X pem
1383   float vx = vpx/vp;
1384   float vy = vpy/vp;
1385   float vz = vpz/vp;
1386
1387   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
1388   float wx = uy*vz - uz*vy;
1389   float wy = uz*vx - ux*vz;
1390   float wz = ux*vy - uy*vx;
1391   float wl = sqrt(wx*wx+wy*wy+wz*wz);
1392   // by construction, (wx,wy,wz) must be a unit vector.
1393   if(fabs(wl - 1.0) > 0.00001) cout << "Calculation error in W vector"<<endl;
1394   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
1395   // should be small if the pair is conversion
1396   //
1397   float cosPhiV = wx*ax + wy*ay;
1398   phiv = acos(cosPhiV);
1399
1400 }
1401
1402 //____________________________________________
1403 bool ana_sgl::pair_cut(void){
1404   bool ret = true;
1405
1406
1407   if(d_flag_phiv==true){
1408     if(d_run>169591 && d_phiv>acos(-1.0)-d_phiv_cut){
1409       ret = false;
1410     }
1411     if(d_run<169591 && d_phiv<d_phiv_cut){
1412       ret = false;
1413     }
1414   }
1415
1416
1417   if(d_flag_emc_cut==true){
1418     if( !(
1419           (d_e1/d_pt1>d_emc_low && d_e1/d_pt1<d_emc_high) ||
1420           (d_e2/d_pt2>d_emc_low && d_e2/d_pt2<d_emc_high)
1421           )
1422         ){
1423       ret = false;
1424     }
1425   }
1426
1427
1428   if(d_flag_pt_cut==true){
1429     if(!(d_pt1>d_pt_cut_low && d_pt2>d_pt_cut_low &&
1430          d_pt1<d_pt_cut_high && d_pt2<d_pt_cut_high)
1431        ){
1432       ret = false;
1433     }
1434   }
1435
1436   return ret;
1437
1438 }
1439
1440 //____________________________________________
1441 void ana_sgl::print_cuts(void){
1442
1443   cout<<" ********** lists of cuts ********** "<<endl;
1444   cout<<d_tpc_dedx_low<<" < TPC dE/dx < "<<d_tpc_dedx_high<<endl;
1445
1446
1447   if(d_flag_kaon_veto==true){
1448     cout<<d_dedx_kaon_veto_low<<" < TPC veto for Kaon < "<<d_dedx_kaon_veto_high<<endl;
1449   }else{
1450     cout<<"  No TPC Kaon Veto "<<endl;
1451   }
1452
1453   if(d_flag_proton_veto==true){
1454     cout<<d_dedx_proton_veto_low<<" < TPC veto for proton < "<<d_dedx_proton_veto_high<<endl;
1455   }else{
1456     cout<<" No TPC Proton Veto "<<endl;
1457   }
1458
1459   if(d_flag_tof_cut==true){
1460     cout<<d_tof_low<<" < kTOFnSigmaEle < "<<d_tof_high<<endl;
1461   }else{
1462     cout<<" No TOF cuts "<<endl;
1463   }
1464   if(d_flag_emc_cut==true){
1465     cout<<d_emc_low<<" < EMC E/p (or for pairs) < "<<d_emc_high<<endl;
1466   }else{
1467     cout<<" No EMC cuts "<<endl;
1468   }
1469   if(d_flag_phiv==true){
1470     cout<<" Phiv cuts "<<d_phiv_cut<<endl;
1471   }else{
1472     cout<<" No phiv cuts "<<endl;
1473   }
1474
1475   if(d_flag_pt_cut==true){
1476     cout<<d_pt_cut_low<<" < pT of electrons < "<<d_pt_cut_high<<endl;
1477   }else{
1478     cout<<" No electron pt cuts "<<endl;
1479   }
1480
1481   if(d_conv_flag==true){
1482     cout<<" remove converson pairs from electron/positron lists "<<endl;
1483   }else{
1484     cout<<" conversion candidates are in the pool."<<endl;
1485   }
1486   cout<<" ********** end of lists of cuts ********** "<<endl;
1487
1488   if(sel_trigger==0){
1489     cout<<" MB "<<endl;
1490   }else if(sel_trigger==1){
1491     cout<<" SemiCentral "<<endl;
1492   }else if(sel_trigger==2){
1493     cout<<" Central "<<endl;
1494   }
1495
1496 }
1497
1498