PSMuCal.cxx File Reference

#include <fstream>
#include <iostream>
#include "TROOT.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "PathLengthFuncs.C"
#include "CDTrackInfo.h"
#include "CDXTalkHitInfo.h"
#include "CDTrackedHitInfo.h"
#include "CDTrackerOptions.h"

Go to the source code of this file.

Functions

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

Function Documentation

int main ( int  argv,
char **  argc 
)

AttenPars

Output file

Definition at line 21 of file PSMuCal.cxx.

References CDTrackInfo::GetAngle(), CDTrackedHitInfo::GetCharge(), CDTrackedHitInfo::GetEnd(), CDTrackedHitInfo::GetPlane(), CDTrackInfo::GetRange(), CDTrackInfo::GetResult(), CDTrackedHitInfo::GetStrip(), CDTrackedHitInfo::GetTransPos(), and CDTrackInfo::GetVertex().

00021                                {
00022 
00023   int NumFiles = argv-2;
00024   char fileID[256];
00025   char runnum[100][256];
00026   
00027   if(NumFiles<=0) {
00028     if(NumFiles==-1) {
00029       std::cout << "Enter tag for output file" << std::endl;
00030       std::cin >> fileID;
00031     }
00032     else sprintf(fileID,"%s",argc[1]);
00033     
00034     std::cout << "Enter Number of Runs to process" << std::endl;      
00035     std::cin >> NumFiles;
00036     std::cout << "Enter each file name and path separated by enter" 
00037               << std::endl;
00038     for(int i=0;i<NumFiles;i++) std::cin >> runnum[i];
00039     
00040   }
00041   else {
00042     sprintf(fileID,"%s",argc[1]);
00043     for(int i=0;i<NumFiles;i++) sprintf(runnum[i],"%s",argc[i+2]);
00044   }
00045   
00046   //gains
00047   float theGains[120][24] = {};
00048   
00050   float attenpars[4][24][9] = {}; 
00051   ifstream inatpar;
00052   inatpar.open("StripAttenParsNewest.dat");
00053   if(inatpar) std::cout << "Opened Attenuation Params" <<std::endl;
00054   int att_crt=0,att_pln=0,att_stp=0;
00055   while(inatpar>>att_crt>>att_pln>>att_stp){    
00056     for(int j=0;j<9;j++){
00057       float tempro = 0;
00058       inatpar>>tempro;
00059       attenpars[(att_crt*2)+att_pln][att_stp][j] = tempro;
00060     }
00061   }
00062 
00063   //Make a bunch of TF1's for doing the attenuation corrections:
00064   TF1 *form[4][24];
00065   char formname[256];
00066   for(int i=0;i<4;i++){
00067     for(int j=0;j<24;j++){
00068       sprintf(formname,"form%i_%i",i,j);
00069       form[i][j] = new TF1(formname,"pol8",-0.5,23.5);
00070       for(int k=0;k<9;k++) form[i][j]->SetParameter(k,attenpars[i][j][k]);
00071     }
00072   }
00073 
00075   char savename[256];
00076   sprintf(savename,"PSMuCal%s.root",fileID);
00077   TFile *save = new TFile(savename,"RECREATE");
00078 
00079   //Strip-End Histos:
00080   TH1F *hist[120][24];
00081   for(int i=0;i<120;i++){
00082     for(int j=0;j<24;j++){
00083       char name[256];
00084       sprintf(name,"stripend%i_%i",i,j);
00085       hist[i][j] = new TH1F(name,name,5000,-0.5,4999.5);
00086     }
00087   }
00088 
00089   //Put all the files into a TChain:
00090   cout << "###########" << endl;
00091   cout << "Adding files to TrackerTree TChain" << endl;
00092   TChain *tree = new TChain("TrackerTree");
00093   for(int i=0;i<NumFiles;i++) tree->Add(runnum[i]);
00094   cout << "Total of " << tree->GetEntries() << " entries" 
00095             << endl;
00096 
00097   //Start to make the strip-end histos for the strip-to-strip calib:
00098   cout << "###########" << endl;
00099   cout << "Making Strip-End Histograms" << endl;
00100 
00101   int ind = 0;
00102   int ent = int(tree->GetEntries());
00103   
00104   //To read the track info and hit info from the tree:
00105   TClonesArray *clonearray = new TClonesArray("CDTrackedHitInfo",1000);
00106   CDTrackInfo *cdti = new CDTrackInfo();
00107   tree->SetBranchAddress("StraightTrackedHitInfo",&clonearray);
00108   tree->SetBranchAddress("TrackInfo",&cdti);
00109 
00110   //Start loop over all events:
00111   while(ind<ent){
00112     
00113     if(ind%50000==0) 
00114       cout << ind << " " << double(ind)/double(ent)*100.0 <<  "%" << endl;
00115     
00116     tree->GetEvent(ind);
00117     float *evenang = cdti->GetAngle(0);
00118     float *oddang = cdti->GetAngle(1);
00119     float *vertex1 = cdti->GetVertex(0);
00120     float *vertex2 = cdti->GetVertex(1);
00121     
00122     //select out tracks with vertical angle > +/- 10deg:
00123     if(cdti->GetResult(0)==1 && cdti->GetResult(1)==1
00124        && fabs(evenang[1])<10. //just use beam/PS 
00125        && oddang[1]<-10. && oddang[1]>-50. //T11
00126        //&& fabs(evenang[1])<10. && fabs(oddang[1])<10. //T7
00127        && cdti->GetRange()>=10. //10 T11, 60 T7
00128        && clonearray
00129        ){
00130       
00131       TClonesArray &c = *clonearray;
00132       int clonesize = clonearray->GetEntries();
00133       int ind2 = 0;
00134       float charge[120][24] = {};
00135       float transpos[120][24] = {};
00136       float numhits = 0;
00137       
00138       //loop over all hits in track in event:
00139       while(ind2<clonesize) {
00140         
00141         CDTrackedHitInfo *cdthi = (CDTrackedHitInfo*) c[ind2];
00142         ind2++;
00143         
00144         numhits+=1;
00145         
00146         int end = cdthi->GetEnd();
00147         int pln = cdthi->GetPlane();
00148         int stp = cdthi->GetStrip();
00149 
00150         theGains[60*(end-1)+pln][stp]=cdthi->GetCharge(0)/cdthi->GetCharge(3);
00151         
00152         int ref = 2*(end-1) + pln%2;
00153         float attencor = form[ref][stp]->Eval(11.5)
00154           /form[ref][stp]->Eval(cdthi->GetTransPos());
00155         
00156         float coradc = cdthi->GetCharge(1);
00157         
00158         charge[60*(end-1)+pln][stp] = coradc*attencor;
00159         transpos[60*(end-1)+pln][stp] = cdthi->GetTransPos();
00160       }
00161       
00162       //loop over stored hits to determine single- and double-ended hits,etc:
00163       for(int loop1=0;loop1<60;loop1++){
00164         
00165         int num0 = 0;
00166         int num1 = 0;
00167         float f0 = 0;
00168         float f1 = 0;
00169         float thetranspos0 = -1;
00170         float thetranspos1 = -1;
00171         
00172         for(int loop2=0;loop2<24;loop2++){
00173           //fill strip-end histos:
00174           f0=charge[loop1][loop2];
00175           f1=charge[loop1+60][loop2];
00176           
00177           if(f0>0 && f1>0) {
00178             //Transpos cut
00179             if((transpos[loop1][loop2]>3.5 &&
00180                 transpos[loop1][loop2]<19.5) &&
00181                (transpos[loop1+60][loop2]>3.5 &&
00182                 transpos[loop1+60][loop2]<19.5)) {
00183               //Double Ended Hit
00184               hist[loop1][loop2]->Fill(f0);
00185               hist[loop1+60][loop2]->Fill(f1);
00186             }
00187           }
00188           else if(f0>0) {
00189             //Transpos cut
00190             if((transpos[loop1][loop2]>3.5 &&
00191                 transpos[loop1][loop2]<19.5)) {
00192               hist[loop1][loop2]->Fill(f0);             
00193               //simple zero correct
00194               hist[loop1+60][loop2]->Fill(0);
00195             }
00196           }
00197           else if(f1>0) {
00198             //Transpos cut
00199             if((transpos[loop1+60][loop2]>3.5 &&
00200                 transpos[loop1+60][loop2]<19.5)) {            
00201               hist[loop1+60][loop2]->Fill(f1);
00202               //simple zero correct
00203               hist[loop1][loop2]->Fill(0);
00204             }
00205           }
00206         }
00207       }
00208     }
00209     ind++;
00210   }
00211   
00212   
00213   delete clonearray;
00214   delete cdti;
00215   
00216   save->cd();
00217   //save strip-end histos:
00218   for(int i=0;i<120;i++){
00219     for(int j=0;j<24;j++){
00220       hist[i][j]->Write();
00221     }
00222   }
00223 
00224   cout << "###########" << endl;
00225   cout << "Making Calibration Constants" << endl;
00226   
00227   char outname[256];
00228   sprintf(outname,"PSCalibConstants_%s.dat",fileID);
00229   ofstream out(outname);
00230   
00231   for(int i=0;i<120;i++){
00232     for(int j=0;j<24;j++){
00233       float gain=theGains[i][j];
00234       out << i << "\t" << j << "\t";
00235       out << hist[i][j]->GetEntries() << "\t";
00236       if(gain>0) {
00237         out << hist[i][j]->GetMean()/gain << "\t";
00238         if(hist[i][j]->GetEntries()>0) out << hist[i][j]->GetRMS()/(gain*sqrt(hist[i][j]->GetEntries())) << "\t";
00239         else out << 0 << "\t";
00240       }
00241       else {
00242         gain=80;
00243         out << hist[i][j]->GetMean()/gain << "\t";
00244         if(hist[i][j]->GetEntries()>0) out << hist[i][j]->GetRMS()/(gain*sqrt(hist[i][j]->GetEntries())) << "\t";
00245         else out << 0 << "\t";
00246       }
00247       out << hist[i][j]->GetMean() << "\t";
00248       if(hist[i][j]->GetEntries()>0) out << hist[i][j]->GetRMS()/sqrt(hist[i][j]->GetEntries()) << "\t";
00249       else out << 0 << "\t";
00250       out << hist[i][j]->GetMean() << "\t";
00251       out << hist[i][j]->GetRMS() << endl;
00252     }
00253   }
00254   
00255   cout << "###########" << endl;
00256   cout << "Writing Files" << endl;
00257   save->Close();
00258   out.close();
00259 
00260   cout << "###########" << endl;
00261   cout << "Calibration Done" << endl;
00262   cout << "###########" << endl;
00263 
00264   return 1;
00265 
00266 }

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

Generated on 8 Jul 2019 for loon by  doxygen 1.6.1