Attachment 'lifetime.c'

Download

   1 #include<TStyle.h>
   2 #include<TFile.h>
   3 #include<TTree.h>
   4 #include<THStack.h>
   5 #include<math.h>
   6 #include<TH1.h>
   7 #include<RooFit.h>
   8 #include<RooPlot.h>
   9 #include<RooAbsData.h>
  10 #include<RooDataHist.h>
  11 #include<RooDataSet.h>
  12 #include<RooRealVar.h>
  13 #include<RooGaussModel.h>
  14 #include<RooDecay.h>
  15 
  16 /******************************************************************/
  17 /* lifetime_rooFit calculates the neutralino lifetime,            */
  18 /* using the reconstructed photonshower and kinematic constraints */
  19 /*                                                                */
  20 /* The fit, a gaussian, convoluted with an exponential            */
  21 /* is realized with rooFit                                        */
  22 /******************************************************************/
  23 
  24 const float c = 299792458;   //speed of light [m/s]
  25 const float E_cms  = 500.0;  //center of mass energy [GeV]
  26 const float E_0    = 250.0;  //beam energy [GeV]
  27 const float M_neu1 = 151.0;  //neutralino mass [GeV]
  28 
  29 const int MAXFILES   = 1; //maximum # cluster per event
  30 const int MAXCLUSTER = 30; //maximum # cluster per event
  31 
  32 using namespace RooFit;
  33 
  34 void lifetime_rooFit2() {
  35 
  36   float E_neu1    = E_cms/2;                                 //neutralino energy [GeV]
  37   float P_0       = sqrt( pow(E_0,2) - pow(M_neu1,2) );      //neutralino momentum [GeV]
  38   float betaGamma = sqrt( pow(E_neu1,2)/pow(M_neu1,2) - 1 ); //beta*gamma
  39 
  40   // give input root tree
  41   TFile *myFile[MAXFILES];
  42   TString fileName[MAXFILES] = {
  43     /*    "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_neu1neu1_tau11ns.root"
  44     "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n1n1aa.root",
  45     "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n2n2aa.root",
  46     "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa.root",
  47     "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_aa.root"*/
  48     "/data/nwattime/LDC_optimization/test.root"
  49   };
  50 
  51   float nExpected[MAXFILES] = {100};//{2001406,70319,70674,3113697};//{216728};//,2001406,70319,70674,3113697};
  52   float efficiency[MAXFILES] ={1};//{1,1,1,1};// {0.5948};//{0.0143,0.0234,0.0233,0.0004};
  53   //0.2ns: 0.5948 //2ns: 0.5368 //11ns: 0.1018
  54 
  55   float _weight;  float _time;
  56 
  57   // helper tree
  58   TFile* _helperFile = new TFile("helperFile","recreate");
  59   TTree* _helperTree = new TTree("helperTree","helperTree");
  60   _helperTree->SetAutoSave();
  61   _helperTree->Branch("time",&_time,"time/F");
  62   _helperTree->Branch("weight",&_weight,"weight/F");
  63 
  64   //DEBUG: plot psi
  65   TH1F *psiHisto = new TH1F("psi","psi",100,0,600);
  66 
  67   for (int fileCounter=0; fileCounter<MAXFILES; fileCounter++) {
  68     myFile[fileCounter] = new TFile(fileName[fileCounter],"READ");
  69 
  70     if ( myFile[fileCounter]->IsOpen() ) {
  71 
  72       int   nCluster;
  73       float _Egamma[MAXCLUSTER];
  74       float _phi[MAXCLUSTER];    float _psi[MAXCLUSTER];
  75       float _pos[MAXCLUSTER][3]; float _axis[MAXCLUSTER][3];
  76       float _cog[MAXCLUSTER];    float _dir[MAXCLUSTER];
  77       float _dist[MAXCLUSTER];
  78 
  79       // open input tree
  80       TTree* _tree = (TTree*)myFile[fileCounter]->Get("tree");
  81       _tree->SetBranchAddress("cluster_nCluster",    &nCluster);
  82       _tree->SetBranchAddress("cluster_energy",      &_Egamma); //photon energy [GeV]
  83       _tree->SetBranchAddress("cluster_clusterAxis", &_axis);   //principle axis of inertia [mm]
  84       _tree->SetBranchAddress("cluster_position",    &_pos);    //centre of gravity [mm]     
  85 
  86       float nEntries = _tree->GetEntriesFast();
  87       _weight = nExpected[fileCounter]*efficiency[fileCounter]/nEntries;
  88       std::cout<<"for file: "<<fileName[fileCounter]<<" found weight "<<_weight<<std::endl;
  89       for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
  90 	_tree->GetEntry(entryCounter);
  91 
  92 	std::cout<<"Found "<<nCluster<<" clusters in event: "<<entryCounter<<std::endl;
  93 
  94 	// calculate reconstructed lifetime
  95 	for (int clusterCounter=0; clusterCounter<nCluster; ++clusterCounter) {
  96 
  97 	  // |centre of gravity|
  98 	  _cog[clusterCounter] = sqrt( pow(_pos[clusterCounter][0],2)
  99 				     + pow(_pos[clusterCounter][1],2)
 100 				     + pow(_pos[clusterCounter][2],2) );
 101 	  // |principle axis|
 102 	  _dir[clusterCounter] = sqrt( pow(_axis[clusterCounter][0],2)
 103 				     + pow(_axis[clusterCounter][1],2)
 104 				     + pow(_axis[clusterCounter][2],2) );
 105 	  std::cout<<"cluster cog: "<<_cog[clusterCounter]<<endl;
 106 	  // <(axis,cog) -> cos(psi)= (axis*cog) / (|axis|*|cog|)
 107 	  _psi[clusterCounter] = acos( ( _pos[clusterCounter][0]   * _axis[clusterCounter][0]
 108 					 + _pos[clusterCounter][1] * _axis[clusterCounter][1]
 109 					 + _pos[clusterCounter][2] * _axis[clusterCounter][2] )
 110 				       / ( _cog[clusterCounter]    * _dir[clusterCounter] ) );
 111 	  // <(neu1,gamma) -> cos(phi) = E0/p0 - m_neu1^2/(2*p0*E_gamma)
 112 
 113 	  if ( _Egamma[clusterCounter]>0 )
 114 	    _phi[clusterCounter] = acos( E_0/P_0
 115 					 - pow(M_neu1,2) / (2*P_0*_Egamma[clusterCounter]) );
 116 	  else
 117 	    _phi[clusterCounter] = 0;
 118 
 119 	  //neutralino flight legth lamda = |cog|*sin(psi)/sin(phi)
 120 	  if ( sin(_phi[clusterCounter]) !=0 )
 121 	    _dist[clusterCounter] = _cog[clusterCounter]
 122 	      * sin(_psi[clusterCounter]) / sin(_phi[clusterCounter]); //[mm]
 123 	  else
 124 	    _dist[clusterCounter] = 0;
 125 
 126 	  _time = _dist[clusterCounter] * 1e6 / (c*betaGamma);
 127 	  _helperTree->Fill();
 128 
 129 	  std::cout<<"flight dist: "<<_dist[clusterCounter]<<endl;
 130 	  std::cout<<"flight time: "<<_time<<endl;
 131 
 132 	  //DEBUG: plot psi
 133 	  psiHisto->Fill(_dist[clusterCounter]);
 134 	}//all reco cluster
 135       } //all entries
 136     }//file open
 137   }//all files
 138 
 139   //DEBUG: plot psi
 140   psiHisto->Draw();
 141 
 142   TFile *tree_file = _helperTree->GetCurrentFile();
 143   tree_file->Write();
 144   
 145   /********* draw neutralino lifetime spectrum *********/
 146 
 147   RooRealVar time("time","decay time [ns]",0,10);
 148   RooRealVar weight("weight","event weight",0,250);
 149   //    RooDataSet *data = new RooDataSet("data","data",_helperTree,RooArgSet(time,weight),WeightVar(weight));
 150   RooDataSet *data = new RooDataSet("data","data",RooArgSet(time,weight),0,WeightVar(weight));
 151   //    data->setWeightVar(weight);
 152     data->Print("v");
 153     
 154     RooPlot* frame = time.frame();
 155     data->plotOn(frame,DataError(RooAbsData::SumW2));
 156 
 157     frame->SetTitle("");
 158     frame->GetYaxis()->SetTitle("# #chi^{0}_{1}");
 159     /*
 160     //Decay model (Exponential convoluted with Gauss)
 161     //name,title,start value,minimum,maximum 
 162     RooRealVar mean("mean","resolution mean",    0.270, 0.270, 0.270);  //.2,0.,5.
 163     RooRealVar sigma("sigma","resolution sigma", 0.185, 0.185, 0.185); //.1,0.,5.
 164     RooRealVar decayConst("decayConst","decay constant", 1.0, 0.01, 15.0); //1.,0.,15.
 165     RooRealVar fsig("fsig","signal fraction", 0.5, 0.1, 1.0); //0.5,0.,1.
 166 
 167     //for an observable time with uncertainty dt for each event
 168     //observable,mean,sigma,observable error
 169     RooGaussModel res("res","det. resol.",time,mean,sigma);
 170     RooDecay      decay("decay","decay model",time,decayConst,res,RooDecay::SingleSided);
 171 
 172     //Gauss model for background: same gaussian as detector resolution
 173     //observable,mean,sigma,observable error
 174     RooGaussModel BgModel("BgModel","background model.",time,mean,sigma);
 175 
 176     //sum signal and background model
 177     //model(x) = fsig*decay(x) + (1-fsig)*BgModel(x)
 178     RooAddPdf  model("model","model",RooArgList(decay,BgModel),fsig);
 179     RooFitResult* r = model->fitTo( *data,
 180 				    Minos(false),
 181 				    Range(0.0, 5.0),
 182 				    Save(true));
 183 
 184 				    model->plotOn(frame);*/
 185     frame->Draw();
 186 
 187 }//lifetime2
 188 

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2010-03-23 14:46:21, 14.1 KB) [[attachment:ClusterWriteEngine.cc]]
  • [get | view] (2010-03-23 14:47:51, 4.6 KB) [[attachment:ClusterWriteEngine.hh]]
  • [get | view] (2010-03-23 14:57:17, 17.4 KB) [[attachment:eventSelection.c]]
  • [get | view] (2010-03-23 15:40:25, 7.3 KB) [[attachment:lifetime.c]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.