//.........这里部分代码省略.........
}
inFile.close();
inFile.open(inputPrefix + "fcal_base_time.txt");
if (inFile.is_open()){
while (getline (inFile, line)){
istringstream iss(line);
iss>>fcal_t_base;
}
}
inFile.close();
inFile.open(inputPrefix + "cdc_base_time.txt");
if (inFile.is_open()){
while (getline (inFile, line)){
istringstream iss(line);
iss>>cdc_t_base;
}
}
inFile.close();
// Do our final step in the timing alignment with tracking
//When the RF is present we can try to simply pick out the correct beam bucket for each of the runs
//First just a simple check to see if we have the appropriate data
bool useRF = false;
double RF_Period = 4.0080161;
TH1I *testHist = Get1DHistogram("HLDetectorTiming", "TAGH_TDC_RF_Compare","Counter ID 001");
if (testHist != NULL){ // Not great since we rely on channel 1 working, but can be craftier later.
useRF = true;
}
ofstream outFile;
TH2I *thisHist;
thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - SC Target Time");
if (useRF) thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - RFBunch Time");
if (thisHist != NULL){
//Statistics on these histograms are really quite low we will have to rebin and do some interpolation
outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::trunc);
outFile.close(); // clear file
outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out | ios::trunc);
outFile.close(); // clear file
int nBinsX = thisHist->GetNbinsX();
int nBinsY = thisHist->GetNbinsY();
TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
for (int i = 1 ; i <= nBinsX; i++){
TH1D *projY = thisHist->ProjectionY("temp", i, i);
// Scan over the histogram
//chose the correct number of bins based on the histogram
float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
float timeWindow = 3; //ns (Full Width)
int binWindow = int(timeWindow / nsPerBin);
double maxEntries = 0;
double maxMean = 0;
for (int j = 1 ; j <= projY->GetNbinsX();j++){
int minBin = j;
int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
double sum = 0, nEntries = 0;
for (int bin = minBin; bin <= maxBin; bin++){
sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
nEntries += projY->GetBinContent(bin);
if (bin == maxBin){
if (nEntries > maxEntries) {
maxMean = sum / nEntries;
maxEntries = nEntries;
int frameStack2_Mall(char* arg){
//Take the arguments and save them into respective strings
std::string infileName, outfileName0, outAllfileName0, outfileName1, outAllfileName1;
std::string inF, outF0, outF1, outAll0, outAll1;
std::string inPrefix, outPrefix;
std::string runs, layers;
std::string runCount;
std::istringstream stm(arg);
inPrefix = "/home/p180f/Do_Jo_Ol_Ma/Analysis/MainProcedure/testMain/rawRoot/";
outPrefix = "/home/p180f/Do_Jo_Ol_Ma/Analysis/MainProcedure/testMain/images/";
outAll0 = "sliceXCuts_allLayers.png";
outAllfileName0 = outPrefix + outAll0;
std::cout << outAll0 << " created\n";
outAll1 = "projYCuts_allLayers.png";
outAllfileName1 = outPrefix + outAll1;
std::cout << outAll1 << " created\n";
const int width=480; //width of the raw image
const int height=640; //height of the raw image
TH2I *frameHistoAll = new TH2I("frameHistoAll","Stacked Frames After Edge Cuts",width/4,0,width,height/4,0,height); //histogram for the stacked images
TH1I *chamber1All = new TH1I("chamber1All","Chamber 1 After Edge Cuts",width/4,0,width);//histogram for chamber 1 (the top one)
TH1I *chamber2All = new TH1I("chamber2All","Chamber 2 After Edge Cuts",width/4,0,width);//histogram for chamber 2
TH1I *chamber3All = new TH1I("chamber3All","Chamber 3 After Edge Cuts",width/4,0,width);//histogram for chamber 3
TH1I *chamber4All = new TH1I("chamber4All","Chamber 4 After Edge Cuts",width/4,0,width);//histogram for chamber 4 (the bottom one)
TCanvas *projCAll = new TCanvas("projCAll","",0,0,800,600);
TCanvas *pc2All = new TCanvas("pc2All", "Stack of 4 Layer Runs", 0, 0, 800, 600);
while (true) {
if (std::getline(stm, layers, ' ')) {
//create the output root file
outF0 = "sliceXCuts_" + layers + "layers.png";
outfileName0 = outPrefix + outF0;
std::cout << outF0 << " created\n";
outF1 = "projYCuts_" + layers + "layers.png";
outfileName1 = outPrefix + outF1;
std::cout << outF1 << " created\n";
//load the input root files
TChain *chain = new TChain("fourChamTree");
for (int i=0; ; i++) {
runCount = std::to_string(i);
inF = "run" + runCount + "_" + layers + "layers.root";
infileName = inPrefix + inF;
ifstream fin;
fin.open(infileName.c_str());
if (!fin.fail()) {
fin.close();
chain->Add(infileName.c_str());
std::cout << "Got " << inF << std::endl;
} else break;
}
int x=-10; //x from file
int y=-10; //y from file
int intensity=-10; //pixle intensity from file
int pNum=0;//the order in which the frame was processed
//the 2d array which will store each frame of image data.
int frame[480][640]={0};
//variables
int UNIXtime=0;
float tdc[2]={-10,-10};
//TTree *T = new TTree("T","TTree of muplus data");
//add the 'branches' to the tree we will now read in
chain->SetBranchAddress("pNum",&pNum); //branch for the frame number
chain->SetBranchAddress("frame",&frame); //branch for frame data
TH2I *frameHisto = new TH2I("frameHisto","Stacked Frames After Edge Cuts",width/4,0,width,height/4,0,height); //histogram for the stacked images
TH1I *chamber1 = new TH1I("chamber1","Chamber 1 After Edge Cuts",width/4,0,width);//histogram for chamber 1 (the top one)
TH1I *chamber2 = new TH1I("chamber2","Chamber 2 After Edge Cuts",width/4,0,width);//histogram for chamber 2
TH1I *chamber3 = new TH1I("chamber3","Chamber 3 After Edge Cuts",width/4,0,width);//histogram for chamber 3
TH1I *chamber4 = new TH1I("chamber4","Chamber 4 After Edge Cuts",width/4,0,width);//histogram for chamber 4 (the bottom one)
//loop over all data in chain
Int_t nevent = chain->GetEntries(); //get the number of entries in the TChain
for (Int_t i=0;i<nevent;i++) {
chain->GetEntry(i);
for(int x=0;x<width;x++){
for(int y=0;y<height;y++){
if(frame[x][y]>0){
frameHisto->Fill(x,y,frame[x][y]);
frameHistoAll->Fill(x,y,frame[x][y]);
if(y>580 && y<610){
chamber1->Fill(x,frame[x][y]);
chamber1All->Fill(x,frame[x][y]);
}
else if(y>400 && y<440){
chamber2->Fill(x,frame[x][y]);
chamber2All->Fill(x,frame[x][y]);
}
//.........这里部分代码省略.........
void plot_pad_size_in_layer(TString digiPar="trd.v13/trd_v13g.digi.par", Int_t nlines=1, Int_t nrows_in_sec=0, Int_t alllayers=1)
{
gStyle->SetPalette(1,0);
gROOT->SetStyle("Plain");
gStyle->SetPadTickX(1);
gStyle->SetPadTickY(1);
gStyle->SetOptStat(kFALSE);
gStyle->SetOptTitle(kFALSE);
Bool_t read = false;
TH2I *fLayerDummy = new TH2I("LayerDummy","",1200,-600,600,1000,-500,500);
fLayerDummy->SetXTitle("x-coordinate [cm]");
fLayerDummy->SetYTitle("y-coordinate [cm]");
fLayerDummy->GetXaxis()->SetLabelSize(0.02);
fLayerDummy->GetYaxis()->SetLabelSize(0.02);
fLayerDummy->GetZaxis()->SetLabelSize(0.02);
fLayerDummy->GetXaxis()->SetTitleSize(0.02);
fLayerDummy->GetXaxis()->SetTitleOffset(1.5);
fLayerDummy->GetYaxis()->SetTitleSize(0.02);
fLayerDummy->GetYaxis()->SetTitleOffset(2);
fLayerDummy->GetZaxis()->SetTitleSize(0.02);
fLayerDummy->GetZaxis()->SetTitleOffset(-2);
TString title;
TString title1, title2, title3;
TString buffer;
TString firstModule = "";
Int_t blockCounter(0), startCounter(0); // , stopCounter(0);
Double_t msX(0), msY(0), mpX(0), mpY(0), mpZ(0), psX(0), psY(0);
Double_t ps1X(0), ps1Y(0), ps2X(0), ps2Y(0), ps3X(0), ps3Y(0);
Int_t modId(0), layerId(0);
Double_t sec1(0), sec2(0), sec3(0);
Double_t row1(0), row2(0), row3(0);
std::map<float, TCanvas*> layerView;// map key is z-position of modules
std::map<float, TCanvas*>::iterator it;
ifstream digipar;
digipar.open(digiPar.Data(), ifstream::in);
while (digipar.good()) {
digipar >> buffer;
//cout << "(" << blockCounter << ") " << buffer << endl;
if (blockCounter == 19)
firstModule = buffer;
if (buffer == (firstModule + ":")){
//cout << buffer << " <===========================================" << endl;
read = true;
}
if (read) {
startCounter++;
if (startCounter == 1) // position of module position in x
{
modId = buffer.Atoi();
layerId = (modId & (15 << 4)) >> 4; // from CbmTrdAddress.h
}
if (startCounter == 5) // position of module position in x
mpX = buffer.Atof();
if (startCounter == 6) // position of module position in y
mpY = buffer.Atof();
if (startCounter == 7) // position of module position in z
mpZ = buffer.Atof();
if (startCounter == 8) // position of module size in x
msX = buffer.Atof();
if (startCounter == 9) // position of module size in y
msY = buffer.Atof();
if (startCounter == 12) // sector 1 size in y
sec1 = buffer.Atof();
if (startCounter == 13) // position of pad size in x - do not take the backslash (@14)
ps1X = buffer.Atof();
if (startCounter == 15) // position of pad size in y
ps1Y = buffer.Atof();
if (startCounter == 17) // sector 2 size in y
sec2 = buffer.Atof();
if (startCounter == 18) // position of pad size in x
{
ps2X = buffer.Atof();
psX = ps2X; // for backwards compatibility - sector 2 is default sector
}
if (startCounter == 19) // position of pad size in y
{
ps2Y = buffer.Atof();
psY = ps2Y; // for backwards compatibility - sector 2 is default sector
}
if (startCounter == 21) // sector 3 size in y
sec3 = buffer.Atof();
if (startCounter == 22) // position of pad size in x
ps3X = buffer.Atof();
if (startCounter == 23) // position of pad size in y
ps3Y = buffer.Atof();
// if (startCounter == 23) // last element
// {
// printf("moduleId : %d, %d\n", modId, layerId);
// printf("pad size sector 1: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps1X, ps1Y, ps1X*ps1Y);
// printf("pad size sector 2: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps2X, ps2Y, ps2X*ps2Y);
// printf("pad size sector 3: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps3X, ps3Y, ps3X*ps3Y);
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论