Attachment 'eventSelection.c'

Download

   1 #include <math.h>
   2 #include <TFile.h>
   3 #include <TTree.h>
   4 #include <TGraph.h>
   5 #include<TString.h>
   6 #include<iostream>
   7 
   8 /******** eventSelection processor ********/
   9 /*                                        */
  10 /*  reads nFiles root trees, and applies  */  
  11 /*  cuts defined in selection criteria    */
  12 /*  result is written to outTree          */
  13 /*                                        */
  14 /******************************************/
  15 
  16 const int nFiles = 1;
  17 const int MAXCLUSTER = 30;
  18 
  19 const float E_cms = 500.0;  //center of mass energy [GeV]
  20 
  21 //selection criteria
  22 const int   _minNhitsEcal = 1500; //minimum # hits in ECAL: 1500
  23 const int   _maxNhitsEcal = 6000; //maximum # hits in ECAL: 6000
  24 const float _EminEcal     = 0;   //minimum energy in ECAL: 80 GeV
  25 const float _EmaxEcal     = 700;  //maximum energy in ECAL: 450 GeV
  26 const float _minEt_miss   = 0;   //minimum transverse energy: 30 GeV
  27 const float _minNgamma    = 0;    //minimum # of cluster: 2
  28 const float _minEcluster  = 0;   //minimum energy of cluster: 20 GeV
  29 const float _maxCosTheta  = 1.00; //maximum cluster theta position: 0.75
  30 
  31 using namespace std;
  32 
  33 void eventSelection()
  34 {
  35   TFile *inFile;
  36   TString inFileNames[nFiles] = {
  37     "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa.root"
  38   };
  39   
  40   //  float nExpected[nFiles] = {216728,3113697,2001046,70319,70674};
  41   float nExpected[nFiles] = {70319};
  42   //500fb =  216728 ee_neu1neu1
  43   //        2001046 ee_n1n1aa
  44   //          70319 ee_n2n2aa
  45   //          70674 ee_n3n3aa
  46   //        3113697 ee_aa    
  47   //         180211 ee_eaea
  48   //       73715739 ee_ee  
  49   float polWeight[nFiles] = {1.};//,1.,1.,1.}; 
  50   //0.02 ep+1.0em-1.0, 0.72 ep-1.0em+1.0
  51   //0.40 ep-1.0em+0.0, 0.10 ep+1.0em+0.0, 0.05 ep+0.0em-1.0, 0.45 ep+0.0em+1.0
  52   float weight[nFiles]; // nExpected*polWeight/nTrue;
  53   
  54   /******** create output tree **********/
  55   struct {
  56     int   iEvt;                        // event number             
  57     int   nCluster;                    // # cluster  per event
  58     int   nGamma;                      // # photons per event
  59     float etMiss;                      // missing transverse energy [GeV] 
  60     float esumCalos;                   // energy sum in calorimeters [GeV] 
  61     int   nHitsCalos;                  // # hits in calorimeters 
  62     float esumECAL;                    // energy sum in ECAL [GeV] 
  63     int   nHitsECAL;                   // # hits in ECAL 
  64     float esumHCAL;                    // energy sum in HCAL [GeV] 
  65     int   nHitsHCAL;                   // # hits in HCAL 
  66     float energy[MAXCLUSTER];          // energy per cluster [GeV]
  67     float clusterPos[MAXCLUSTER][3];   // cluster center of gravity 
  68     float clusterAxis[MAXCLUSTER][3];  // cluster axis of inertia
  69     float thetaDir[MAXCLUSTER];        // cluster angle theta [degree]
  70     float phiDir[MAXCLUSTER];	       // cluster angle phi [degree]
  71     float cosThetaPos[MAXCLUSTER];     // cosine of cluster position theta
  72     float phiPos[MAXCLUSTER];	       // cluster position phi [degree]
  73     float dca[MAXCLUSTER];	       // distance of closest approach to IP [mm]
  74     int   nHitsPerCluster[MAXCLUSTER]; // #hits per cluster
  75     int   trueNcluster;                // true # cluster / event
  76     float trueDist[MAXCLUSTER];        // true distance of MCparticle to IP [mm]
  77     float trueGammaEnergy[MAXCLUSTER]; // true energy of MCparticle gamma [GeV]
  78     float trueGammaDirX[MAXCLUSTER];   // true x direction of MCparticle gamma [mm]
  79     float trueGammaDirY[MAXCLUSTER];   // true y direction of MCparticle gamma [mm]
  80     float trueGammaDirZ[MAXCLUSTER];   // true z direction of MCparticle gamma [mm]
  81     float trueGammaPosX[MAXCLUSTER];   // true x position of MCparticle gamma [mm]
  82     float trueGammaPosY[MAXCLUSTER];   // true y position of MCparticle gamma [mm]
  83     float trueGammaPosZ[MAXCLUSTER];   // true z position of MCparticle gamma [mm]
  84     float trueNeutralinoEnergy[MAXCLUSTER]; // true energy of MCparticle neu1 [GeV]
  85     float trueNeutralinoDirX[MAXCLUSTER];   // true x direction of MCparticle neu1 [mm]
  86     float trueNeutralinoDirY[MAXCLUSTER];   // true y direction of MCparticle neu1 [mm]
  87     float trueNeutralinoDirZ[MAXCLUSTER];   // true z direction of MCparticle neu1 [mm]
  88     float trueThetaDir[MAXCLUSTER];    // true gamma angle theta [degree]
  89     float truePhiDir[MAXCLUSTER];      // true gamma angle phi [degree]
  90   } _clusterFill;
  91 
  92   TFile *outFile = new TFile("/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa_Nhits.root","RECREATE");
  93   TTree *outTree = new TTree("selectedTree","selectedTree");
  94 
  95   outTree->Branch("eventNumber",      &_clusterFill.iEvt,       "eventNumber/I");
  96   outTree->Branch("ECAL_energySum",   &_clusterFill.esumECAL,   "ECAL_energySum/F");
  97   outTree->Branch("ECAL_nHits",       &_clusterFill.nHitsECAL,  "ECAL_nHits/I");
  98   outTree->Branch("HCAL_energySum",   &_clusterFill.esumHCAL,   "HCAL_energySum/F");
  99   outTree->Branch("HCAL_nHits",       &_clusterFill.nHitsHCAL,  "HCAL_nHits/I");
 100   outTree->Branch("calos_energySum",  &_clusterFill.esumCalos,  "calos_energySum/F");
 101   outTree->Branch("calos_nHits",      &_clusterFill.nHitsCalos, "calos_nHits/I");
 102   outTree->Branch("cluster_nCluster", &_clusterFill.nCluster,   "cluster_nCluster/I");
 103   outTree->Branch("cluster_nGamma",   &_clusterFill.nGamma,     "cluster_nGamma/I");
 104   outTree->Branch("cluster_cosThetaPos", &_clusterFill.cosThetaPos,"cluster_cosThetaPos[cluster_nCluster]/F");
 105   outTree->Branch("cluster_dca",      &_clusterFill.dca,        "cluster_dca[cluster_nCluster]/F");
 106   outTree->Branch("cluster_energy",   &_clusterFill.energy,     "cluster_energy[cluster_nCluster]/F");
 107   outTree->Branch("cluster_etMiss",   &_clusterFill.etMiss,     "cluster_etMiss/F");
 108   outTree->Branch("cluster_nHits",    &_clusterFill.nHitsPerCluster, "cluster_nHits[cluster_nCluster]/I");
 109   outTree->Branch("cluster_phiDir",   &_clusterFill.phiDir,     "cluster_phiDir[cluster_nCluster]/F");
 110   outTree->Branch("cluster_phiPos",   &_clusterFill.phiPos,     "cluster_phiPos[cluster_nCluster]/F");
 111   outTree->Branch("cluster_position", &_clusterFill.clusterPos, "cluster_position[cluster_nCluster][3]/F");
 112   outTree->Branch("cluster_clusterAxis", &_clusterFill.clusterAxis,"cluster_axisOfInertia[cluster_nCluster][3]/F");
 113   outTree->Branch("cluster_thetaDir", &_clusterFill.thetaDir,   "cluster_thetaDir[cluster_nCluster]/F");
 114   outTree->Branch("cluster_trueNcluster", &_clusterFill.trueNcluster, "cluster_trueNcluster/I");
 115   outTree->Branch("cluster_trueDist", &_clusterFill.trueDist,   "cluster_trueDist[cluster_nCluster]/F");
 116   outTree->Branch("cluster_trueGammaEnergy", &_clusterFill.trueGammaEnergy, "cluster_trueGammaEnergy[cluster_nCluster]/F");
 117   outTree->Branch("cluster_trueGammaDirX",   &_clusterFill.trueGammaDirX,   "cluster_trueGammaDirX[cluster_nCluster]/F");
 118   outTree->Branch("cluster_trueGammaDirY",   &_clusterFill.trueGammaDirY,   "cluster_trueGammaDirY[cluster_nCluster]/F");
 119   outTree->Branch("cluster_trueGammaDirZ",   &_clusterFill.trueGammaDirZ,   "cluster_trueGammaDirZ[cluster_nCluster]/F");
 120   outTree->Branch("cluster_trueGammaPosX",   &_clusterFill.trueGammaPosX,   "cluster_trueGammaPosX[cluster_nCluster]/F");
 121   outTree->Branch("cluster_trueGammaPosY",   &_clusterFill.trueGammaPosY,   "cluster_trueGammaPosY[cluster_nCluster]/F");
 122   outTree->Branch("cluster_trueGammaPosZ",   &_clusterFill.trueGammaPosZ,   "cluster_trueGammaPosZ[cluster_nCluster]/F");
 123   outTree->Branch("cluster_trueNeutralinoEnergy", &_clusterFill.trueNeutralinoEnergy, "cluster_trueNeutralinoEnergy[cluster_nCluster]/F");
 124   outTree->Branch("cluster_trueNeutralinoDirX",   &_clusterFill.trueNeutralinoDirX,   "cluster_trueNeutralinoDirX[cluster_nCluster]/F");
 125   outTree->Branch("cluster_trueNeutralinoDirY",   &_clusterFill.trueNeutralinoDirY,   "cluster_trueNeutralinoDirY[cluster_nCluster]/F");
 126   outTree->Branch("cluster_trueNeutralinoDirZ",   &_clusterFill.trueNeutralinoDirZ,   "cluster_trueNeutralinoDirZ[cluster_nCluster]/F");
 127   outTree->Branch("cluster_truePhiDir", &_clusterFill.truePhiDir, "cluster_truePhiDir[cluster_nCluster]/F");
 128   outTree->Branch("cluster_trueThetaDir", &_clusterFill.trueThetaDir, "cluster_trueThetaDir[cluster_nCluster]/F");
 129 
 130   for (int fileCounter = 0; fileCounter<nFiles; ++fileCounter) {
 131 
 132     inFile = new TFile(inFileNames[fileCounter],"READ");
 133 
 134     if ( inFile->IsOpen() ) {
 135       
 136       cout<<"opened file "<<fileCounter<<": "<<inFileNames[fileCounter]<<endl;
 137       
 138       /******** read input tree **********/
 139       int   _iEvt;
 140       int   _nCluster;  int   _nGamma;    int _trueNcluster;
 141       float _etMiss;
 142       float _esumCalos; float _esumECAL;  float _esumHCAL; 
 143       int _nHitsCalos;  int   _nHitsECAL; int   _nHitsHCAL;
 144       float _energy[MAXCLUSTER]; int _nHitsPerCluster[MAXCLUSTER];
 145       float _clusterPos[MAXCLUSTER][3];   float _clusterAxis[MAXCLUSTER][3]; 
 146       float _thetaDir[MAXCLUSTER]; float _cosThetaPos[MAXCLUSTER];
 147       float _phiDir[MAXCLUSTER];   float _phiPos[MAXCLUSTER];
 148       float _dca[MAXCLUSTER];
 149       float _trueDist[MAXCLUSTER]; 
 150       float _trueGammaEnergy[MAXCLUSTER]; float _trueNeutralinoEnergy[MAXCLUSTER];
 151       float _trueGammaDirX[MAXCLUSTER]; float _trueGammaDirY[MAXCLUSTER]; float _trueGammaDirZ[MAXCLUSTER];
 152       float _trueGammaPosX[MAXCLUSTER]; float _trueGammaPosY[MAXCLUSTER]; float _trueGammaPosZ[MAXCLUSTER];
 153       float _trueNeutralinoDirX[MAXCLUSTER]; float _trueNeutralinoDirY[MAXCLUSTER]; float _trueNeutralinoDirZ[MAXCLUSTER];
 154       float _trueThetaDir[MAXCLUSTER]; float _truePhiDir[MAXCLUSTER];
 155 
 156       TTree* inTree = (TTree*)inFile->Get("tree");
 157       inTree->SetBranchAddress("eventNumber",     &_iEvt);
 158       inTree->SetBranchAddress("ECAL_energySum",  &_esumECAL);
 159       inTree->SetBranchAddress("ECAL_nHits",      &_nHitsECAL);
 160       inTree->SetBranchAddress("HCAL_energySum",  &_esumHCAL);
 161       inTree->SetBranchAddress("HCAL_nHits",      &_nHitsHCAL);
 162       inTree->SetBranchAddress("calos_energySum", &_esumCalos);
 163       inTree->SetBranchAddress("calos_nHits",     &_nHitsCalos);
 164       inTree->SetBranchAddress("cluster_cosThetaPos",  &_cosThetaPos);
 165       inTree->SetBranchAddress("cluster_dca",          &_dca);
 166       inTree->SetBranchAddress("cluster_energy",       &_energy);
 167       inTree->SetBranchAddress("cluster_etMiss",       &_etMiss);
 168       inTree->SetBranchAddress("cluster_nCluster",     &_nCluster);
 169       inTree->SetBranchAddress("cluster_nGamma",       &_nGamma);
 170       inTree->SetBranchAddress("cluster_nHits",        &_nHitsPerCluster);
 171       inTree->SetBranchAddress("cluster_phiDir",       &_phiDir);
 172       inTree->SetBranchAddress("cluster_phiPos",       &_phiPos);
 173       inTree->SetBranchAddress("cluster_position",     &_clusterPos);
 174       inTree->SetBranchAddress("cluster_clusterAxis",  &_clusterAxis);
 175       inTree->SetBranchAddress("cluster_thetaDir",     &_thetaDir);
 176       inTree->SetBranchAddress("cluster_trueNcluster", &_trueNcluster);
 177       inTree->SetBranchAddress("cluster_trueDist",     &_trueDist);
 178       inTree->SetBranchAddress("cluster_trueGammaEnergy", &_trueGammaEnergy);
 179       inTree->SetBranchAddress("cluster_trueGammaDirX",   &_trueGammaDirX);
 180       inTree->SetBranchAddress("cluster_trueGammaDirY",   &_trueGammaDirY);
 181       inTree->SetBranchAddress("cluster_trueGammaDirZ",   &_trueGammaDirZ);
 182       inTree->SetBranchAddress("cluster_trueGammaPosX",   &_trueGammaPosX);
 183       inTree->SetBranchAddress("cluster_trueGammaPosY",   &_trueGammaPosY);
 184       inTree->SetBranchAddress("cluster_trueGammaPosZ",   &_trueGammaPosZ);
 185       inTree->SetBranchAddress("cluster_trueNeutralinoEnergy", &_trueNeutralinoEnergy);
 186       inTree->SetBranchAddress("cluster_trueNeutralinoDirX",   &_trueNeutralinoDirX);
 187       inTree->SetBranchAddress("cluster_trueNeutralinoDirY",   &_trueNeutralinoDirY);
 188       inTree->SetBranchAddress("cluster_trueNeutralinoDirZ",   &_trueNeutralinoDirZ);
 189       inTree->SetBranchAddress("cluster_truePhiDir",   &_truePhiDir);
 190       inTree->SetBranchAddress("cluster_trueThetaDir", &_trueThetaDir);
 191 
 192       int nEntries = inTree->GetEntriesFast();  
 193       cout<<" found input tree with "<<nEntries<<" entries"<<endl;
 194 
 195       //weight events according to expected crosssection
 196       //      polWeight[fileCounter] = 1;
 197       weight[fileCounter] = nExpected[fileCounter]*polWeight[fileCounter]/nEntries;
 198       cout<<" and weight "<<weight[fileCounter]<<endl;
 199 
 200       /******** apply selection criteria chosen above **********/
 201       int selectedEvents = 0;
 202       int notEcluster    = 0;
 203       int notCosTheta    = 0;
 204       int notNhits       = 0;
 205       int notEtmiss      = 0;
 206       int notEsum        = 0;
 207       int notNgamma      = 0;
 208       int isGood[nEntries];
 209       int isCosTheta[nEntries];
 210 
 211       for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
 212 	inTree->GetEntry(entryCounter);
 213 	isGood[entryCounter] = 0;
 214 	isCosTheta[entryCounter] = 0;
 215 	if ( _esumECAL > _EminEcal && _esumECAL < _EmaxEcal) {
 216 	  if ( _etMiss > _minEt_miss) { 
 217 	    if ( _nHitsECAL > _minNhitsEcal && _nHitsECAL < _maxNhitsEcal) {
 218 	      if ( _nCluster >= _minNgamma ) {
 219 		for (int clusterCounter=0; clusterCounter<_nCluster; clusterCounter++) {	      
 220 		  if ( (_cosThetaPos[clusterCounter] < _maxCosTheta) && (_cosThetaPos[clusterCounter] > -1*_maxCosTheta) ) {
 221 		    isCosTheta[entryCounter]++;
 222 		    if ( _energy[clusterCounter] > _minEcluster ) {
 223 		      isGood[entryCounter]++;
 224 		    }
 225 		    else notEcluster++; 
 226 		  }
 227 		  else  notCosTheta++; 
 228 		}//all cluster	  
 229 	      }
 230 	      else notNgamma++;
 231 	    }
 232 	    else notNhits++; 
 233 	  }
 234 	  else notEtmiss++; 
 235 	}
 236 	else notEsum++; 
 237 	
 238       }//all entries
 239       cout<<" looped over input tree "<<endl;
 240       
 241       /******** fill output tree **********/
 242       for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
 243 	inTree->GetEntry(entryCounter);
 244 
 245 	if ( isGood[entryCounter] >= _minNgamma && isGood[entryCounter] > 0 && isCosTheta[entryCounter] == _nCluster ) {
 246 
 247 	  //	  cout<<isCosTheta[entryCounter]<<" of "<<_nCluster<<" cluster with |cos(theta)| > "<<_maxCosTheta<<endl;      
 248 
 249 	  _clusterFill.iEvt         = _iEvt;
 250 	  _clusterFill.esumECAL     = _esumECAL;
 251 	  _clusterFill.esumHCAL     = _esumHCAL;
 252 	  _clusterFill.esumCalos    = _esumCalos;
 253 	  _clusterFill.nHitsECAL    = _nHitsECAL;
 254 	  _clusterFill.nHitsHCAL    = _nHitsHCAL;
 255 	  _clusterFill.nHitsCalos   = _nHitsCalos;
 256 	  _clusterFill.etMiss       = _etMiss;
 257 	  _clusterFill.nCluster     = _nCluster;
 258 	  _clusterFill.nGamma       = _nGamma;
 259 	  _clusterFill.trueNcluster = _trueNcluster;
 260 
 261 	  for (int clusterCounter=0; clusterCounter<_nCluster; ++clusterCounter) {	      	  	    
 262 	    _clusterFill.energy[clusterCounter]          = _energy[clusterCounter];
 263 	    _clusterFill.nHitsPerCluster[clusterCounter] = _nHitsPerCluster[clusterCounter];
 264 	    _clusterFill.clusterPos[clusterCounter][0]   = _clusterPos[clusterCounter][0];
 265 	    _clusterFill.clusterPos[clusterCounter][1]   = _clusterPos[clusterCounter][1];
 266 	    _clusterFill.clusterPos[clusterCounter][2]   = _clusterPos[clusterCounter][2];
 267 	    _clusterFill.clusterAxis[clusterCounter][0]  = _clusterAxis[clusterCounter][0];
 268 	    _clusterFill.clusterAxis[clusterCounter][1]  = _clusterAxis[clusterCounter][1];
 269 	    _clusterFill.clusterAxis[clusterCounter][2]  = _clusterAxis[clusterCounter][2];
 270 	    _clusterFill.thetaDir[clusterCounter]        = _thetaDir[clusterCounter];
 271 	    _clusterFill.phiDir[clusterCounter]          = _phiDir[clusterCounter];
 272 	    _clusterFill.cosThetaPos[clusterCounter]     = _cosThetaPos[clusterCounter];
 273 	    _clusterFill.phiPos[clusterCounter]          = _phiPos[clusterCounter];
 274 	    _clusterFill.dca[clusterCounter]             = _dca[clusterCounter];
 275 	    _clusterFill.trueDist[clusterCounter]        = _trueDist[clusterCounter];
 276 	    _clusterFill.trueGammaEnergy[clusterCounter] = _trueGammaEnergy[clusterCounter];
 277 	    _clusterFill.trueGammaDirX[clusterCounter]   = _trueGammaDirX[clusterCounter];
 278 	    _clusterFill.trueGammaDirY[clusterCounter]   = _trueGammaDirY[clusterCounter];
 279 	    _clusterFill.trueGammaDirZ[clusterCounter]   = _trueGammaDirZ[clusterCounter];
 280 	    _clusterFill.trueGammaPosX[clusterCounter]   = _trueGammaPosX[clusterCounter];
 281 	    _clusterFill.trueGammaPosY[clusterCounter]   = _trueGammaPosY[clusterCounter];
 282 	    _clusterFill.trueGammaPosZ[clusterCounter]   = _trueGammaPosZ[clusterCounter];
 283 	    _clusterFill.trueNeutralinoEnergy[clusterCounter] = _trueNeutralinoEnergy[clusterCounter];
 284 	    _clusterFill.trueNeutralinoDirX[clusterCounter]   = _trueNeutralinoDirX[clusterCounter];
 285 	    _clusterFill.trueNeutralinoDirY[clusterCounter]   = _trueNeutralinoDirY[clusterCounter];
 286 	    _clusterFill.trueNeutralinoDirZ[clusterCounter]   = _trueNeutralinoDirZ[clusterCounter];
 287 	    _clusterFill.trueThetaDir[clusterCounter]    = _trueThetaDir[clusterCounter];
 288 	    _clusterFill.truePhiDir[clusterCounter]      = _truePhiDir[clusterCounter];
 289 	  }//all cluster
 290 
 291 	  outTree->SetWeight(weight[fileCounter]);
 292 	  outTree->Fill();
 293 	  selectedEvents++;    		    
 294 	}//is good
 295       }//all entries
 296       
 297       cout<<"writing output tree with "<<selectedEvents<<" selected events"<<endl;
 298       cout<<"\t found "<<notEcluster<<" cluster with less then "<<_minEcluster<<" GeV per cluster"<<endl;
 299       cout<<"\t found "<<notCosTheta<<" cluster with cos(theta) more then "<<_maxCosTheta<<endl;
 300       cout<<"\t found "<<notNhits<<" events with less then "<<_minNhitsEcal<<" hits in ECAL"<<endl;
 301       cout<<"\t found "<<notEtmiss<<" events with less then "<<_minEt_miss<<" GeV transverse energy"<<endl;
 302       cout<<"\t found "<<notEsum<<" events with less then "<<_EminEcal<<" GeV in ECAL"<<endl;
 303       cout<<"\t found "<<notNgamma<<" events with less then "<<_minNgamma<<" photons"<<endl;
 304 
 305     }//file is open
 306     inFile->Close();
 307   }//all files
 308   
 309   outFile->Write();
 310 
 311 }//eventSelection
 312 

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.