MUON/masse.C

   1 void masse() {
   2 //////////////////////////////////////////////////////////
   3 //   This file has been automatically generated
   4 //     (Wed Jun 16 15:35:08 1999 by ROOT version 2.21/07)
   5 //   from TTree ntuple/Reconst ntuple
   6 //   found on file: reconst.root
   7 //////////////////////////////////////////////////////////
   8
   9 //Reset ROOT and connect tree file
  10 gROOT->Reset();
  11
  12 TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("reconst.root");
  13 if (!f) {
  14      f = new TFile("reconst.root");
  15 }
  16
  17 TTree *ntuple = (TTree*)gDirectory->Get("ntuple");
  18
  19 TPostScript ps("mass.ps",112);
  20 //ps->Open("mass.ps");
  21 //Declaration of leaves types
  22 Int_t           ntrack=500;
  23
  24 Int_t           ievr;
  25 Int_t           ntrackr;
  26 Int_t           istatr[ntrack];
  27 Int_t           isignr[ntrack];
  28 Float_t         pxr[ntrack];
  29 Float_t         pyr[ntrack];
  30 Float_t         pzr[ntrack];
  31 Float_t         zvr[ntrack];
  32 Float_t         chi2r[ntrack];
  33
  34 //Set branch addresses
  35 ntuple->SetBranchAddress("ievr",&ievr);
  36 ntuple->SetBranchAddress("ntrackr",&ntrackr);
  37 ntuple->SetBranchAddress("istatr",&istatr);
  38 ntuple->SetBranchAddress("isignr",&isignr);
  39 ntuple->SetBranchAddress("pxr",&pxr);
  40 ntuple->SetBranchAddress("pyr",&pyr);
  41 ntuple->SetBranchAddress("pzr",&pzr);
  42 ntuple->SetBranchAddress("zvr",&zvr);
  43 ntuple->SetBranchAddress("chi2r",&chi2r);
  44
  45 //     This is the loop skeleton
  46 //       To read only selected branches, Insert statements like:
  47 // ntuple->SetBranchStatus("*",0);  // disable all branches
  48 // ntuple->SetBranchStatus("branchname",1);  // activate branchname
  49
  50 TH1F *h1 = new TH1F("h1","mass",240,0.,12.);
  51 TH1F *h2 = new TH1F("h2","mass good muon",240,0.,12.);
  52 // Upsilon  mass cut
  53 TH1F *h5 = new TH1F("h5"," Mass Upsilon",70,8.,11.);
  54 // J/psi  mass cut
  55 //TH1F *h5 = new TH1F("h5"," Mass J/psi",120,1.,5.);
  56 TH1F *h3 = new TH1F("h3","Pt",100,0.,20.);
  57 TH1F *h4 = new TH1F("h4","Pt bon muon",100,0.,20.);
  58
  59 gStyle->SetOptFit();
  60 //gStyle->SetOptStat(0000);
  61 gStyle->SetOptStat(10001110);
  62
  63 Float_t mass = 1.;
  64 Float_t amassmu = 0.10566;
  65 Float_t Chi2Cut = 100.;
  66 Float_t PtCut = 0.;
  67 //Float_t MassMin = 8.96;
  68 //Float_t MassMax = 9.96;
  69 // Upsilon
  70 Float_t MassMin = 8.5;
  71 Float_t MassMax = 10.5;
  72 // J/psi
  73 //Float_t MassMin = 2.5;
  74 //Float_t MassMax = 3.8;
  75
  76 //Float_t Chi2Cut = 999.;
  77
  78 Int_t nentries = ntuple->GetEntries();
  79
  80 Int_t nbytes = 0;
  81 Int_t nbtrack[50]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  82
  83 Int_t compt = 0;
  84 Int_t Nres = 0;
  85
  86 //    Good muons
  87 for (Int_t iev=0; iev<nentries; iev++) {
  88     nbytes += ntuple->GetEvent(iev);
  89 //    ntuple->SetBranchStatus("main",1);
  90     ntuple->SetBranchStatus("*",1);
  91
  92     //printf("\n");
  93     //  printf("ntrackr=%d\n",ntrackr);
  94
  95     if(ntrackr < 50)   nbtrack[ntrackr]=nbtrack[ntrackr]+1 ;
  96
  97
  98     for (Int_t i=0; i<(ntrackr-1); i++) {
  99
 100       Float_t px1 =  pxr[i] ;
 101       Float_t py1 =  pyr[i] ;
 102       Float_t pz1 =  pzr[i] ;
 103       Float_t Pt1 =  TMath::Sqrt(px1*px1+py1*py1) ;
 104
 105 //      if (chi2r[i] < Chi2Cut && istatr[i] == 1 && Pt1 > PtCut) {
 106       h4->Fill(Pt1);
 107
 108       for (Int_t j=i+1; j<ntrackr; j++) {
 109         Float_t px2 =  pxr[j] ;
 110         Float_t py2 =  pyr[j] ;
 111         Float_t pz2 =  pzr[j] ;
 112         Float_t Pt2 =  TMath::Sqrt(px2*px2+py2*py2) ;
 113         if (chi2r[j]<Chi2Cut && istatr[j] == 1 && Pt2 > PtCut ) {
 114             if (isignr[i]!=isignr[j]) {
 115                 Float_t p12 = px1*px1+py1*py1+pz1*pz1;
 116                 Float_t p22 = px2*px2+py2*py2+pz2*pz2;
 117                 Float_t e1 = TMath::Sqrt(amassmu*amassmu+p12);
 118                 Float_t e2 = TMath::Sqrt(amassmu*amassmu+p22);
 119                 Float_t amass = TMath::Sqrt((e1+e2)*(e1+e2)-(px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
 120 //              printf("amass=%f\n",amass);
 121                 h2->Fill(amass);
 122                 if (amass<11. && amass>8.) compt++;
 123             }
 124         }
 125       }
 126     }
 127
 128
 129
 130
 131     //   All tracks
 132     for (Int_t i=0; i<(ntrackr-1); i++) {
 133       Float_t px1 =  pxr[i] ;
 134       Float_t py1 =  pyr[i] ;
 135       Float_t pz1 =  pzr[i] ;
 136       Float_t Pt1 =  TMath::Sqrt(px1*px1+py1*py1) ;
 137       if (chi2r[i] < Chi2Cut && Pt1 > PtCut) {
 138         h3->Fill(Pt1);
 139         for (Int_t j=i+1; j<ntrackr; j++) {
 140           Float_t px2 =  pxr[j] ;
 141           Float_t py2 =  pyr[j] ;
 142           Float_t pz2 =  pzr[j] ;
 143           Float_t Pt2 =  TMath::Sqrt(px2*px2+py2*py2) ;
 144           if (chi2r[j]<Chi2Cut && Pt2 > PtCut) {
 145             if (isignr[i]!=isignr[j]) {
 146               Float_t p12 = px1*px1+py1*py1+pz1*pz1;
 147               Float_t p22 = px2*px2+py2*py2+pz2*pz2;
 148               Float_t e1 = TMath::Sqrt(amassmu*amassmu+p12);
 149               Float_t e2 = TMath::Sqrt(amassmu*amassmu+p22);
 150               Float_t amass = TMath::Sqrt((e1+e2)*(e1+e2)-(px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
 151               //                        printf("amass=%f\n",amass);
 152               h1->Fill(amass);
 153               if (amass > MassMin && amass < MassMax) {
 154                 //              if (ntrackr==2) h5->Fill(amass);
 155                 h5->Fill(amass);
 156 //              if (amass < 9.715 && amass > 9.205)
 157                 if (amass < 9.766 && amass > 9.154)
 158                     {
 159                         printf("amass=%f\n",amass);
 160                         Nres++;
 161                     }
 162               }
 163             }
 164           }
 165         }
 166       }
 167     }
 168 }
 169
 170 //    printf("compt= %d\n",compt);
 171
 172 for (Int_t i = 0 ; i<50 ; i++)
 173 {
 174
 175   printf(" nb of events= %d with %d tracks \n",nbtrack[i], i) ;
 176
 177 }
 178 g1= new TF1("g1","gaus",9.30,9.75) ;
 179 //g1= new TF1("g1","gaus",8.7,10.2) ;
 180 //g1= new TF1("g1","gaus",2.95,3.25) ;
 181 h5->Fit("g1","R0")      ;
 182 h5->SetXTitle("Mass (GeV/c^2!)");
 183 h5->Draw();
 184 printf("Nombre de resonances : %d\n",Nres);
 185
 186 ps->Close();
 187 }
 188
 189
 190
 191
 192
 193
 194
 195