StRoot: StEmcPool/StEmcOfflineCalibrationMaker/macros/electron

00001 #include <iostream>
00002 #include <fstream>
00003 using namespace std;
00004 
00005 #include "TFile.h"
00006 #include "TH1.h"
00007 #include "TH2.h"
00008 #include "TF1.h"
00009 #include "TMath.h"
00010 #include "TPostScript.h"
00011 #include "TCanvas.h"
00012 #include "TStyle.h"
00013 #include "TLatex.h"
00014 #include "TTree.h"
00015 #include "TLine.h"
00016 
00017 #include "CalibrationHelperFunctions.h"
00018 #include "CalibrationHelperFunctions.cxx"
00019 
00020 void drawTower(TH1* h, TF1* f, int id, CalibrationHelperFunctions* helper);
00021 double fit_function(double *x, double *par);
00022 double fit_function2(double *x, double *par);
00023 double background_only_fit(double *x, double *par);
00024 int lookup_crate(float,float);
00025 
00026 void electron_master(const char* file_list="electrons.list",const char* output="electronmaster.root", const char* dbDate="1999-01-01 00:11:00", const char* geantfile="geant_func.root", const char* gfname="mip.gains.final", const char* ngname="electron.gains", const char* rgname="matrel.gains"){
00027   //**********************************************//
00028   //Load Libraries                                //
00029   //**********************************************//
00030   gROOT->Macro("LoadLogger.C");
00031   gROOT->Macro("loadMuDst.C");
00032   gSystem->Load("StTpcDb");
00033   gSystem->Load("StDaqLib");
00034   gSystem->Load("StDetectorDbMaker");
00035   gSystem->Load("St_db_Maker");
00036   gSystem->Load("StDbUtilities");
00037   gSystem->Load("StEmcRawMaker");
00038   gSystem->Load("StMcEvent");
00039   gSystem->Load("StMcEventMaker");//***
00040   gSystem->Load("StEmcSimulatorMaker");//***
00041   gSystem->Load("StEmcADCtoEMaker");
00042   gSystem->Load("StEpcMaker");
00043   gSystem->Load("StDbBroker");
00044   gSystem->Load("StEEmcUtil");
00045   gSystem->Load("StAssociationMaker");
00046   gSystem->Load("StEmcTriggerMaker");
00047   gSystem->Load("StTriggerUtilities");
00048   gSystem->Load("StEmcOfflineCalibrationMaker");
00049 
00050   cout<<"input filelist: "<<file_list<<endl;
00051   cout<<"output:          "<<output<<endl;
00052   cout<<"db Date:        "<<dbDate<<endl;
00053   cout<<"geant curve:    "<<geantfile<<endl;
00054 
00055 
00056   //**********************************************//
00057   //Initialize Stuff                              //
00058   //**********************************************//
00059 
00060   CalibrationHelperFunctions *helper = new CalibrationHelperFunctions();
00061         
00062   StBemcTablesWriter *bemctables = new StBemcTablesWriter();
00063   bemctables->loadTables(dbDate,"sim");
00064 
00065   //chain all input files together
00066   char file[300];
00067   TChain* tree = new TChain("skimTree");
00068   ifstream filelist(file_list);
00069   while(1){
00070     filelist >> file;
00071     if(!filelist.good()) break;
00072     cout<<file<<endl;
00073     tree->Add(file);
00074   }
00075 
00076   const int ntowers = 4800;
00077   const int nrings = 40;
00078   const int ncrates = 30;
00079   float gains[ntowers];
00080   int status[ntowers];
00081   float peaks[ntowers];
00082   float peakerr[ntowers];
00083   float gainerr[ntowers];
00084 
00085   
00086   ifstream rgfile(rgname);
00087   while(1){
00088     int id;
00089     float rg,rgerr;
00090     rgfile >> id >> rg >> rgerr;
00091     if(!rgfile.good())break;
00092     gains[id-1] = 0;
00093     if(rg > 0)gains[id-1] = rg;
00094     gainerr[id-1] = rgerr;
00095     status[id-1] = 0;
00096     if(rg != 0 && rgerr < 1)status[id-1] = 1;
00097   }
00098 
00099   ifstream gainfile(gfname);
00100   while(1){
00101     int id,stat;
00102     float peak,err,gain,eta,theta;
00103     gainfile >> id >> peak >> err >> stat;
00104     if(!gainfile.good())break;
00105     //eta = helper->getEta(id);
00106     //theta = helper->getTheta(id);
00107     //gain = 0.264*(1+0.056*eta*eta)/(sin(theta)*peak);
00108     //peaks[id-1] = peak;
00109     //gains[id-1] = gain;
00110     if(status[id-1] != 1 || stat != 1) status[id-1] = stat;
00111     //peakerr[id-1] = err;
00112     //cout<<id<<" "<<gain;
00113   }
00114 
00115 
00116   
00117   TFile* geant_file = new TFile(geantfile,"READ");
00118   TF1* geant_fit = (TF1*)geant_file->Get("geant_fit");
00119 
00120   TFile outfile(output,"RECREATE");
00121 
00122   TH1 *crate_histo[ncrates];
00123   TF1 *crate_fit[ncrates];
00124   TH1 *electron_histo[ntowers];
00125   TF1 *fit[ntowers];
00126   TH1 *prs_histo[ntowers];
00127   TH1 *ring_histo[nrings];
00128   TH1 *ring2_histo[nrings/2];
00129   TF1 *ring2fit[nrings/2];
00130 
00131   TF1 *ringfit[nrings];
00132   TH2F* ring_pve[nrings];
00133   TH2F* jan_pve[6];
00134   float eta;
00135   int etaindex;
00136   TH1F* ringprec = new TH1F("ringprec","",40,-1.0,1.0);
00137   TH1F* ringprec2 = new TH1F("ringprec2","",40,-1.0,1.0);
00138   ringprec->SetYTitle("E/p");
00139   ringprec->SetXTitle("#eta");
00140   ringprec2->SetYTitle("E/p");
00141   ringprec2->SetXTitle("#eta");
00142   TH1F* ring2prec = new TH1F("ring2prec","",nrings/2,-1.0,1.0);
00143   TH1F* ring2prec2 = new TH1F("ring2prec2","",nrings/2,-1.0,1.0);
00144   ring2prec->SetYTitle("E/p");
00145   ring2prec->SetXTitle("#eta");
00146   ring2prec2->SetYTitle("E/p");
00147   ring2prec2->SetXTitle("#eta");
00148   double ew[nrings];
00149 
00150   TH1F* crateprec = new TH1F("crateprec","",30,0.5,30.5);
00151   crateprec->SetXTitle("Crate");
00152   crateprec->SetYTitle("E/p");
00153 
00154   TH2F *energyleak = new TH2F("energyleak","",20,0.0,0.03,20,0.0,1.0);
00155   TH2F *findbg = new TH2F("findbg","",20,0.0,0.03,30,0.0,5.0);
00156   TH1F *energymean = new TH1F("energymean","",20,0.0,0.03);
00157   TH1F *leakmean = new TH1F("leakmean","",20,0.0,0.03);
00158   energyleak->SetXTitle("#DeltaR");
00159   energyleak->SetYTitle("leaked energy / total energy");
00160   findbg->SetXTitle("#DeltaR");
00161   findbg->SetYTitle("Total energy / track momentum");
00162   TH2F *tevsp = new TH2F("tevsp","",50,0.0,20.0,50,0.0,30.0);
00163   TH1F *pmean = new TH1F("pmean","",20,0.0,15.0);
00164   tevsp->SetXTitle("track momentum (GeV)");
00165   tevsp->SetYTitle("Total Energy (GeV)");
00166   TH2F *tevspcent = new TH2F("tevspcent","",20,0.0,15.0,20,0.0,15.0);
00167   tevspcent->SetXTitle("track momentum (GeV)");
00168   tevspcent->SetYTitle("Energy in central tower (GeV)");
00169   TH1F *cmean = new TH1F("cmain","",20,0.0,15.0);
00170   TH2F *sistertracks = new TH2F("sistertracks","",20,0.0,8.0,20,0.0,8.0);
00171   sistertracks->SetXTitle("Track momentum (GeV)");
00172   sistertracks->SetYTitle("Neighbor momentum (GeV)");
00173   TH2F* dEdxvsp = new TH2F("dEdxvsp","",100,0.15,1.3,100,-5.7,-5.);
00174   dEdxvsp->SetXTitle("Log(p)");
00175   dEdxvsp->SetYTitle("Log(dE/dx)");
00176   TH2F* dEdxvsp_east = new TH2F("dEdxvsp_east","",100,0.15,1.3,100,-5.7,-5.0);
00177   dEdxvsp_east->SetXTitle("Log(p)");
00178   dEdxvsp_east->SetYTitle("Log(dE/dx)");
00179   TH2F* dEdxvsp_west = new TH2F("dEdxvsp_west","",100,0.15,1.3,100,-5.7,-5.0);
00180   dEdxvsp_west->SetXTitle("Log(p)");
00181   dEdxvsp_west->SetYTitle("Log(dE/dx)");
00182   TH2F* energyleak2 = new TH2F("energyleak2","",20,0.0,0.03,20,0.0,1.0);
00183   TH1F* energymean2 = new TH1F("energymean2","",20,0.0,0.03);
00184   TH1F* towermult = new TH1F("towermult","",9,0.0,9.0);
00185   energyleak2->SetXTitle("#DeltaR");
00186   energyleak2->SetYTitle("leaked energy/total energy");
00187   towermult->SetXTitle("Neighbors with energy");
00188   TH2F* multvsp = new TH2F("multvsp","",20,0.0,20.0,9,0.0,9.0);
00189   multvsp->SetXTitle("Track momentum (GeV)");
00190   multvsp->SetYTitle("Neighbors with energy");
00191   TH1F* multmean = new TH1F("multmean","",20,0.0,20.0);
00192   TH2F* tep3 = new TH2F("tep3","2 < p < 3",20,0.0,0.03,20,0.0,4.0);
00193   TH2F* tep5 = new TH2F("tep5","3 < p < 5",20,0.0,0.03,20,0.0,4.0);
00194   TH2F* tep10 = new TH2F("tep10","5 < p < 10",20,0.0,0.03,20,0.0,4.0);
00195   tep3->SetXTitle("DeltaR");
00196   tep3->SetYTitle("Total energy / track momentum");
00197   tep5->SetXTitle("DeltaR");
00198   tep5->SetYTitle("Total energy / track momentum");
00199   tep10->SetXTitle("DeltaR");
00200   tep10->SetYTitle("Total energy / track momentum");
00201   TH1F* tep3mean = new TH1F("tep3mean","",20,0,0.03);
00202   TH1F* tep5mean = new TH1F("tep5mean","",20,0,0.03);
00203   TH1F* tep10mean = new TH1F("tep10mean","",20,0,0.03);
00204   TH1F* multen = new TH1F("multen","",40,0.0,1.0);
00205   multen->SetXTitle("Energy in neighbor towers (GeV)");
00206   TH1F* east_histo = new TH1F("east_histo","Electron E/p in East",40,0.0,2.0);
00207   TH1F* west_histo = new TH1F("west_histo","Electron E/p in West",40,0.0,2.0);
00208   TH1F* all_histo = new TH1F("all_histo","Electron E/p",40,0.0,2.0);
00209   TH2F* pvsep = new TH2F("pvsep","Electron p vs E/p",120,0,3.0,20,0,20.0);
00210   pvsep->SetYTitle("p (Gev)");
00211   pvsep->SetXTitle("E/p");
00212   TH2F* pvsep0 = new TH2F("pvsep0","Electron p vs E/p",120,0,3.0,20,0,20.0);
00213   pvsep0->SetYTitle("p (Gev)");
00214   pvsep0->SetXTitle("E/p");
00215   TH2F* evsep = new TH2F("evsep","Electron E vs E/p",120,0,3.0,20,0,20.0);
00216   evsep->SetYTitle("E (GeV)");
00217   evsep->SetXTitle("E/p");
00218 
00219   TH1F* bsmde = new TH1F("bsmde","BSMDE ADC TOT",1500,0,1500);
00220   TH1F* bsmdp = new TH1F("bsmdp","BSMDP ADC TOT",1500,0,1500);
00221 
00222   TH1F* bsmde_central = new TH1F("bsmde_central","BSMDE ADC TOT",100,0,1500);
00223   TH1F* bsmde_mid = new TH1F("bsmde_mid","BSMDE ADC TOT",100,0,1500);
00224   TH1F* bsmde_outer = new TH1F("bsmde_outer","BSMDE ADC TOT",100,0,1500);
00225 
00226   TH2F* bsmdep = new TH2F("bsmdep","BSMDE v BSMDP",100,0,1500,100,0,1500);
00227 
00228   TH2F* bsmdevp = new TH2F("bsmdevp","BSMDE v p",100,1.0,15.0,100,0,1500);
00229   TH2F* bsmdpvp = new TH2F("bsmdpvp","BSMDP v p",100,1.0,15.0,100,0,1500);
00230   TH2F* bsmdevep = new TH2F("bsmdevep","BSMDE v E/p",100,0.0,2.0,100,0,1500);
00231   TH2F* bsmdpvep = new TH2F("bsmdpvep","BSMDP v E/p",100,0.0,2.0,100,0,1500);
00232 
00233   TH2F* bsmdeve = new TH2F("bsmdeve","BSMDE v E",100,1.0,30.0,100,0,1500);
00234   TH2F* bsmdpve = new TH2F("bsmdpve","BSMDP v E",100,1.0,30.0,100,0,1500);
00235 
00236   TH1F* httrig = new TH1F("httrig","HT Trigger",5,-0.5,4.5);
00237 
00238   TH1F* pplus = new TH1F("pplus","e+ p",100,0,20);
00239   TH1F* pminus = new TH1F("pminus","e- p",100,0,20);
00240 
00241   TH1F* posep = new TH1F("posep","e+ E/p",60,0,3.0);
00242   TH1F* negep = new TH1F("negep","e- E/p",60,0,3.0);
00243 
00244 
00245   double mipgains[40] = {0.0153384,0.0199435,0.0213435,0.0195892,0.0200876,0.0182992,0.0178393,0.0174514,0.0174794,0.0164953,0.0158151,0.0166528,0.0162374,0.0146685,0.0149648,0.0153254,0.0146326,0.0151221,0.0145512,0.0148793,0.0150701,0.0145633,0.0150305,0.0154839,0.0142929,0.0147277,0.0150759,0.0158363,0.0156365,0.0160211,0.0164128,0.0172054,0.017725,0.0172588,0.0177851,0.0187099,0.0200016,0.0200538,0.0169604,0.014336};
00246 
00247 
00248    //create the tower histograms
00249   char name1[100];
00250   /*
00251   for(int i=0; i<ntowers; i++){
00252     sprintf(name,"electron_histo_%i",i+1);
00253     electron_histo[i] = new TH1D(name,"",30,0.,140.);
00254     electron_histo[i]->SetXTitle("#frac{ADC}{GeV * Sin(#theta)}");
00255     sprintf(name1,"prs_histo_%i",i+1);
00256     prs_histo[i] = new TH1D(name1,"",60,0.,500.);
00257     prs_histo[i]->SetXTitle("ADC");
00258   }
00259   */
00260   for(int i = 0; i < ncrates; i++){
00261     char name[100];
00262     char title[100];
00263     sprintf(name,"crate_%i",i+1);
00264     sprintf(title,"E/p for Crate %i",i+1);
00265     crate_histo[i] = new TH1F(name,title,60,0.,3.0);
00266     crate_histo[i]->SetXTitle("E/p");
00267   }
00268 
00269   for(int i = 0; i < nrings/2; i++){
00270     char name[100];
00271     sprintf(name,"ring2_histo_%i",i);
00272     ring2_histo[i] = new TH1F(name,"",60,0.,3.0);
00273     ring2_histo[i]->SetXTitle("E/p");
00274   }
00275 
00276   char name2[100];
00277   for(int i=0; i<nrings;i++)
00278     {
00279       sprintf(name2,"ring_histo_%i",i);
00280       //ring_histo[i] = new TH1D(name2,"",30,0.,140.0);
00281       //ring_histo[i]->SetXTitle("ADC / GeV Sin(#theta)");
00282       ring_histo[i] = new TH1D(name2,"",60,0.,3.0);
00283       ring_histo[i]->SetXTitle("E/p");
00284       char namerpve[100];
00285       sprintf(namerpve,"ring_pve_%i",i);
00286       ring_pve[i] = new TH2F(namerpve,"",20,0,20.0,20,0,20.0);
00287       ring_pve[i]->SetXTitle("E (GeV)");
00288       ring_pve[i]->SetYTitle("p (GeV)");
00289     }
00290 
00291   for(int i = 0; i < 6; i++){
00292     char jname[100];
00293     sprintf(jname,"jan_pve_%i",i);
00294     jan_pve[i] = new TH2F(jname,"",120,0,3.0,20,0,20.0);
00295     jan_pve[i]->SetXTitle("E/p");
00296     jan_pve[i]->SetYTitle("p (GeV)");
00297   }
00298   //global graphics functions
00299   gStyle->SetOptStat("oue");
00300   gStyle->SetOptFit(111);
00301   gStyle->SetCanvasColor(10);
00302   gStyle->SetCanvasBorderMode(0);
00303   gStyle->SetOptTitle(0);
00304   gStyle->SetPalette(1);
00305   gStyle->SetStatColor(0);
00306 
00307   StEmcOfflineCalibrationTrack* track = new StEmcOfflineCalibrationTrack();
00308   StEmcOfflineCalibrationCluster* cluster = new StEmcOfflineCalibrationCluster();
00309   tree->SetBranchAddress("clusters",&cluster);
00310 
00311   //**********************************************//
00312   //Loop Over Tracks, Fill Histograms             //
00313   //**********************************************//
00314 
00315   int nentries = tree->GetEntries();
00316   cout<<nentries<<endl;
00317   int ngoodhit = 0;
00318   int nplt10 = 0;
00319   int nnosis = 0;
00320   int nfinal = 0;
00321   int nbsmdgood = 0;
00322   int nnottrig = 0;
00323   int nfidu = 0;
00324   int nenterexit = 0;
00325   for(int j=0; j<nentries; j++){
00326     tree->GetEntry(j);
00327     track = &(cluster->centralTrack);
00328     TClonesArray *tracks = cluster->tracks;
00329     //if(1) cout<<"reading "<<j<<" of "<<nentries<<endl;
00330 
00331     httrig->Fill((float)track->htTrig);
00332 
00333     if(track->charge > 0)pplus->Fill(track->p);
00334     if(track->charge < 0)pminus->Fill(track->p);
00335 
00336     int bsmdeadctot = 0;
00337     int bsmdpadctot = 0;
00338     for(int i = 0; i < 11; i++){
00339       if(track->smde_adc[i] > track->smde_pedestal[i])bsmdeadctot += track->smde_adc[i] - track->smde_pedestal[i];
00340       if(track->smdp_adc[i] > track->smdp_pedestal[i])bsmdpadctot += track->smdp_adc[i] - track->smdp_pedestal[i];
00341     }
00342 
00343     double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);
00344 
00345 
00346     double scaled_adc = (track->tower_adc[0] - track->tower_pedestal[0]) / track->p;
00347     double geant_scale = geant_fit->Eval(dR);
00348     scaled_adc *= geant_scale;
00349     //cout<<scaled_adc<<endl;
00350     int index = track->tower_id[0];
00351 
00352     //figure out eta and etaindex
00353     eta = helper->getEta(index);
00354     if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
00355     etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
00356 
00357     double tgain = bemctables->calib(1,track->tower_id[0])*gains[index-1];
00358     //double tgain = gains[index-1];
00359 
00360 
00361     if((track->tower_adc[0] - track->tower_pedestal[0]) < 2.5 * track->tower_pedestal_rms[0])continue;
00362     ngoodhit++;
00363     dEdxvsp->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
00364     if(track->tower_id[0] <= 2400)dEdxvsp_west->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
00365     if(track->tower_id[0] > 2400)dEdxvsp_east->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
00366 
00367     if(track->tower_id[0] != track->tower_id_exit)continue;
00368     nenterexit++;
00369 
00370     pvsep0->Fill(scaled_adc*tgain,track->p);
00371     if(track->p > 6)continue;   
00372     nplt10++;           
00373 
00374     //if(track->p < 1.5)continue;
00375     //if(track->p < 3.0)continue;
00376     if(track->htTrig == 2)continue;
00377     nnottrig++;
00378     //change the fiducial cut to a square with diagonal = 0.06 in deta, dphi space
00379     float squarefid = 0.02;//0.03/TMath::Sqrt(2.0);
00380     //if(TMath::Abs(track->deta) > squarefid || TMath::Abs(track->dphi) > squarefid)continue;
00381 
00382     //calculate geant scaled, pedestal subtracted adc
00383     //if(dR > 0.0125)continue;
00384     if(dR > 0.025)continue;
00385     nfidu++;
00386     //if(track->p > 6.0)continue;                               
00387     //if(track->p > 15)continue;
00388     //cout<<track->dEdx<<endl;
00389 
00390 
00391     //if(track->dEdx*1000000 > 4.5 || track->dEdx*1000000 < 3.5)continue;
00392 
00393     //cout<<track->htTrig<<endl;
00394 
00395 
00396     if(track->dEdx < 3.4e-6)continue;
00397     if((bsmdeadctot > -1 && bsmdeadctot < 50) && (bsmdpadctot < 50 && bsmdpadctot > -1))continue;
00398     nbsmdgood++;
00399     //if(bsmdeadctot < 500) continue;
00400 
00401     //if(bsmdeadctot < 84.*track->p)continue;
00402     //if(bsmdeadctot > 200.*track->p + 1500)continue;
00403 
00404     //if(bsmdpadctot < 800)continue;
00405 
00406     int numsis = tracks->GetEntries();
00407     if(numsis > 0)continue;
00408 
00409     float totalbtow = 0;
00410     float maxEt = 0;
00411     int maxId = -1;
00412     for(int i = 0; i < 9; i++){
00413       if(track->tower_adc[i] - track->tower_pedestal[i] < 0)continue;
00414       float theta = helper->getTheta(track->tower_id[i]);
00415       float nextEt = (track->tower_adc[i] - track->tower_pedestal[i]) * bemctables->calib(1,track->tower_id[i])*sin(theta);
00416       totalbtow += nextEt;
00417       if(nextEt > maxEt){
00418         maxEt = nextEt;
00419         maxId = i;
00420       }
00421     }
00422     nnosis++;
00423     if(maxId != 0) continue;
00424 
00425 
00426     nfinal++;
00427 
00428     //if(!track->nonhtTrig)continue;
00429     //if(track->nHits < 25)continue;
00430     //if(status[index-1]==1)ring_histo[etaindex]->Fill(scaled_adc*gains[index-1]);
00431 
00432     //scaled_adc = totalbtow/track->p;
00433     //tgain = 1.0;
00434 
00435     ring_histo[etaindex]->Fill(scaled_adc*tgain);
00436     if(etaindex == 0 || etaindex == 39){
00437       //cout<<etaindex<<" "<<tgain<<" "<<track->p<<" "<<scaled_adc*tgain*track->p<<" "<<scaled_adc*tgain<<endl;
00438     }
00439     //cout<<index<<" "<<gains[index-1]<<" "<<scaled_adc*track->p<<" "<<track->p<<" "<<status[index-1]<<endl;
00440     ring_pve[etaindex]->Fill(scaled_adc*tgain*track->p,track->p);
00441 
00442     dEdxvsp->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
00443 
00444     float abseta = TMath::Abs(eta);
00445     if(abseta > 0.95){
00446       jan_pve[5]->Fill(scaled_adc*tgain,track->p);
00447     }else if(abseta > 0.9){
00448       jan_pve[4]->Fill(scaled_adc*tgain,track->p);
00449     }else if(abseta > 0.6){
00450       jan_pve[3]->Fill(scaled_adc*tgain,track->p);
00451     }else if(abseta > 0.3){
00452       jan_pve[2]->Fill(scaled_adc*tgain,track->p);
00453     }else if(abseta > 0.05){
00454       jan_pve[1]->Fill(scaled_adc*tgain,track->p);
00455     }else{
00456       jan_pve[0]->Fill(scaled_adc*tgain,track->p);
00457     }
00458 
00459     all_histo->Fill(scaled_adc*tgain);
00460     pvsep->Fill(scaled_adc*tgain,track->p);
00461     evsep->Fill(scaled_adc*tgain,track->p*scaled_adc*tgain);
00462     tevsp->Fill(track->p,scaled_adc*tgain*track->p);
00463 
00464     if(track->charge > 0)posep->Fill(scaled_adc*tgain);
00465     if(track->charge < 0)negep->Fill(scaled_adc*tgain);
00466 
00467     //if(scaled_adc*tgain < 0.7 || scaled_adc*tgain > 5.0)continue;
00468     bsmde->Fill(bsmdeadctot);
00469     bsmdp->Fill(bsmdpadctot);
00470     bsmdep->Fill(bsmdeadctot,bsmdpadctot);
00471 
00472     if(abseta > 0.6){
00473       bsmde_outer->Fill(bsmdeadctot);
00474     }else if(abseta > 0.3){
00475       bsmde_mid->Fill(bsmdeadctot);
00476     }else{
00477       bsmde_central->Fill(bsmdeadctot);
00478     }
00479 
00480     bsmdevp->Fill(track->p,bsmdeadctot);
00481     bsmdpvp->Fill(track->p,bsmdpadctot);
00482     bsmdevep->Fill(scaled_adc*tgain,bsmdeadctot);
00483     bsmdpvep->Fill(scaled_adc*tgain,bsmdpadctot);
00484 
00485     bsmdeve->Fill(scaled_adc*tgain*track->p,bsmdeadctot);
00486     bsmdpve->Fill(scaled_adc*tgain*track->p,bsmdpadctot);
00487 
00488     if(dR > 0.015)continue;
00489     if(track->tower_id[0] <= 2400){
00490       west_histo->Fill(scaled_adc*tgain);
00491     }
00492     if(track->tower_id[0] > 2400){
00493       east_histo->Fill(scaled_adc*tgain);
00494     }
00495 
00496   }
00497   cout<<"processed electron tree"<<endl;
00498   cout<<"ngoodhit: "<<ngoodhit<<endl;
00499   cout<<"nenterexit: "<<nenterexit<<endl;
00500   cout<<"n not trig: "<<nnottrig<<endl;
00501   cout<<"n p < 10: "<<nplt10<<endl;
00502   cout<<"n in fidu: "<<nfidu<<endl;
00503   cout<<"nbsmdhit: "<<nbsmdgood<<endl;
00504   cout<<"n no sis: "<<nnosis<<endl;
00505   cout<<"n final: "<<nfinal<<endl;
00506 
00507   double ew[21];
00508   for(int h=0;h<21;h++){
00509     TH1D* projection = energyleak->ProjectionY("projection",h,h);
00510     float mean = projection->GetMean();
00511     energymean->SetBinContent(h,mean);
00512     TH1D* projection1 = findbg->ProjectionY("projection1",h,h);
00513     float mean1 = projection1->GetMean();
00514     leakmean->SetBinContent(h,mean1);
00515     TH1D* projection2 = tevsp->ProjectionY("projection2",h,h);
00516     float mean2 = projection2->GetMean();
00517     pmean->SetBinContent(h,mean2);
00518     ew[h] = projection2->GetRMS();
00519     TH1D* projection3 = tevspcent->ProjectionY("projection3",h,h);
00520     float mean3 = projection3->GetMean();
00521     cmean->SetBinContent(h,mean3);
00522     TH1D* projection4 = energyleak2->ProjectionY("projection4",h,h);
00523     float mean4 = projection4->GetMean();
00524     energymean2->SetBinContent(h,mean4);
00525     TH1D* projection5 = multvsp->ProjectionY("projection5",h,h);
00526     float mean5 = projection5->GetMean();
00527     multmean->SetBinContent(h,mean5);
00528     TH1D* projection6 = tep3->ProjectionY("projection6",h,h);
00529     float mean6 = projection6->GetMean();
00530     tep3mean->SetBinContent(h,mean6);
00531     TH1D* projection7 = tep3->ProjectionY("projection7",h,h);
00532     float mean7 = projection7->GetMean();
00533     tep5mean->SetBinContent(h,mean7);
00534     TH1D* projection8 = tep3->ProjectionY("projection8",h,h);
00535     float mean8 = projection8->GetMean();
00536     tep10mean->SetBinContent(h,mean8);
00537   }
00538 
00539   TF1* fitleak = new TF1("fitleak","[0]",0,0.03);
00540   fitleak->SetLineWidth(0.1);
00541   leakmean->Fit(fitleak,"rq");
00542 
00543 
00544   //**********************************************//
00545   //Fit Tower Histograms                          //
00546   //**********************************************//
00547   /*
00548   for(int i=0; i<ntowers; i++){
00549           
00550     if(i%600 == 0) cout<<"fitting tower "<<i+1<<" of "<<ntowers<<endl;
00551                 
00552     sprintf(name,"fit_%i",i+1);
00553                 
00554                 //this fit is for the electron tree
00555     fit[i] = new TF1(name,fit_function,0.,140.,6);
00556     fit[i]->SetParameter(1,65.);
00557     fit[i]->SetParameter(2,10.);
00558     fit[i]->SetParameter(3,10.); //relative height of peak to bg
00559     fit[i]->SetParameter(4,10.);
00560     fit[i]->SetParameter(5,3.);
00561     fit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma");
00562                 
00563     fit[i]->SetLineColor(kGreen);
00564     fit[i]->SetLineWidth(0.6);
00565                 
00566     electron_histo[i]->Fit(fit[i],"rq");
00567   }
00568   */
00569   //**********************************************//
00570   //Fit Ring Histograms                           //
00571   //**********************************************//
00572 
00573   for(int i=0; i<nrings; i++){
00574           
00575     //cout<<"fitting ring "<<i+1<<" of "<<nrings<<endl;
00576                 
00577     sprintf(name,"ring_fit_%i",i);
00578     /*
00579     ringfit[i] = new TF1(name,fit_function,0.,140.,6);
00580     ringfit[i]->SetParameter(1,1.);
00581     ringfit[i]->SetParameter(2,0.2);
00582     ringfit[i]->SetParameter(3,1.5); //relative height of peak to bg
00583     ringfit[i]->SetParameter(4,0.15);
00584     ringfit[i]->SetParameter(5,0.8);
00585     ringfit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma");
00586     */          
00587     //TF1* ffff = new TF1("ffff","expo(0) + gaus(2)",0.4,1.7);
00588 
00589     ring_histo[i]->Sumw2();
00590     //ring_histo[i]->Rebin(3);
00591 
00592     ringfit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);
00593     ringfit[i]->SetParLimits(0,0,10.0*ring_histo[i]->GetBinContent(1));
00594     ringfit[i]->SetParLimits(1,-10000,0);
00595     ringfit[i]->SetParLimits(2,0,10.0*ring_histo[i]->GetMaximum());
00596     ringfit[i]->SetParLimits(3,0,10);
00597     ringfit[i]->SetParameter(0,ring_histo[i]->GetBinContent(1));
00598     ringfit[i]->SetParameter(1,-ring_histo[i]->GetBinContent(1)/6.0);
00599     ringfit[i]->SetParameter(2,ring_histo[i]->GetMaximum());
00600     ringfit[i]->SetParameter(3,0.95);
00601     ringfit[i]->SetParameter(4,0.15);
00602     ringfit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
00603     /*
00604     ringfit[i] = new TF1(name,fit_function2,0.1,2.5,8);
00605     ringfit[i]->SetParameter(1,1.);
00606     ringfit[i]->SetParameter(2,0.2);
00607     ringfit[i]->SetParameter(3,1.5); //relative height of peak to bg
00608     ringfit[i]->SetParameter(4,0.25);
00609     ringfit[i]->SetParameter(5,0.15);
00610     ringfit[i]->SetParameter(7,0.8);
00611     ringfit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma","Bg2 constant","Bg2 decay");
00612         
00613     */
00614     ringfit[i]->SetLineColor(kBlue);
00615     ringfit[i]->SetLineWidth(0.6);
00616     
00617     ring_histo[i]->Fit(ringfit[i],"rql","",0.3,1.7);
00618                 
00619     ringprec->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
00620     ringprec->SetBinError(i+1,ringfit[i]->GetParameter(4));
00621     ringprec2->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
00622     ringprec2->SetBinError(i+1,ringfit[i]->GetParError(3));
00623     //ew[i] = 4066/(60*(fit[i]->GetParameter(2))*(fit[i]->GetParameter(2)));
00624 
00625     float mean = ringfit[i]->GetParameter(3);
00626     float merr = ringfit[i]->GetParError(3);
00627     cout<<"ring "<<i<<" "<<mean<<" "<<merr/mean<<" "<<ring_histo[i]->GetEntries()<<endl;
00628   }
00629 
00630   for(int i = 0; i < nrings/2; i++){
00631     ring2_histo[i]->Add(ring_histo[2*i]);
00632     ring2_histo[i]->Add(ring_histo[2*i+1]);
00633     sprintf(name,"ring2_fit_%i",i);
00634     ring2fit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);
00635 
00636     ring2_histo[i]->Rebin(3);
00637 
00638     ring2fit[i]->SetParLimits(0,0,10.0*ring2_histo[i]->GetBinContent(1));
00639     ring2fit[i]->SetParLimits(1,-10000,0);
00640     ring2fit[i]->SetParLimits(2,0,10.0*ring2_histo[i]->GetMaximum());
00641     ring2fit[i]->SetParLimits(3,0,10);
00642     ring2fit[i]->SetParLimits(4,0.17,0.175);
00643     ring2fit[i]->SetParameter(0,ring2_histo[i]->GetBinContent(1));
00644     ring2fit[i]->SetParameter(1,-ring2_histo[i]->GetBinContent(1)/6.0);
00645     ring2fit[i]->SetParameter(2,ring2_histo[i]->GetMaximum());
00646     ring2fit[i]->SetParameter(3,0.95);
00647     ring2fit[i]->SetParameter(4,0.11245);
00648     ring2fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
00649 
00650     ring2_histo[i]->Fit(ring2fit[i],"rql","",0.3,1.7);
00651 
00652     ring2prec->SetBinContent(i+1,(ring2fit[i]->GetParameter(3)));
00653     ring2prec->SetBinError(i+1,ring2fit[i]->GetParameter(4));
00654     ring2prec2->SetBinContent(i+1,(ring2fit[i]->GetParameter(3)));
00655     ring2prec2->SetBinError(i+1,ring2fit[i]->GetParError(3));
00656 
00657     cout<<"ring2 "<<i<<" "<<ring2fit[i]->GetParameter(3)<<" "<<ring2fit[i]->GetParError(3)<<endl;
00658   }
00659 
00660   ofstream newgain(ngname);
00661   
00662 
00663   float gains2[ntowers];
00664   float gerr2[ntowers];
00665   for(int i = 0; i < ntowers; i++){
00666     float eta = helper->getEta(i+1);
00667     if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
00668     int etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
00669     float adjust = ring2fit[(int)etaindex/2]->GetParameter(3);
00670     //cout<<etaindex<<" "<<(int)etaindex/2<<" "<<adjust<<endl;
00671     float ng = 0;
00672     float ne = 0;
00673     if(status[i] == 1){
00674       float og = bemctables->calib(1,i+1)*gains[i];
00675       ng = og/adjust;
00676       float aerr = ring2fit[(int)etaindex/2]->GetParError(3);
00677       ne = sqrt(pow(og*aerr/(adjust*adjust),2));
00678     }
00679     newgain << i+1 << " " << ng << " " << ne << " " << status[i] << endl;
00680     gains2[i] = ng;
00681     gerr2[i] = ne;
00682 
00683 
00684   }
00685 
00686   newgain.close();
00687 
00689    //Using new gains regenerate by crate//
00691 
00692 
00693   /*
00694   for(int j=0; j<nentries; j++){
00695     tree->GetEntry(j);
00696     track = &(cluster->centralTrack);
00697     TClonesArray *tracks = cluster->tracks;
00698 
00699     int bsmdeadctot = 0;
00700     int bsmdpadctot = 0;
00701     for(int i = 0; i < 11; i++){
00702       if(track->smde_adc[i] > track->smde_pedestal[i])bsmdeadctot += track->smde_adc[i] - track->smde_pedestal[i];
00703       if(track->smdp_adc[i] > track->smdp_pedestal[i])bsmdpadctot += track->smdp_adc[i] - track->smdp_pedestal[i];
00704     }
00705 
00706     double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);
00707 
00708 
00709     double scaled_adc = (track->tower_adc[0] - track->tower_pedestal[0]) / track->p;
00710     double geant_scale = geant_fit->Eval(dR);
00711     scaled_adc *= geant_scale;
00712     //cout<<scaled_adc<<endl;
00713     int index = track->tower_id[0];
00714 
00715     //figure out eta and etaindex
00716     eta = helper->getEta(index);
00717     if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
00718     etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
00719 
00720     double tgain = bemctables->calib(1,track->tower_id[0])*gains[index-1]*gains2[index-1];
00721     //double tgain = gains[index-1];
00722 
00723 
00724     if((track->tower_adc[0] - track->tower_pedestal[0]) < 2.5 * track->tower_pedestal_rms[0])continue;
00725 
00726     if(track->tower_id[0] != track->tower_id_exit)continue;
00727 
00728     if(track->htTrig == 2)continue;
00729     if(track->p > 6)continue;   
00730 
00731     if(dR > 0.025)continue;
00732 
00733     if(track->dEdx < 3.4e-6)continue;
00734     if((bsmdeadctot > -1 && bsmdeadctot < 50) && (bsmdpadctot < 50 && bsmdpadctot > -1))continue;
00735 
00736     int numsis = tracks->GetEntries();
00737     if(numsis > 0)continue;
00738 
00739     float totalbtow = 0;
00740     float maxEt = 0;
00741     int maxId = -1;
00742     for(int i = 0; i < 9; i++){
00743       if(track->tower_adc[i] - track->tower_pedestal[i] < 0)continue;
00744       float theta = helper->getTheta(track->tower_id[i]);
00745       float nextEt = (track->tower_adc[i] - track->tower_pedestal[i]) * bemctables->calib(1,track->tower_id[i])*sin(theta);
00746       totalbtow += nextEt;
00747       if(nextEt > maxEt){
00748         maxEt = nextEt;
00749         maxId = i;
00750       }
00751     }
00752     if(maxId != 0) continue;
00753 
00754     eta = helper->getEta(index);
00755     float phi = helper->getPhi(index);
00756     int crate = lookup_crate(eta,phi);
00757     if(status[index-1]==1)crate_histo[crate-1]->Fill(scaled_adc*tgain);
00758   }
00759 
00760   for(int i = 0; i < ncrates; i++){
00761     sprintf(name,"crate_fit_%i",i);
00762     crate_histo[i]->Sumw2();    
00763     crate_histo[i]->Rebin(4);   
00764 
00765     crate_fit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);
00766     crate_fit[i]->SetParLimits(0,0,10.0*crate_histo[i]->GetBinContent(1));
00767     crate_fit[i]->SetParLimits(1,-10000,0);
00768     crate_fit[i]->SetParLimits(2,0,10.0*crate_histo[i]->GetMaximum());
00769     crate_fit[i]->SetParLimits(3,0,10);
00770     crate_fit[i]->SetParameter(0,crate_histo[i]->GetBinContent(1));
00771     crate_fit[i]->SetParameter(1,-crate_histo[i]->GetBinContent(1)/6.0);
00772     crate_fit[i]->SetParameter(2,-crate_histo[i]->GetMaximum());
00773     crate_fit[i]->SetParameter(3,0.929);
00774     crate_fit[i]->SetParameter(4,0.156);
00775     crate_fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
00776     crate_fit[i]->SetLineColor(kBlue);
00777     crate_fit[i]->SetLineWidth(0.6);
00778 
00779 
00780     crate_histo[i]->Fit(crate_fit[i],"rql","",0.3,1.8);
00781     float mean = crate_histo[i]->GetFunction(name)->GetParameter(3);
00782     float merr = crate_histo[i]->GetFunction(name)->GetParError(3);
00783     crateprec->SetBinContent(i+1,mean);
00784     crateprec->SetBinError(i+1,merr);
00785     cout<<"crate "<<i+1<<" "<<mean<<" "<<merr/mean<<endl;
00786   }
00787   
00788   ofstream newgain(ngname);
00789 
00790   float gains3[ntowers];
00791   float gerr3[ntowers];
00792   for(int i = 0; i < ntowers; i++){
00793     float eta = helper->getEta(i+1);
00794     float phi = helper->getPhi(i+1);
00795     int crate = lookup_crate(eta,phi);
00796     float adjust = crate_fit[crate-1]->GetParameter(3);
00797     float og = bemctables->calib(1,i)*gains[i]*gains2[i];
00798     float ng = og;
00799     float aerr = crate_fit[crate-1]->GetParError(3);
00800     float ne = sqrt(pow(gains[i]*gerr2[i],2) + pow(gains2[i]*gainerr[i],2));
00801     if(fabs(adjust-1)/aerr > 1.5){
00802       ne = sqrt(pow(ne/(adjust),2)+pow(og*aerr/(adjust*adjust),2));
00803       ng /= adjust;
00804     }
00805     newgain << i+1 << " " << ng << " " << ne << " " << status[i] << endl;
00806     gains3[i] = ng;
00807     gerr3[i] = ne;
00808   }
00809 
00810   newgain.close();
00811   */
00812   outfile.Write();
00813   outfile.Close();
00814 
00815 }
00816 
00817 //constrained double-Gaussian to fit the background
00818 double background_only_fit(double *x, double *par){
00819   double par3 = par[0]/1.5;
00820   double par4 = par[1] + 10.;
00821   double par5 = par[2] * 6.5;
00822         
00823   double fitval = 0;
00824   if(par[2] != 0){
00825     double arg1 = (x[0]-par[1])/par[2];
00826     double arg2 = (x[0]-par4)/par5;
00827     fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2);
00828   }
00829         
00830   return fitval;
00831 }
00832 
00833 double fit_function(double *x, double *par){
00834         //6-parameter fit includes
00835         //3 param electron Gaussian
00836         //par3 = relative height of peak/bg ~ 10
00837         //par4 = mean of main bg Gaussian ~ 10
00838         //par5 = width of main bg Gaussian ~ 3
00839         
00840   double par3 = 0;
00841   if(par[3] != 0)par3 = par[0] / par[3];
00842   double par6 = par3/1.5;
00843   double par7 = par[4] + 10.;
00844   double par8 = par[5] * 6.5;
00845         
00846   double fitval = 0;
00847   if(par[2] != 0 && par[5] != 0){
00848     double arg1 = (x[0]-par[1])/par[2];
00849     double arg2 = (x[0]-par[4])/par[5];
00850     double arg3 = (x[0]-par7)/par8;
00851     fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2) + par6*TMath::Exp(-0.5*arg3*arg3);
00852   }
00853         
00854   return fitval;
00855 }
00856 
00857 double fit_function2(double *x, double *par){
00858         //6-parameter fit includes
00859         //3 param electron Gaussian
00860         //par3 = relative height of peak/bg ~ 10
00861         //par4 = mean of main bg Gaussian ~ 10
00862         //par5 = width of main bg Gaussian ~ 3
00863         
00864   double par3 = 0;
00865   if(par[3] != 0)par3 = par[0] / par[3];
00866   double par6 = par[6];
00867   double par7 = par[4] + 10.;
00868   double par8 = par[7];
00869         
00870   double fitval = 0;
00871   if(par[2] != 0 && par[5] != 0 && par8 != 0){
00872     double arg1 = (x[0]-par[1])/par[2];
00873     double arg2 = (x[0]-par[4])/par[5];
00874     double arg3 = (x[0])/par8;
00875     fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2) + par6*TMath::Exp(-arg3);
00876   }
00877   return fitval;
00878 }
00879 
00880 int lookup_crate(float eta,float phi)
00881 {
00882     if (eta < 0){
00883       if (phi < -2.72) return 5;
00884       else if (phi < -2.30) return 6;
00885       else if (phi < -1.88) return 7;
00886       else if (phi < -1.46) return 8;
00887       else if( phi < -1.04) return 9;
00888       else if( phi < -0.62) return 10;
00889       else if( phi < -0.20) return 11;
00890       else if( phi <  0.22) return 12;
00891       else if( phi <  0.64) return 13;
00892       else if( phi <  1.06) return 14;
00893       else if( phi <  1.48) return 15;
00894       else if( phi <  1.90) return 1;
00895       else if( phi <  2.32) return 2;
00896       else if( phi <  2.74) return 3;
00897       else             return 4;
00898     }else{
00899       if   (phi < -2.72) return 20;
00900       else if( phi < -2.30) return 21;
00901       else if( phi < -1.88) return 22;
00902       else if( phi < -1.46) return 23;
00903       else if( phi < -1.04) return 24;
00904       else if( phi < -0.62) return 25;
00905       else if( phi < -0.20) return 26;
00906       else if( phi <  0.22) return 27;
00907       else if( phi <  0.64) return 28;
00908       else if( phi <  1.06) return 29;
00909       else if( phi <  1.48) return 30;
00910       else if( phi <  1.90) return 16;
00911       else if( phi <  2.32) return 17;
00912       else if( phi <  2.74) return 18;
00913       else             return 19;
00914     }
00915 }
StEmcPool/StEmcOfflineCalibrationMaker/macros/electron_master.C
