#include #include #include #include #include #include #include void plot_dedx(int firstrun, int lastrun){ // Plot dE/dx for sets of 10 subruns vs. time for an inital run - final run // For cosmic triggers, negative tracks with nodes>35 // To run: // root // .L stability.C // plot_dedx(4979,5115) // for runs 4979-5115 // // Path to files may need to change, depending on run. // Check you can access the files // Ouput is 3 plots vs. time // dE/dx fitted mean value. Errors set to the standard error on the mean // Written to dedx_runXXXX_runYYYY_cosmic.gif // dE/dx RMS/mean. Errors set from # of entries used in histogram // Written to dedx_sigma_runXXXX_runYYYY_cosmic.gif // First run's histogram, for debugging, not saved to a gif gStyle->SetOptStat(111111111); SetupStyle(); gStyle->SetOptTitle(1); char title[256]; char filename[256]; char plotcommand[256]; char plotcommandtime[256]; char cuts[256]; TH1F *dedxrun[6000][12]; TH1F *dedxhist[6000][100][12]; TH1F *dedxtemp[12]; TGraphErrors *gdEdx[12]; TGraphErrors *gdEdx_sigma[12]; TH1F *dedxrun[6000][12]; Double_t truelastsub[6000]; Double_t time[6000]; double meanval[12]={600.,375.,375.,375.,375.,400.,400.,410,410.,425.,425.,500.}; // dE/dx approx mean values read off of an old plot==> would be better from a calc fvsrun = new TCanvas("fvsrun","dE/dx vs run, cosmics",100,100,1200,800); fvsrun->Divide(3,4); fvsrun->cd(1); fvssigma = new TCanvas("fvssigma","RMS(dE/dx) vs run, cosmics",100,100,1200,800); fvssigma->Divide(3,4); fvssigma->cd(1); fvsab = new TCanvas("fvsab","dE/dx first run",100,100,1200,800); fvsab->Divide(3,4); fvsab->cd(1); int firstsub = 0; int lastsub = 100; int batchsize = 20; int batch = (lastsub-firstsub)/batchsize; for(int ctr=0;ctr<=(lastrun-firstrun)*batch+lastsub/batchsize;ctr++){ truelastsub[ctr]=0.; for(int mom=0;mom<12;mom++){ sprintf(title,"dedxrun_%d_%d",ctr,mom); dedxrun[ctr][mom] = new TH1F(title,title,100,0.1,1000); } } int batchctr = 0; for(int run=firstrun;run<=lastrun;run++){ batchctr=(run-firstrun)*batch; for(int sub=firstsub;sub<=lastsub;sub++){ int ctr=batchctr+sub/batchsize; // Loop over the relevant files sprintf(filename,"/neut/data11/kendallm/irods/first_pass/34/%d/dq-tpc-bc_%08d_%04d.root",run,run,sub); // sprintf(filename,"/t2k/dataquality/dq-tpc-bc_%08d_%04d.root",run,sub); if(FileExists(filename)){ TFile *pg = new TFile(filename); // if(!pg){ //did not help the Error file not found message if(!pg->IsZombie()){ // catches cases where file did not close properly if(pg->GetNkeys()!=0){ // check that there are keys in the file if(pg->GetNkeys()==0){ // std::cout << "No keys: " << filename <IsOpen()){ TChain *c = (TChain*)pg.Get("bctpc"); if(c->GetEntries()!=0){ for(int mom=0;mom<12;mom++){ sprintf(title,"dedx_%d_%d_%d",run,sub,mom); dedxhist[run][sub][mom] = new TH1F(title,title,100,0.1,1000); sprintf(title,"dedxtemp_%d_%d_%d",mom); dedxtemp[mom] = new TH1F(title,title,100,0.1,1000); Double_t lowmom =0.; Double_t highmom = 1500.; // set up momentum bin delimeters if(mom<9&&mom>0){ lowmom = 200.+100*mom; highmom = lowmom+100.; }elseif(mom==0){ lowmom = 0; highmom = 200.; }elseif(mom==9){ lowmom = 1000.; highmom = lowmom+200.; }elseif(mom==10){ lowmom = 1200.; highmom = lowmom+300.; }elseif(mom==11){ lowmom = 1500.; highmom = 2000.; } for(int trk=0;trk<2;trk++){ sprintf(cuts,"momentum[%d]>%d&&momentum[%d]<=%d&&ctrig&&ntrack>=%d&&nnodes>35&&rho[%d]<0.",trk,lowmom,trk,highmom,trk,trk); sprintf(plotcommand,"dedx[%d]>>dedx_%d_%d_%d",trk,run,sub,mom); c->Draw(plotcommand,cuts); dedxrun[ctr][mom]->Add(dedxhist[run][sub][mom]); // add the histograms for 10 subrun files together } } //momentum loop sprintf(title,"time_%d_%d",run,sub); sprintf(plotcommandtime,"timestamp>>time_%d_%d",run,sub); fvsab->cd(1); c->Draw("timestamp>>temp"); TH1F *temp = (TH1F*)gDirectory->Get("temp"); // std::cout << "temp run:" << run << " sub " << sub << " " << temp->GetMean() << std::endl; //std::cout << "batchctr" << batchctr << std::endl; time[ctr]=temp->GetMean()+time[ctr]; truelastsub[ctr]++; //temp->Draw(); } // non zero entries in tree } // open file } // file has keys } // filezombie } // file exists } //subrun } //run for(int ctr=0;ctr<=(lastrun-firstrun)*batch+lastsub/batchsize;ctr++){ Double_t numbersub=truelastsub[ctr]; if(numbersub>0.1){ time[ctr]=time[ctr]/numbersub; // average time for # subevents used // std::cout << "#sub " << truelastsub[ctr] << std::endl; } } for(int mom=0;mom<12;mom++){ gdEdx[mom] = new TGraphErrors(); sprintf(title,"gdEdx_%d",mom); gdEdx[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx[mom]->SetMarkerStyle(24); gdEdx_sigma[mom] = new TGraphErrors(); sprintf(title,"gdEdx_sigma_%d",mom); gdEdx_sigma[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx_sigma[mom]->SetMarkerStyle(24); int ctr=0; for(int brun=0;brun<=(lastrun-firstrun)*batch+lastsub/batchsize;brun++){ if(dedxrun[brun][mom]){ if(dedxrun[brun][mom]->GetEntries()>50.){ // only fill if >50 entries in histogram // currently about 100 events / subrun if(truelastsub[brun]>0){ dedxrun[brun][mom]->Fit("gaus","Q","",dedxrun[brun][mom]->GetMean()-dedxrun[brun][mom]->GetRMS()*2.,dedxrun[brun][mom]->GetMean()+dedxrun[brun][mom]->GetRMS()*1.); // fit dE/dx TF1 *func = dedxrun[brun][mom]->GetFunction("gaus"); gdEdx[mom]->SetPoint(ctr,time[brun],func->GetParameter(1)); // Extract dE/dx fitted mean gdEdx[mom]->SetPointError(ctr,10.,func->GetParameter(2)/(TMath::Sqrt(dedxrun[brun][mom]->GetEntries()))); // fitted RMS/sqrt(entries), standard error on the mean gdEdx_sigma[mom]->SetPoint(ctr,time[brun],func->GetParameter(2)/func->GetParameter(1)); // fitted RMS/ mean gdEdx_sigma[mom]->SetPointError(ctr,0.,func->GetParameter(2)/(func->GetParameter(1)*(TMath::Sqrt(dedxrun[brun][mom]->GetEntries())))); // std::cout << " error? " << func->GetParameter(2)/(TMath::Sqrt(dedxrun[brun][mom]->GetEntries())) <entries } // TH1 exists } // batch ctr } // momentum for(int mom=0;mom<12;mom++){ p=fvsrun->cd(mom+1); if(mom==0) gdEdx[mom]->SetTitle("Mean dE/dx value for momentum<200"); if(mom==1) gdEdx[mom]->SetTitle("Mean dE/dx value for 200>momentum>300"); if(mom==2) gdEdx[mom]->SetTitle("Mean dE/dx value for 300>momentum>400"); if(mom==3) gdEdx[mom]->SetTitle("Mean dE/dx value for 400>momentum>500"); if(mom==4) gdEdx[mom]->SetTitle("Mean dE/dx value for 500>momentum>600"); if(mom==5) gdEdx[mom]->SetTitle("Mean dE/dx value for 600>momentum>700"); if(mom==6) gdEdx[mom]->SetTitle("Mean dE/dx value for 700>momentum>800"); if(mom==7) gdEdx[mom]->SetTitle("Mean dE/dx value for 800>momentum>900"); if(mom==8) gdEdx[mom]->SetTitle("Mean dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx[mom]->SetTitle("Mean dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx[mom]->SetTitle("Mean dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx[mom]->SetTitle("Mean dE/dx value for 1500>momentum>2000"); // Set plotting time axis gdEdx[mom]->GetXaxis()->SetTimeDisplay(1); gdEdx[mom]->GetXaxis()->SetTimeFormat("#splitline{%b-%d}{%H:%M}"); gdEdx[mom]->GetXaxis()->SetLabelOffset(0.05); gdEdx[mom]->GetXaxis()->SetLabelSize(0.06); gdEdx[mom]->GetYaxis()->SetLabelSize(0.06); // Set maximum and minimum plot range gdEdx[mom]->SetMaximum(meanval[mom]*1.3); gdEdx[mom]->SetMinimum(meanval[mom]*0.7); gdEdx[mom]->Draw("ap"); // Draw lines at approx expected mean value of dE/dx and +/- 10% TLine *line1 = new TLine(time[0],meanval[mom],time[(lastrun-firstrun)*batch],meanval[mom]); TLine *linep = new TLine(time[0],meanval[mom]*1.1,time[(lastrun-firstrun)*batch],meanval[mom]*1.1); TLine *linem = new TLine(time[0],meanval[mom]*0.9,time[(lastrun-firstrun)*batch],meanval[mom]*0.9); line1->SetLineColor(2); linep->SetLineColor(2); linem->SetLineColor(2); linep->SetLineStyle(2); linem->SetLineStyle(2); line1->Draw("same"); linep->Draw("same"); linem->Draw("same"); } // momentum fvsrun->Print(Form("dedx_run%d_run%d_cosmic.gif",firstrun,lastrun)); for(int mom=0;mom<12;mom++){ pp=fvssigma->cd(mom+1); // gdEdx_sigma[mom]->Print(); TLine *line1 = new TLine(time[0],0.1,time[(lastrun-firstrun)*batch],0.1); line1->SetLineColor(2); line1->SetLineStyle(2); if(mom==0) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for momentum<200"); if(mom==1) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 200>momentum>300"); if(mom==2) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 300>momentum>400"); if(mom==3) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 400>momentum>500"); if(mom==4) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 500>momentum>600"); if(mom==5) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 600>momentum>700"); if(mom==6) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 700>momentum>800"); if(mom==7) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 800>momentum>900"); if(mom==8) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1500>momentum>2000"); // Set plotting time axis gdEdx_sigma[mom]->GetXaxis()->SetTimeDisplay(1); gdEdx_sigma[mom]->GetXaxis()->SetTimeFormat("#splitline{%b-%d}{%H:%M}"); gdEdx_sigma[mom]->GetXaxis()->SetLabelOffset(0.05); gdEdx_sigma[mom]->GetXaxis()->SetLabelSize(0.06); gdEdx_sigma[mom]->GetYaxis()->SetLabelSize(0.06); // Set maximum and minimum plot range gdEdx_sigma[mom]->SetMaximum(0.2); gdEdx_sigma[mom]->SetMinimum(0.); gdEdx_sigma[mom]->Draw("ap"); line1->Draw("same"); } fvssigma->Print(Form("dedx_sigma_run%d_run%d_cosmic.gif",firstrun,lastrun)); for(int mom=0;mom<12;mom++){ ppp=fvsab->cd(mom+1); dedxrun[0][mom]->Draw(); } } void plot_calib_histo_cosmic(int firstrun, int lastrun){ // Plot calib constants for sets of 10 subruns vs. time for an inital // run - final run // For cosmic triggers // Plot the calibration constants from the orginal histograms // To run: // root // .L stability.C // plot_dedx(4979,5115) // for runs 4979-5115 // // Path to files may need to change, depending on run. // Check you can access the files // Ouput is plots of the cosmic, beam constants vs. time // gStyle->SetOptStat(111111111); SetupStyle(); gStyle->SetOptTitle(1); char title[256]; char histogram[256]; char filename[256]; char plotcommand[256]; char plotcommandtime[256]; char cuts[256]; TH1F *calib[6000][72]; TGraphErrors *calib_cosmic[72]; Double_t time[6000]; Double_t truelastsub[6000]; Double_t maxcalib = 500.; Double_t mincalib = 100.; fvsrun10 = new TCanvas("fvsrun10","TPC1 calib (cosmics)",100,100,1200,800); fvsrun10->Divide(3,4); fvsrun10->cd(1); fvsrun20 = new TCanvas("fvsrun20","TPC2 calib (cosmics)",100,100,1200,800); fvsrun20->Divide(3,4); fvsrun20->cd(1); fvsrun30 = new TCanvas("fvsrun30","TPC3 calib (cosmics)",100,100,1200,800); fvsrun30->Divide(3,4); fvsrun30->cd(1); fvsab = new TCanvas("fvsab","first cosmic hist and fit",100,100,1200,800); fvsab->Divide(3,4); fvsab->cd(1); int firstsub = 0; int lastsub = 100; int batchsize = 100; int batch = (lastsub-firstsub)/batchsize; // Initialize subrun counter and calib TH1F array for(int ctr=0;ctr<=(lastrun-firstrun)*batch+lastsub/batchsize;ctr++){ truelastsub[ctr]=0.; for(int mm=0;mm<72;mm++){ sprintf(title,"calib_%d_%d",ctr,mm); calib[ctr][mm] = new TH1F(title,title,100,0,2000); } } int batchctr=0; for(int run=firstrun;run<=lastrun;run++){ batchctr=(run-firstrun)*batch; for(int sub=firstsub;sub<=lastsub;sub++){ int ctr=batchctr+sub/batchsize; // Loop over the relevant files sprintf(filename,"/neut/data11/kendallm/irods/first_pass/34/%d/dq-tpc-bc_%08d_%04d.root",run,run,sub); // on the soff machine // sprintf(filename,"/t2k/dataquality/dq-tpc-bc_%08d_%04d.root",run,sub); if(FileExists(filename)){ TFile *pg = new TFile(filename); // if(!pg){ //did not help the Error file not found message if(!pg->IsZombie()){ // catches cases where file did not close properly if(pg->GetNkeys()!=0){ // check that there are keys in the file if(pg->GetNkeys()==0){ // std::cout << "No keys: " << filename <IsOpen()){ TChain *c = (TChain*)pg.Get("bctpc"); if(c->GetEntries()!=0){ sprintf(title,"time_%d_%d",run,sub); sprintf(plotcommandtime,"timestamp>>time_%d_%d",run,sub); fvsab->cd(1); c->Draw("timestamp>>tempp"); TH1F *tempp = (TH1F*)gDirectory->Get("tempp"); time[ctr]=tempp->GetMean()+time[ctr]; for(int i=0;i<72;i++){ sprintf(histogram,"CosmicCalib_Mod_%d",i); TH1F *temp = (TH1F*)gDirectory->Get(histogram); if(i==10) std::cout << "run, subrun" << run << " " << sub << " main counter" << ctr << " calib "<< calib[ctr][i]->GetEntries() << std::endl; calib[ctr][i]->Add(temp); temp->Reset(); } // loop MM truelastsub[ctr]++; } // chain !=0 } // open file } // file has keys } // filezombie } // file exists } //subrun } //run for(int ctr=0;ctr<=(lastrun-firstrun)*batch+lastsub/batchsize;ctr++){ Double_t numbersub=truelastsub[ctr]; if(numbersub>0.1){ time[ctr]=time[ctr]/numbersub; // average time for # subevents used } } for(int mm=0;mm<72;mm++){ calib_cosmic[mm] = new TGraphErrors(); sprintf(title,"calib_cosmic_%d",mm); calib_cosmic[mm]->SetName(title); // if Name is not set, will not appear in output file calib_cosmic[mm]->SetMarkerStyle(24); int ctr=0; for(int brun=0;brun<=(lastrun-firstrun)*batch+lastsub/batchsize;brun++){ // cosmic constants // std::cout << "counter " << brun << std::endl; if(calib[brun][mm]){ if(calib[brun][mm]->GetEntries()>50.){ // only fill if >50 entries in histogram // 10 subruns appears to be sufficienty for cosmics if(truelastsub[brun]>0){ calib[brun][mm]->Fit("landau","Q","",0,2000); // Fit Landau distribution TF1 *func = calib[brun][mm]->GetFunction("landau"); calib_cosmic[mm]->SetPoint(ctr,time[brun],func->GetParameter(1)); // calib_cosmic[mm]->SetPointError(ctr,10.,func->GetParError(1)); // ctr++; if(func->GetParameter(1)>maxcalib){ std::cout << "Point " << ctr << " above max range on graph for MM " << mm << std::endl; } if(func->GetParameter(1)GetParameter(1)<200.){ std::cout << "Point " << brun << " has a low Landau MM " << mm << " with" << func->GetParameter(1) << " mean " << calib[brun][mm]->GetMean() << std::endl; } // if(mm==49) std::cout << " mm " << mm << " ctr " << brun << " #sub " << truelastsub[brun] << " " << calib[brun][mm]->GetMean() << std::endl; // seems to have non zero entries }else{ // seems to have non zero entries, but these do not correspond to proper subruns. // if(mm==65) std::cout << " mm " << mm << " ctr " << brun << " #sub " << truelastsub[brun] << " " << calib[brun][mm]->GetMean() << std::endl; // seems to have non zero entries } } // >entries } // TH1 exists } // batch ctr } // mm for(int mm=0;mm<12;mm++){ p=fvsrun10->cd(mm+1); sprintf(title,"Calib for cosmic trig, TPC 1 MM %d",mm); calib_cosmic[mm]->SetTitle(title); calib_cosmic[mm]->SetMarkerColor(1); //RP0 black, circle calib_cosmic[mm+12]->SetMarkerColor(2); //RP1 red calib_cosmic[mm+12]->SetMarkerStyle(27); //RP1 open diamond // Set plotting time axis calib_cosmic[mm]->GetXaxis()->SetTimeDisplay(1); calib_cosmic[mm]->GetXaxis()->SetTimeFormat("#splitline{%b-%d}{%H:%M}"); calib_cosmic[mm]->GetXaxis()->SetLabelOffset(0.05); calib_cosmic[mm]->GetXaxis()->SetLabelSize(0.06); calib_cosmic[mm]->GetYaxis()->SetLabelSize(0.06); // Set maximum and minimum plot range calib_cosmic[mm]->SetMaximum(maxcalib); calib_cosmic[mm]->SetMinimum(mincalib); calib_cosmic[mm+12]->SetMaximum(maxcalib); calib_cosmic[mm+12]->SetMinimum(mincalib); calib_cosmic[mm]->Draw("ap"); calib_cosmic[mm+12]->Draw("same,p"); // Need to check if there are points outside range } // mmentum fvsrun10->Print(Form("calib_tpc1_run%d_run%d_cosmic.gif",firstrun,lastrun)); for(int mm=0;mm<12;mm++){ p=fvsrun20->cd(mm+1); Int_t offset = 24; sprintf(title,"Calib for cosmic trig, TPC 2 MM %d",mm); calib_cosmic[mm+offset]->SetTitle(title); calib_cosmic[mm+offset]->SetMarkerColor(1); //RP0 black, circle calib_cosmic[mm+offset+12]->SetMarkerColor(2); //RP1 red calib_cosmic[mm+offset+12]->SetMarkerStyle(27); //RP1 ?? // Set plotting time axis calib_cosmic[mm+offset]->GetXaxis()->SetTimeDisplay(1); calib_cosmic[mm+offset]->GetXaxis()->SetTimeFormat("#splitline{%b-%d}{%H:%M}"); calib_cosmic[mm+offset]->GetXaxis()->SetLabelOffset(0.05); calib_cosmic[mm+offset]->GetXaxis()->SetLabelSize(0.06); calib_cosmic[mm+offset]->GetYaxis()->SetLabelSize(0.06); // Set maximum and minimum plot range calib_cosmic[mm+offset]->SetMaximum(maxcalib); calib_cosmic[mm+offset]->SetMinimum(mincalib); calib_cosmic[mm+offset+12]->SetMaximum(maxcalib); calib_cosmic[mm+offset+12]->SetMinimum(mincalib); calib_cosmic[mm+offset]->Draw("ap"); calib_cosmic[mm+offset+12]->Draw("same,p"); // Need to check if there are points outside range } // 1 RP of MM (12) fvsrun20->Print(Form("calib_tpc2_run%d_run%d_cosmic.gif",firstrun,lastrun)); for(int mm=0;mm<12;mm++){ p=fvsrun30->cd(mm+1); Int_t offset = 48; sprintf(title,"Calib for cosmic trig, TPC 3 MM %d",mm); calib_cosmic[mm+offset]->SetTitle(title); calib_cosmic[mm+offset]->SetMarkerColor(1); //RP0 black, circle calib_cosmic[mm+offset+12]->SetMarkerColor(2); //RP1 red calib_cosmic[mm+offset+12]->SetMarkerStyle(27); //RP1 ?? // Set plotting time axis calib_cosmic[mm+offset]->GetXaxis()->SetTimeDisplay(1); calib_cosmic[mm+offset]->GetXaxis()->SetTimeFormat("#splitline{%b-%d}{%H:%M}"); calib_cosmic[mm+offset]->GetXaxis()->SetLabelOffset(0.05); calib_cosmic[mm+offset]->GetXaxis()->SetLabelSize(0.06); calib_cosmic[mm+offset]->GetYaxis()->SetLabelSize(0.06); // Set maximum and minimum plot range calib_cosmic[mm+offset]->SetMaximum(maxcalib); calib_cosmic[mm+offset]->SetMinimum(mincalib); calib_cosmic[mm+offset+12]->SetMaximum(maxcalib); calib_cosmic[mm+offset+12]->SetMinimum(mincalib); calib_cosmic[mm+offset]->Draw("ap"); calib_cosmic[mm+offset+12]->Draw("same,p"); // Need to check if there are points outside range } // 1 RP of MM (12) fvsrun30->Print(Form("calib_tpc3_run%d_run%d_cosmic.gif",firstrun,lastrun)); /* for(int mm=0;mm<12;mm++){ ppp=fvsab->cd(mm+1); calib[0][mm]->Draw(); } */ fvsab->cd(1); calib[10][3]->Draw(); fvsab->cd(2); calib[1][49]->Draw(); fvsab->cd(3); calib[33][49]->Draw(); fvsab->cd(4); calib[127][49]->Draw(); return; } void plot_dedx_run(int firstrun, int lastrun){ // TChain c("bc"); // declare TChain, must the the right tree to know to add it // warning for the file size gStyle->SetOptStat(111111111); gStyle->SetOptTitle(1); char title[256]; char filename[256]; char plotcommand[256]; char cuts[256]; // int run = 4980; TH1F *dedxrun[6000][12]; TH1F *dedxhist[6000][30][12]; TGraphErrors *gdEdx[12]; TGraphErrors *gdEdx_sigma[12]; TH1F *dedxrun[6000][12]; TH1F *dedxhist[6000][30][12]; TGraphErrors *gdEdx[12]; TGraphErrors *gdEdx_sigma[12]; TH1F *dedxrunb[6000][12]; TH1F *dedxhistb[6000][30][12]; TGraphErrors *gdEdxb[12]; TGraphErrors *gdEdxb_sigma[12]; double meanval[12]={600.,375.,375.,375.,375.,400.,400.,410,410.,425.,425.,500.}; // dE/dx approx mean values read off of an old plot==> would be better from a calc fvsrun = new TCanvas("fvsrun","dE/dx vs run, cosmics",100,100,1200,800); fvsrun->Divide(3,4); fvsrun->cd(1); fvssigma = new TCanvas("fvssigma","RMS(dE/dx) vs run, cosmics",100,100,1200,800); fvssigma->Divide(3,4); fvssigma->cd(1); fvsrunb = new TCanvas("fvsrunb","dE/dx vs run, beam",100,100,1200,800); fvsrunb->Divide(3,4); fvsrunb->cd(1); fvssigmab = new TCanvas("fvssigmab","RMS(dE/dx) vs run, beam",100,100,1200,800); fvssigmab->Divide(3,4); fvssigmab->cd(1); fvsab = new TCanvas("fvsab","dE/dx first run",100,100,1200,800); fvsab->Divide(3,4); fvsab->cd(1); int firstsub = 0; int lastsub = 9; // need to check for the done file first. for(int run=firstrun;run<=lastrun;run++){ for(int mom=0;mom<12;mom++){ sprintf(title,"dedxrun_%d_%d",run,mom); dedxrun[run][mom] = new TH1F(title,title,100,0.1,1000); sprintf(title,"dedxrunb_%d_%d",run,mom); dedxrunb[run][mom] = new TH1F(title,title,100,0.1,1000); } for(int sub=firstsub;sub<=lastsub;sub++){ // c.Add(Form("/t2k/dataquality/dq-tpc-bc_%08d_%04d.root",run,sub)); // pad with 0s 8 wide for run, pad with subrun 4 wide // open file and save a histogram for all momentum and then momentum bins // TFile *pg = (TFile*)gROOT->GetListOfFiles()->FindObject(Form("dq-tpc-bc_%08d_%04d.root",run,sub)); // did not help the Error file not found message // if(!pg){ did not help the Error file not found message sprintf(filename,"dq-tpc-bc_%08d_%04d.root",run,sub); TFile *pg = new TFile(filename); // TFile *pg = new TFile(Form("dq-tpc-bc_%08d_%04d.root",run,sub)); std::cout << " after TFile" << std::endl; if(pg->IsOpen()){ std::cout << " after Open" << std::endl; if(!pg->IsZombie()){ // made no difference as well, error is from TFile creation TChain *c = (TChain*)pg.Get("bc"); c->Draw("timestamp"); for(int mom=0;mom<12;mom++){ sprintf(title,"dedx_%d_%d_%d",run,sub,mom); dedxhist[run][sub][mom] = new TH1F(title,title,100,0.1,1000); sprintf(title,"dedxb_%d_%d_%d",run,sub,mom); dedxhistb[run][sub][mom] = new TH1F(title,title,100,0.1,1000); // cint interpreter title is equivalent to the variable name Double_t lowmom =0.; Double_t highmom = 1500.; if(mom<9&&mom>0){ lowmom = 200.+100*mom; highmom = lowmom+100.; }elseif(mom==0){ lowmom = 0; highmom = 200.; }elseif(mom==9){ lowmom = 1000.; highmom = lowmom+200.; }elseif(mom==10){ lowmom = 1200.; highmom = lowmom+300.; }elseif(mom==11){ lowmom = 1500.; highmom = 2000.; } // sprintf(cuts,"momentum>%d&&momentum<=%d",lowmom,highmom); sprintf(cuts,"momentum>%d&&momentum<=%d&&ctrig",lowmom,highmom); // c->Draw(plotcommand,"momentum>0&&momentum<=1500"); sprintf(plotcommand,"dedx>>dedx_%d_%d_%d",run,sub,mom); // c->Draw(plotcommand,"","title");// this creates a new histogram w/ that title c->Draw(plotcommand,cuts); // sprintf(cuts,"momentum>%d&&momentum<=%d",lowmom,highmom); sprintf(cuts,"momentum>%d&&momentum<=%d&&btrig",lowmom,highmom); sprintf(plotcommand,"dedx>>dedxb_%d_%d_%d",run,sub,mom); c->Draw(plotcommand,cuts); dedxrun[run][mom]->Add(dedxhist[run][sub][mom]); dedxrunb[run][mom]->Add(dedxhistb[run][sub][mom]); } //momentum loop } // file loop } // zombie } // subrun } // run // } // if!pg // also do as a function of time for(int mom=0;mom<12;mom++){ gdEdx[mom] = new TGraphErrors(); sprintf(title,"gdEdx_%d",mom); gdEdx[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx[mom]->SetMarkerStyle(24); gdEdxb[mom] = new TGraphErrors(); sprintf(title,"gdEdxb_%d",mom); gdEdxb[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdxb[mom]->SetMarkerStyle(24); gdEdx_sigma[mom] = new TGraphErrors(); sprintf(title,"gdEdx_sigma_%d",mom); gdEdx_sigma[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx_sigma[mom]->SetMarkerStyle(24); gdEdxb_sigma[mom] = new TGraphErrors(); sprintf(title,"gdEdx_sigma_%d",mom); gdEdxb_sigma[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdxb_sigma[mom]->SetMarkerStyle(24); int ctr=0; int ctrb=0; for(int run=firstrun;run<=lastrun;run++){ if(dedxrun[run][mom]){ if(dedxrun[run][mom]->GetEntries()>20.){ // only fill if >20 entries in histogram // currently about 100 events / subrun dedxrun[run][mom]->Fit("gaus","Q","",dedxrun[run][mom]->GetMean()-dedxrun[run][mom]->GetRMS()*2.,dedxrun[run][mom]->GetMean()+dedxrun[run][mom]->GetRMS()*1.); TF1 *func = dedxrun[run][mom]->GetFunction("gaus"); /* debug std::cout << " p0 norm" << func->GetParameter(0) << std::endl; std::cout << " p1 m" << func->GetParameter(1) << std::endl; std::cout << " p2 s " << func->GetParameter(2) << std::endl; */ gdEdx[mom]->SetPoint(ctr,run,func->GetParameter(1)); // dE/dx fitted mean gdEdx[mom]->SetPointError(ctr,0.5,func->GetParameter(2)/(TMath::Sqrt(dedxrun[run][mom]->GetEntries()))); // fitted RMS/sqrt(entries) gdEdx_sigma[mom]->SetPoint(ctr,run,func->GetParameter(2)/func->GetParameter(1)); // fitted RMS/ mean gdEdx_sigma[mom]->SetPointError(ctr,0.,0.); ctr++; } } // TH1F exists } // run } // momentum for(int mom=0;mom<12;mom++){ p=fvsrun->cd(mom+1); // Double_t meanval[mom] = 440.; gdEdx[mom]->GetXaxis()->SetTitle("Run number"); if(mom==0) gdEdx[mom]->SetTitle("Mean dE/dx value for momentum<200"); if(mom==1) gdEdx[mom]->SetTitle("Mean dE/dx value for 200>momentum>300"); if(mom==2) gdEdx[mom]->SetTitle("Mean dE/dx value for 300>momentum>400"); if(mom==3) gdEdx[mom]->SetTitle("Mean dE/dx value for 400>momentum>500"); if(mom==4) gdEdx[mom]->SetTitle("Mean dE/dx value for 500>momentum>600"); if(mom==5) gdEdx[mom]->SetTitle("Mean dE/dx value for 600>momentum>700"); if(mom==6) gdEdx[mom]->SetTitle("Mean dE/dx value for 700>momentum>800"); if(mom==7) gdEdx[mom]->SetTitle("Mean dE/dx value for 800>momentum>900"); if(mom==8) gdEdx[mom]->SetTitle("Mean dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx[mom]->SetTitle("Mean dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx[mom]->SetTitle("Mean dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx[mom]->SetTitle("Mean dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdx[mom]->SetMaximum(meanval[mom]*1.3); gdEdx[mom]->SetMinimum(meanval[mom]*0.7); gdEdx[mom]->Draw("ap"); // Draw lines at approx expected mean value of dE/dx and +/- 10% TLine *line1 = new TLine(firstrun,meanval[mom],lastrun,meanval[mom]); TLine *linep = new TLine(firstrun,meanval[mom]*1.1,lastrun,meanval[mom]*1.1); TLine *linem = new TLine(firstrun,meanval[mom]*0.9,lastrun,meanval[mom]*0.9); line1->SetLineColor(2); linep->SetLineColor(2); linem->SetLineColor(2); linep->SetLineStyle(2); linem->SetLineStyle(2); line1->Draw("same"); linep->Draw("same"); linem->Draw("same"); } // fvsrun->Print(Form("dedx_run%d_run%d_cosmic.gif",firstrun,lastrun)); for(int mom=0;mom<12;mom++){ pp=fvssigma->cd(mom+1); // gdEdx_sigma[mom]->Print(); TLine *line1 = new TLine(firstrun,0.1,lastrun,0.1); line1->SetLineColor(2); line1->SetLineStyle(2); gdEdx_sigma[mom]->GetXaxis()->SetTitle("Run number"); if(mom==0) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for momentum<200"); if(mom==1) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 200>momentum>300"); if(mom==2) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 300>momentum>400"); if(mom==3) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 400>momentum>500"); if(mom==4) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 500>momentum>600"); if(mom==5) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 600>momentum>700"); if(mom==6) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 700>momentum>800"); if(mom==7) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 800>momentum>900"); if(mom==8) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdx_sigma[mom]->SetMaximum(0.2); gdEdx_sigma[mom]->SetMinimum(0.); gdEdx_sigma[mom]->Draw("ap"); line1->Draw("same"); } // fvssigma->Print(Form("dedx_sigma_run%d_run%d_cosmic.gif",firstrun,lastrun)); /* // beam triggers are too low stats, would need to compare batches of 10 runs at a time for(int mom=0;mom<12;mom++){ p=fvsrunb->cd(mom+1); // Double_t meanval[mom] = 440.; gdEdxb[mom]->GetXaxis()->SetTitle("Run number"); if(mom==0) gdEdxb[mom]->SetTitle("Mean dE/dx value for momentum<200"); if(mom==1) gdEdxb[mom]->SetTitle("Mean dE/dx value for 200>momentum>300"); if(mom==2) gdEdxb[mom]->SetTitle("Mean dE/dx value for 300>momentum>400"); if(mom==3) gdEdxb[mom]->SetTitle("Mean dE/dx value for 400>momentum>500"); if(mom==4) gdEdxb[mom]->SetTitle("Mean dE/dx value for 500>momentum>600"); if(mom==5) gdEdxb[mom]->SetTitle("Mean dE/dx value for 600>momentum>700"); if(mom==6) gdEdxb[mom]->SetTitle("Mean dE/dx value for 700>momentum>800"); if(mom==7) gdEdxb[mom]->SetTitle("Mean dE/dx value for 800>momentum>900"); if(mom==8) gdEdxb[mom]->SetTitle("Mean dE/dx value for 900>momentum>1000"); if(mom==9) gdEdxb[mom]->SetTitle("Mean dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdxb[mom]->SetTitle("Mean dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdxb[mom]->SetTitle("Mean dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdxb[mom]->SetMaximum(meanval[mom]*1.3); gdEdxb[mom]->SetMinimum(meanval[mom]*0.7); gdEdxb[mom]->Draw("ap"); // Draw lines at approx expected mean value of dE/dx and +/- 10% TLine *line1 = new TLine(firstrun,meanval[mom],lastrun,meanval[mom]); TLine *linep = new TLine(firstrun,meanval[mom]*1.1,lastrun,meanval[mom]*1.1); TLine *linem = new TLine(firstrun,meanval[mom]*0.9,lastrun,meanval[mom]*0.9); line1->SetLineColor(2); linep->SetLineColor(2); linem->SetLineColor(2); linep->SetLineStyle(2); linem->SetLineStyle(2); line1->Draw("same"); linep->Draw("same"); linem->Draw("same"); } fvsrunb->Print(Form("dedx_run%d_run%d_beam.gif",firstrun,lastrun)); for(int mom=0;mom<12;mom++){ pp=fvssigmab->cd(mom+1); TLine *line1 = new TLine(firstrun,0.1,lastrun,0.1); line1->SetLineColor(2); line1->SetLineStyle(2); gdEdxb_sigma[mom]->GetXaxis()->SetTitle("Run number"); if(mom==0) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for momentum<200"); if(mom==1) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 200>momentum>300"); if(mom==2) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 300>momentum>400"); if(mom==3) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 400>momentum>500"); if(mom==4) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 500>momentum>600"); if(mom==5) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 600>momentum>700"); if(mom==6) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 700>momentum>800"); if(mom==7) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 800>momentum>900"); if(mom==8) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 900>momentum>1000"); if(mom==9) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdxb_sigma[mom]->SetTitle("RMS dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdxb_sigma[mom]->SetMaximum(0.2); gdEdxb_sigma[mom]->SetMinimum(0.); gdEdxb_sigma[mom]->Draw("ap"); line1->Draw("same"); } fvssigmab->Print(Form("dedx_sigma_run%d_run%d_beam.gif",firstrun,lastrun)); */ for(int mom=0;mom<12;mom++){ ppp=fvsab->cd(mom+1); dedxrun[firstrun][mom]->Draw(); } } // Plot dE/dx vs. subruns for cosmic triggers // Plots mean, RMS from histogram, not fit value. void plot_dedx_subrun(int run, int firstsub, int lastsub){ double meanval[12]={600.,375.,375.,375.,375.,400.,400.,410,410.,425.,425.,500.}; // dE/dx approx mean values read off of an old plot==> would be better from a calc char title[256]; char plotcommand[256]; char cuts[256]; TH1F *dedxhist[100][12]; TGraphErrors *gdEdx[12]; TGraphErrors *gdEdx_sigma[12]; fvssub = new TCanvas("fvssub","dE/dx vs subrun, cosmic",100,100,1200,800); fvssub->Divide(3,4); fvssub->cd(1); fvssubsigma = new TCanvas("fvssub","RMS(dE/dx) vs subrun, cosmic",100,100,1200,800); fvssubsigma->Divide(3,4); fvssubsigma->cd(1); fvsab = new TCanvas("fvsab","dE/dx first run",100,100,1200,800); fvsab->Divide(3,4); fvsab->cd(1); for(int sub=firstsub;sub<=lastsub;sub++){ // open file and save a histogram for all momentum and then momentum bins TFile *pg = new TFile(Form("dq-tpc-bc_%08d_%04d.root",run,sub)); std::cout << "needs zombie protectioN!" << std::endl; if(pg->IsOpen()){ TChain *c = (TChain*)pg->Get("bc"); c->Draw("timestamp"); for(int mom=0;mom<12;mom++){ sprintf(title,"dedx_%d_%d_%d",run,sub,mom); dedxhist[sub][mom] = new TH1F(title,title,100,0.1,1000); // cint interpreter title is equivalent to the variable name Double_t lowmom =0.; Double_t highmom = 1500.; if(mom<9&&mom>0){ lowmom = 200.+100*mom; highmom = lowmom+100.; }elseif(mom==0){ lowmom = 0; highmom = 200.; }elseif(mom==9){ lowmom = 1000.; highmom = lowmom+200.; }elseif(mom==10){ lowmom = 1200.; highmom = lowmom+300.; }elseif(mom==11){ lowmom = 1500.; highmom = 2000.; } sprintf(cuts,"momentum>%d&&momentum<=%d&&ctrig",lowmom,highmom); sprintf(plotcommand,"dedx>>dedx_%d_%d_%d",run,sub,mom); c->Draw(plotcommand,cuts); } //momentum loop } // file loop } //subrun for(int mom=0;mom<12;mom++){ gdEdx[mom] = new TGraphErrors(); // Graph of event vs. (total) missed spills (need total) sprintf(title,"gdEdx_%d",mom); gdEdx[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx[mom]->SetMarkerStyle(24); gdEdx_sigma[mom] = new TGraphErrors(); sprintf(title,"gdEdx_sigma_%d",mom); gdEdx_sigma[mom]->SetName(title); // if Name is not set, will not appear in output file gdEdx_sigma[mom]->SetMarkerStyle(24); int ctr=0; for(int sub=firstsub;sub<=lastsub;sub++){ if(dedxhist[sub][mom]){ if(dedxhist[sub][mom]->GetEntries()>20.){ // only fill if >20 entries in histogram // currently about 100 events / subrun dedxhist[sub][mom]->Fit("gaus","Q","",dedxhist[sub][mom]->GetMean()-dedxhist[sub][mom]->GetRMS()*3.,dedxhist[sub][mom]->GetMean()+dedxhist[sub][mom]->GetRMS()*3.); TF1 *func = dedxhist[sub][mom]->GetFunction("gaus"); // Fitting looks terrible for subruns, so just do mean/sigma /* gdEdx[mom]->SetPoint(ctr,sub,func->GetParameter(1)); // dE/dx fitted mean gdEdx[mom]->SetPointError(ctr,0.5,func->GetParameter(2)/(TMath::Sqrt(dedxhist[sub][mom]->GetEntries()))); // fitted RMS/sqrt(entries) gdEdx_sigma[mom]->SetPoint(ctr,sub,func->GetParameter(2)/func->GetParameter(1)); // fitted RMS/ mean gdEdx_sigma[mom]->SetPointError(ctr,0.,0.); */ gdEdx[mom]->SetPoint(ctr,sub,dedxhist[sub][mom]->GetMean()); // dE/dx fitted mean gdEdx[mom]->SetPointError(ctr,0.5,dedxhist[sub][mom]->GetRMS()/(TMath::Sqrt(dedxhist[sub][mom]->GetEntries()))); // fitted RMS/sqrt(entries) gdEdx_sigma[mom]->SetPoint(ctr,sub,dedxhist[sub][mom]->GetRMS()/dedxhist[sub][mom]->GetMean()); // fitted RMS/ mean gdEdx_sigma[mom]->SetPointError(ctr,0.,0.); ctr++; } // entries > 20 } //Histo exists } //subrun loop } // momentum loop for(int mom=0;mom<12;mom++){ p=fvssub->cd(mom+1); gdEdx[mom]->GetXaxis()->SetTitle("Subrun number"); if(mom==0) gdEdx[mom]->SetTitle("Mean dE/dx value for momentum<200"); if(mom==1) gdEdx[mom]->SetTitle("Mean dE/dx value for 200>momentum>300"); if(mom==2) gdEdx[mom]->SetTitle("Mean dE/dx value for 300>momentum>400"); if(mom==3) gdEdx[mom]->SetTitle("Mean dE/dx value for 400>momentum>500"); if(mom==4) gdEdx[mom]->SetTitle("Mean dE/dx value for 500>momentum>600"); if(mom==5) gdEdx[mom]->SetTitle("Mean dE/dx value for 600>momentum>700"); if(mom==6) gdEdx[mom]->SetTitle("Mean dE/dx value for 700>momentum>800"); if(mom==7) gdEdx[mom]->SetTitle("Mean dE/dx value for 800>momentum>900"); if(mom==8) gdEdx[mom]->SetTitle("Mean dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx[mom]->SetTitle("Mean dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx[mom]->SetTitle("Mean dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx[mom]->SetTitle("Mean dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdx[mom]->SetMaximum(meanval[mom]*1.5); gdEdx[mom]->SetMinimum(meanval[mom]*0.5); gdEdx[mom]->Draw("ap"); // Draw lines at approx expected mean value of dE/dx and +/- 10% TLine *line1 = new TLine(firstsub,meanval[mom],lastsub,meanval[mom]); TLine *linep = new TLine(firstsub,meanval[mom]*1.1,lastsub,meanval[mom]*1.1); TLine *linem = new TLine(firstsub,meanval[mom]*0.9,lastsub,meanval[mom]*0.9); line1->SetLineColor(2); linep->SetLineColor(2); linem->SetLineColor(2); linep->SetLineStyle(2); linem->SetLineStyle(2); line1->Draw("same"); linep->Draw("same"); linem->Draw("same"); } fvssub->Print(Form("dedx_run%d_subrun_cosmic.gif",run)); for(int mom=0;mom<12;mom++){ pp=fvssubsigma->cd(mom+1); TLine *line1 = new TLine(firstsub,0.1,lastsub,0.1); line1->SetLineColor(2); line1->SetLineStyle(2); gdEdx_sigma[mom]->GetXaxis()->SetTitle("Run number"); if(mom==0) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for momentum<200"); if(mom==1) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 200>momentum>300"); if(mom==2) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 300>momentum>400"); if(mom==3) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 400>momentum>500"); if(mom==4) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 500>momentum>600"); if(mom==5) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 600>momentum>700"); if(mom==6) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 700>momentum>800"); if(mom==7) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 800>momentum>900"); if(mom==8) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 900>momentum>1000"); if(mom==9) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1000>momentum>1200"); if(mom==10) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1200>momentum>1500"); if(mom==11) gdEdx_sigma[mom]->SetTitle("RMS dE/dx value for 1500>momentum>2000"); // Set maximum and minimum plot range gdEdx_sigma[mom]->SetMaximum(0.3); gdEdx_sigma[mom]->SetMinimum(0.); gdEdx_sigma[mom]->Draw("ap"); line1->Draw("same"); } fvssubsigma->Print(Form("dedx_sigma_run%d_subrun_cosmic.gif",run)); for(int mom=0;mom<12;mom++){ ppp=fvsab->cd(mom+1); dedxhist[firstsub][mom]->Draw(); } } /* root [1] bc->Print() ****************************************************************************** *Tree :bc : TPC dE/dx * *Entries : 3135 : Total = 190438 bytes File Size = 80794 * * : : Tree compression factor = 1.00 * ****************************************************************************** *Br 0 :eventid : eventid/I * *Entries : 3135 : Total Size= 13182 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 1 :spillno : spillno/I * *Entries : 3135 : Total Size= 13182 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 2 :btrig : btrig/O * *Entries : 3135 : Total Size= 3759 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 3 :ctrig : ctrig/O * *Entries : 3135 : Total Size= 3759 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 4 :run : run/I * *Entries : 3135 : Total Size= 13158 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 5 :subrun : subrun/I * *Entries : 3135 : Total Size= 13176 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 6 :event : event/I * *Entries : 3135 : Total Size= 13170 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 7 :timestamp : timestamp/I * *Entries : 3135 : Total Size= 13194 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 8 :ntrack : ntrack/I * *Entries : 3135 : Total Size= 13176 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 9 :nnodes : nnodes/I * *Entries : 3135 : Total Size= 13176 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 10 :momentum : momentum/D * *Entries : 3135 : Total Size= 25736 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 11 :errmomentum : errmomentum/D * *Entries : 3135 : Total Size= 25754 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* *Br 12 :dedx : dedx/D * *Entries : 3135 : Total Size= 25712 bytes One basket in memory * *Baskets : 0 : Basket Size= 32000 bytes Compression= 1.00 * *............................................................................* */ bool FileExists(const char *name){ struct stat stFileInfo; int intStat; // Attempt to get the file attributes intStat = stat(name,&stFileInfo); if(intStat == 0) return true; else{ return false; } } void SetupStyle(){ gROOT->SetStyle("Pub"); gStyle->SetTitleFont(132,"xyzt"); gStyle->SetTitleFontSize(0.075); // gStyle->SetTitleYOffset(1.5); // doesn't work? int nbCol=36; int offset=0; int iJul; TColor **colJul; colJul = new TColor*[36]; for (iJul=0;iJul<36;iJul++) colJul[iJul] = new TColor(999-iJul,0.0,0.7,0.7); colJul[0]->SetRGB(1.0,0.9,1.0); colJul[1]->SetRGB(1.0,0.8,1.0); colJul[2]->SetRGB(1.0,0.6,1.0); colJul[3]->SetRGB(1.0,0.0,1.0); colJul[4]->SetRGB(0.88,0.0,1.0); colJul[5]->SetRGB(0.7,0.0,1.0); colJul[6]->SetRGB(0.55,0.0,1.0); colJul[7]->SetRGB(0.0,0.0,0.98); colJul[8]->SetRGB(0.0,0.5,1.0); colJul[9]->SetRGB(0.0,0.65,1.0); colJul[10]->SetRGB(0.0,0.8,1.0); colJul[11]->SetRGB(0.0,0.9,1.0); colJul[12]->SetRGB(0.0,1.0,1.0); colJul[13]->SetRGB(0.0,0.97,0.85); colJul[14]->SetRGB(0.0,0.95,0.7); colJul[15]->SetRGB(0.0,0.9,0.58); colJul[16]->SetRGB(0.0,0.87,0.0); colJul[17]->SetRGB(0.5,0.9,0.0); colJul[18]->SetRGB(0.65,0.93,0.0); colJul[19]->SetRGB(0.80,0.96,0.0); colJul[20]->SetRGB(0.9,0.98,0.0); colJul[21]->SetRGB(1.0,1.0,0.0); colJul[22]->SetRGB(1.0,0.95,0.0); colJul[23]->SetRGB(1.0,0.90,0.0); colJul[24]->SetRGB(1.0,0.85,0.0); colJul[25]->SetRGB(1.0,0.78,0.0); colJul[26]->SetRGB(1.0,0.68,0.0); colJul[27]->SetRGB(1.0,0.58,0.0); colJul[28]->SetRGB(1.0,0.48,0.0); colJul[29]->SetRGB(1.0,0.0,0.0); colJul[30]->SetRGB(0.92,0.0,0.0); colJul[31]->SetRGB(0.85,0.0,0.0); colJul[32]->SetRGB(0.74,0.0,0.0); colJul[33]->SetRGB(0.62,0.0,0.0); colJul[34]->SetRGB(0.5,0.0,0.0); colJul[35]->SetRGB(0.0,0.0,0.0); int tabColJul[36]; for (iJul=0;iJul<36;iJul++) tabColJul[iJul]=999-iJul; gStyle->SetPalette(nbCol,tabColJul+offset); printf("\nCAUTION : defining colors 964 to 999.\n"); printf("CAUTION : changing palette ; now has %d colors.\n\n",nbCol); return; }