StopMuRange.cxx File Reference

#include <iostream>
#include <fstream>
#include "TROOT.h"
#include "TClonesArray.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TF1.h"
#include "TString.h"
#include "CalDetTracker/CDPIDInfo.h"
#include "CalDetTracker/CDTrackInfo.h"
#include "CalDetTracker/CDTrackedHitInfo.h"

Go to the source code of this file.

Functions

TROOT simple ("simple","simple", 0)
int main (int argc, char **argv)

Function Documentation

int main ( int  argc,
char **  argv 
)

COMMAND LINE

CALIB CONSTANTS

BEAM RUN KEY

Cerenkov ADC cuts:

ATTENUATION PARAMS

FILE + TREE

Definition at line 20 of file StopMuRange.cxx.

References MuELoss::e, exit(), CDTrackInfo::GetAngle(), CDTrackedHitInfo::GetCharge(), CDTrackedHitInfo::GetEnd(), CDPIDInfo::GetKovADC3(), CDTrackInfo::GetNHits(), CDTrackedHitInfo::GetPlane(), CDTrackInfo::GetRange(), CDTrackInfo::GetResult(), CDTrackedHitInfo::GetStrip(), CDTrackedHitInfo::GetTime(), CDTrackedHitInfo::GetTransPos(), CDTrackInfo::GetTriggerTime(), and CDTrackInfo::GetVertex().

00020                                {
00021   
00022   int beamRunKey = 0;
00023   char *ccTag;
00024   char cctag[256];
00025   int usePID = 0;  //false, don't use PID by default
00026 
00028   //read in command line options or ask for some options
00029 
00030   if(argc==3) {
00031     beamRunKey = atoi(argv[1]);
00032     ccTag = argv[2];
00033   }
00034   else if(argc==4) {
00035     beamRunKey = atoi(argv[1]);
00036     ccTag = argv[2];
00037     usePID = atoi(argv[3]);
00038   }
00039   else {
00040     char charbeamrunkey[30];
00041     cout << "Enter beamRunKey" << endl;
00042     cin >> charbeamrunkey;
00043     beamRunKey = atoi(charbeamrunkey);
00044     
00045     cout << "Enter CalibConstants file" << endl;
00046     cin >> cctag;
00047     ccTag = cctag;
00048 
00049     cout << "Use PID? 0=no, 1=yes" << endl;
00050     cin >> usePID;
00051   }
00052 
00053   //some constants:
00054   double den_iron = 7.87;   // g/cm**3
00055   double den_scint = 1.032; // g/cm**3
00056 
00057 
00059   //Read in calibration constants
00060   float theMeans[120][24] = {0};
00061   float temp_nent = 0;
00062   int temp_plane = 0;
00063   int temp_strip = 0;
00064   float temp_mn = 0;
00065   float temp_rms = 0;
00066   float temp_mean = 0;
00067   float temp_meanerr = 0;
00068   float temp_non = 0;
00069   ifstream inmean;
00070   inmean.open(ccTag);
00071   if(inmean) cout << "Opened CalibConstants file: " << ccTag << endl;
00072   if(!inmean) {
00073     cerr << "File not open: check that file exists!." << endl;
00074     exit(0);
00075   }
00076   double totalMean[4]={0};
00077   double aveMean[4]={0};
00078   double totalNent[4]={0};
00079   char isItANan[20];
00080   while(inmean>>temp_plane>>temp_strip>>temp_nent) {
00081     if(temp_nent<2) {
00082       inmean >> isItANan >> isItANan >> isItANan;
00083       inmean >> isItANan >> isItANan >> isItANan;
00084       cout << "Missing: " << temp_plane << "\t" << temp_strip << endl;
00085       theMeans[temp_plane][temp_strip] =-1;
00086     }
00087     else {
00088       inmean>>temp_non>>temp_non>>temp_mean>>temp_meanerr>>temp_mn>>temp_rms;
00089       theMeans[temp_plane][temp_strip] = temp_mean;
00090       int index=(temp_plane%2)+2*(temp_plane>=60);
00091       totalMean[index]+=temp_mean;
00092       totalNent[index]++;
00093     }
00094   }
00095   inmean.close();
00096   
00097   for(int i=0;i<4;i++) {
00098     if(totalNent[i]>0) {
00099       aveMean[i]=totalMean[i]/totalNent[i];
00100     }
00101     else {
00102       aveMean[i]=1.;
00103     }
00104   }
00105   
00106   for(int pl=0;pl<120;pl++) {
00107     for(int st=0;st<24;st++) {
00108       int index=(pl%2)+2*(pl>=60);
00109       if(theMeans[pl][st]<=0) {
00110         theMeans[pl][st]=aveMean[index];
00111       }
00112     }
00113   }
00114   cout << "Calib constants read in" << endl;
00115 
00116 
00118   //read in file and read in appropriate info
00119   ifstream BeamRunKeyFile("beamRunKey.dat");
00120   if(!BeamRunKeyFile) {
00121     cerr << "Couldn't open file: beamRunKey.dat" << endl;
00122     exit(0);
00123   }
00124   else cout << "Opened file: beamRunKey.dat" << endl;
00125 
00126   char runDir[256];
00127   int cnt1 = 0;
00128   while(cnt1<beamRunKey){
00129     BeamRunKeyFile>>runDir;
00130     cnt1+=1;
00131   }
00132   cout << "Using directory " << runDir 
00133        << " for beamRunKey=" << beamRunKey <<endl;
00134 
00135 
00137   //need beam energy in order to set stopping range upper and lower limits
00138   //use energies: 0.8,1.0,1.2,1.4,1.6,1.8,2.0
00139   int numberOfEnergies = 7;
00140   double *beamEnergy = new double[numberOfEnergies];
00141   beamEnergy[0] = 0.8;
00142   beamEnergy[1] = 1.0;
00143   beamEnergy[2] = 1.2;
00144   beamEnergy[3] = 1.4;
00145   beamEnergy[4] = 1.6;
00146   beamEnergy[5] = 1.8;
00147   beamEnergy[6] = 2.0;
00148   int *beamEnergyFactor10 = new int[numberOfEnergies];
00149   float *lowPlaneRangeCut = new float[numberOfEnergies];
00150   float *highPlaneRangeCut = new float[numberOfEnergies];
00151   for(int i=0;i<numberOfEnergies;i++) {
00152     beamEnergyFactor10[i] = int(10.0*(beamEnergy[i]+0.01));
00153     lowPlaneRangeCut[i] = 0;
00154     highPlaneRangeCut[i] = 0;
00155   }
00156   
00157   {
00158     ifstream BeamEnergyAndRange("theMuonRanges.dat");
00159     if(!BeamEnergyAndRange) {
00160       cerr << "Couldn't open file: theMuonRanges.dat" << endl;
00161       exit(0);
00162     }
00163     else cout << "Opened file: theMuonRanges.dat" << endl;
00164     
00165     float tempEnergy, lowPlane, highPlane;
00166     while(BeamEnergyAndRange >> tempEnergy >> lowPlane >> highPlane) {
00167       highPlane-=1;
00168       int tempIntEnergy=int(10.0*(tempEnergy+0.01));
00169       for(int i=0;i<numberOfEnergies;i++){
00170         if(tempIntEnergy==beamEnergyFactor10[i]){ 
00171           lowPlaneRangeCut[i]= lowPlane;
00172           highPlaneRangeCut[i] = highPlane;
00173           break;
00174         }
00175       }
00176     }
00177   }
00178   
00179   for(int i=0;i<numberOfEnergies;i++)
00180     cout << "For " << beamEnergy[i] << " GeV, using range cuts: " 
00181          << lowPlaneRangeCut[i] << " to " << highPlaneRangeCut[i] << endl;
00182 
00183 
00185   // Need cerenkov ADC upper and lower values for cutting on muons
00186   // This is only applicable to 1.8, 2.0 GeV runs since T7,T11 Cerenkovs 
00187   // could not be pumped high enough to see muons at lower energies. 
00188   // At higher energies, muons do not stop in detector.
00189   float *lowADCKovCut = new float[numberOfEnergies];
00190   float *highADCKovCut = new float[numberOfEnergies];
00191   lowADCKovCut[0] = 0.; highADCKovCut[0] = 100.;
00192   lowADCKovCut[1] = 0.; highADCKovCut[1] = 100.;
00193   lowADCKovCut[2] = 0.; highADCKovCut[2] = 100.;
00194   lowADCKovCut[3] = 0.; highADCKovCut[3] = 100.;
00195   lowADCKovCut[4] = 0.; highADCKovCut[4] = 100.;
00196   lowADCKovCut[5] = 100.; highADCKovCut[5] = 1000.;
00197   lowADCKovCut[6] = 600.; highADCKovCut[6] = 1500.;
00198 
00199 
00201   //read in params for attenuation correction
00202   float attenpars[4][24][9] = {};
00203 
00204   ifstream inatpar;
00205   char attname[256];
00206   sprintf(attname,"StripAttenParsNewest.dat");
00207   inatpar.open(attname);
00208   if(!inatpar) {
00209     cerr << "Couldn't open file: StripAttenParsNewest.dat" << endl;
00210     exit(0);
00211   }
00212   else cout << "Opened file: StripAttenParsNewest.dat" << endl;
00213 
00214   int att_crt=0,att_pln=0,att_stp=0;
00215   while(inatpar>>att_crt>>att_pln>>att_stp){    
00216     for(int j=0;j<9;j++){
00217       float tempro = 0;
00218         inatpar>>tempro;
00219         attenpars[(att_crt*2)+att_pln][att_stp][j] = tempro;
00220     }
00221   }
00222   
00223   TF1 *form[4][24];
00224   char formname[256];
00225   for(int i=0;i<4;i++){
00226     for(int j=0;j<24;j++){
00227       sprintf(formname,"form%i_%i",i,j);
00228       form[i][j] = new TF1(formname,"pol8",-0.5,23.5);
00229       for(int k=0;k<9;k++) form[i][j]->SetParameter(k,attenpars[i][j][k]);
00230     }
00231   }
00232 
00233 
00235   //Get tag out of calibconstants filename:
00236   TString ccStrg(ccTag);
00237   ccStrg.Chop(); // get
00238   ccStrg.Chop(); // rid
00239   ccStrg.Chop(); // of
00240   ccStrg.Chop(); // ".dat"
00241   //get rid of everything before 1st "_"
00242   if(ccStrg.Index("_")>0) ccStrg.Remove(0,ccStrg.Index("_")+1); 
00243   cout << "Calib constants filename tag: " << ccStrg.Data() << endl;
00244 
00245   //open output file
00246   char savename[256];
00247   if(usePID) sprintf(savename,"StopMu_RunKey%i_%s_PID.root",
00248                      beamRunKey,ccStrg.Data());
00249   else sprintf(savename,"StopMu_RunKey%i_%s.root",beamRunKey,ccStrg.Data());
00250   cout << "Opening output file: " << savename << endl;
00251   TFile *save = new TFile(savename,"RECREATE");
00252   
00253   cout << "Setting up output tree" << endl;
00254   //range tree 
00255   Double_t treeRange,treeTotPE,treeTotMIP,treeDzDs,treeBeamEnergy,treeTrkMom;
00256   int treeLTP;
00257   TTree *rangeTree = new TTree("rangeTree","Tree of Muon Ranges");
00258   rangeTree->Branch("beamEnergy",&treeBeamEnergy,"beamEnergy/D");
00259   rangeTree->Branch("range",&treeRange,"range/D");
00260   rangeTree->Branch("momentum",&treeTrkMom,"momentum/D");
00261   rangeTree->Branch("totPE",&treeTotPE,"totPE/D");
00262   rangeTree->Branch("totMIP",&treeTotMIP,"totMIP/D");
00263   rangeTree->Branch("dzds",&treeDzDs,"dzds/D");
00264   rangeTree->Branch("LTP",&treeLTP,"LTP/I");
00265 
00266   cout << "Starting to process files" << endl;
00267   
00268   int start_i = 0;
00269   if(usePID) start_i = 5; //if we are using PID then only look at 1.8,2.0GeV
00270 
00271   for(int i=start_i;i<numberOfEnergies;i++){
00272     for(int j=0;j<2;j++){
00273 
00274       TChain *tree = new TChain("TrackerTree");
00275       char theDir[256];
00276       if(j==0) {
00277         sprintf(theDir,"%s/negatives/%1.1f/*",runDir,beamEnergy[i]);
00278         treeBeamEnergy = -beamEnergy[i];
00279       }
00280       else {
00281         sprintf(theDir,"%s/positives/%1.1f/*",runDir,beamEnergy[i]);
00282         treeBeamEnergy = beamEnergy[i];
00283       }
00284 
00285       int numfiles = tree->Add(theDir);
00286       if(numfiles==0) {
00287         cout << "No files for " << treeBeamEnergy << " GeV" << endl;
00288         continue;
00289       }
00290 
00291       int ent = int(tree->GetEntries());
00292       int ind = 0;
00293       
00294       TClonesArray *clonearray = new TClonesArray("CDTrackedHitInfo",1000);
00295       CDTrackInfo *cdti = new CDTrackInfo();
00296       CDPIDInfo *pidi = new CDPIDInfo();
00297       
00298       tree->SetBranchAddress("TrackedHitInfo",&clonearray);
00299       tree->SetBranchAddress("TrackInfo",&cdti);
00300       if(usePID) tree->SetBranchAddress("PIDInfo",&pidi);      
00301       
00302       cout << "Beam energy " << treeBeamEnergy 
00303            << " GeV... Processing " << numfiles << " Files, ("  
00304            <<  tree->GetEntries() << " events) :"<< endl; 
00305       
00306       int timeErrCnt = 0; //counts number of events with out of time hits
00307 
00308       while(ind<ent){
00309         
00310         if(ind%10000==0) cout<<double(ind)/double(ent)*100.0<<"%" << endl;
00311         
00312         tree->GetEvent(ind);
00313         float *evenang = cdti->GetAngle(0);
00314         float *oddang = cdti->GetAngle(1);
00315         float *vertex1 = cdti->GetVertex(0);
00316         float *vertex2 = cdti->GetVertex(1);
00317         float trkRange =cdti->GetRange();
00318         float nHitsEven = cdti->GetNHits(0);
00319         float nHitsOdd = cdti->GetNHits(1);
00320         
00321         treeLTP = int(vertex2[0]);
00322 
00323         if(cdti->GetResult(0)==1&&cdti->GetResult(1)==1//tracked in both views
00324            &&fabs(evenang[1])<5&&fabs(oddang[1])<5   //small angle to beamline
00325            &&(vertex1[0]==0||vertex1[0]==1)        //first hit plane is 0 or 1
00326            &&!(0.8*trkRange>(nHitsEven+nHitsOdd))//#hits consistent with range
00327            && clonearray){
00328 
00329           TClonesArray &c = *clonearray;
00330           int clonesize = clonearray->GetEntries();
00331           int ind2 = 0;
00332           float mips[120][24] = {0.0};
00333           float pes[120][24] = {0.0};
00334           float adcs[120][24] = {0.0};
00335           float numhits = 0;
00336           
00337           bool useEvent = true;
00338           
00339           bool fillRangeTree =false;
00340           if(usePID) {
00341             if(pidi->GetKovADC3()>lowADCKovCut[i]
00342                &&pidi->GetKovADC3()<highADCKovCut[i]&&trkRange<62.0) 
00343               fillRangeTree = true;
00344           }
00345           else if(trkRange>=lowPlaneRangeCut[i]&&trkRange<=highPlaneRangeCut[i]
00346                   &&trkRange<62.0) {//latter indicates a tracking problem
00347             fillRangeTree = true;
00348           }
00349           else useEvent=false;
00350           
00351           
00352           treeTotMIP=0;
00353           treeTotPE=0;
00354           treeRange=trkRange;
00355           float range_gmm=(trkRange-1)*den_iron*2.54 + trkRange*den_scint*1.0;
00356           treeTrkMom = 0.105658*exp(1.0363*log(range_gmm)-4.383);
00357 
00358           treeDzDs = 1.0/sqrt(1.0 + pow(TMath::Tan(evenang[1]
00359                                                    *TMath::Pi()/180.0),2)
00360                               + pow(TMath::Tan(oddang[1]
00361                                                *TMath::Pi()/180.0),2));
00362           
00363           while(ind2<clonesize) {
00364             CDTrackedHitInfo *cdthi = (CDTrackedHitInfo*) c[ind2];
00365             ind2++;
00366             
00367             numhits+=1;
00368             
00369             int end = cdthi->GetEnd();
00370             int pln = cdthi->GetPlane();
00371             int stp = cdthi->GetStrip();
00372             
00373             if(pln<20&&(stp<8||stp>16)) useEvent=false; //use only events that
00374                                                         //look like beam muons
00375 
00376             //make sure time is consistent with trigger time for all hits
00377             if((cdthi->GetTime())/1.5625e-9 
00378                - cdti->GetTriggerTime()>-200 
00379                && (cdthi->GetTime())/1.5625e-9
00380                - cdti->GetTriggerTime()<100
00381                && useEvent==true) useEvent=true;
00382             else {
00383               if(useEvent) 
00384                 timeErrCnt+=1;
00385               useEvent=false;
00386             }
00387            
00388             if(!useEvent) break;
00389  
00390             //prepare for doing an attenuation correction
00391             int ref = 2*(end-1) + pln%2;
00392             float attencor = form[ref][stp]->Eval(11.5)
00393               /form[ref][stp]->Eval(cdthi->GetTransPos());
00394             attencor=1;
00395             
00396             float ds=1;
00397             float mip = cdthi->GetCharge(1)/theMeans[60*(end-1)+pln][stp];
00398             float pe = cdthi->GetCharge(3);
00399             float adc = cdthi->GetCharge(0);
00400             
00401             mip*=attencor;
00402             mip/=ds;
00403             pe*=attencor;
00404             pe/=ds;
00405             adc*=attencor;
00406             adc/=ds;
00407         
00408             treeTotPE+=pe;
00409             treeTotMIP+=mip;
00410             
00411             mips[60*(end-1)+pln][stp]=mip;
00412             pes[60*(end-1)+pln][stp]=pe;
00413             adcs[60*(end-1)+pln][stp]=adc;
00414             
00415           }
00416 
00417           if(useEvent){
00418             if(fillRangeTree) rangeTree->Fill();
00419           }
00420         }       
00421         ind++;
00422       }
00423 
00424       delete clonearray;
00425       delete cdti;
00426       delete pidi;
00427       
00428       cout << "Number of events with timing errors: " << timeErrCnt << endl;
00429 
00430     }
00431   }
00432  
00433   save->cd();
00434   save->Write();
00435   save->Close();
00436   
00437   return 0;
00438   
00439 }

TROOT simple ( "simple"  ,
"simple"  ,
 
)

Generated on 8 Jul 2019 for loon by  doxygen 1.6.1