test/Validation/HcalCompare.C

0001 #include "HcalCompare.h"
0002
0003 void compare(char* hardware_file, char* software_file)
0004 {
0005   TStyle *mystyle = new TStyle("mystyle","my style");
0006   initStyle(mystyle);
0007   mystyle->cd();
0008   index_map::index_map get_index;
0009   TCanvas *c1 = new TCanvas("c1","comparison");
0010   TFile *file_hard = TFile::Open(hardware_file);
0011   if(!file_hard)
0012     {
0013       cout << "No hardware file found!";
0014       return;
0015     }
0016   TFile *file_soft = TFile::Open(software_file);
0017   if(!file_soft)
0018     {
0019       cout << "No software file found!";
0020       return;
0021     }
0022
0023   TTree* tree_hard = (TTree*)file_hard->Get("TPGntuple");
0024   if(!tree_hard)
0025     {
0026       cout << "No tree found for hardware!";
0027       return;
0028     }
0029
0030   TTree* tree_soft = (TTree*)file_soft->Get("TPGntuple");
0031   if(!tree_soft)
0032     {
0033       cout << "No tree found for software!";
0034       return;
0035     }
0036
0037   float tpg_energy_hard[4176], tpg_energy_soft[4176];
0038   int ieta_hard[4176], iphi_hard[4176], ieta_soft[4176], iphi_soft[4176], index_hard[4176], index_soft[4176];
0039
0040   tree_hard->SetBranchAddress("tpg_energy",tpg_energy_hard);
0041   tree_hard->SetBranchAddress("ieta",ieta_hard);
0042   tree_hard->SetBranchAddress("iphi",iphi_hard);
0043   tree_hard->SetBranchAddress("tpg_index",index_hard);
0044
0045   tree_soft->SetBranchAddress("tpg_energy",tpg_energy_soft);
0046   tree_soft->SetBranchAddress("ieta",ieta_soft);
0047   tree_soft->SetBranchAddress("iphi",iphi_soft);
0048   tree_soft->SetBranchAddress("tpg_index",index_soft);
0049
0050   TH1F *comparison = new TH1F("comparison","tpg_hardware - tpg_software",400,-200,200);
0051   TH2F *mean_2d = new TH2F("res2D_mean","Mean of TP_{hardware}-TP_{software}",65,-32,33,72,1,73);
0052   TH2F *rms_2d = new TH2F("res2D_rms","RMS of TP_{hardware}-TP_{software}",65,-32,33,72,1,73);
0053
0054
0055   TH1F *hardware = new TH1F("hardware","hardware",200,0,200);
0056   TH1F *software = new TH1F("software","software",200,0,200);
0057
0058   TObjArray hard_minus_soft(4176);
0059   TObjArray hard_v_soft(4176);
0060
0061   char name[20], title[20];
0062
0063
0064   int true_index[4176];
0065   int limit;
0066   //  cout << "^^^^^\nNumber of Entries: hardware = " << tree_hard->GetEntries() << " software = " << tree_soft->GetEntries() << "\n";
0067   if ((int)tree_hard->GetEntries() <= (int)tree_soft->GetEntries())
0068     {
0069       limit = (int)tree_hard->GetEntries();
0070     }
0071   else
0072     {
0073       limit = (int)tree_soft->GetEntries();
0074     }
0075
0076   map<int,float> tpg_hard_map, tpg_soft_map;
0077   for(int i=0; i < limit; i++)
0078     {
0079       tree_hard->GetEntry(i);
0080       tree_soft->GetEntry(i);
0081       //cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$\n";
0082       tpg_hard_map.clear();
0083       tpg_soft_map.clear();
0084       for(int j=0; j < 4176; j++)
0085     {
0086       true_index[j] = get_index.ntpg(j);
0087       if (index_hard[j] != 0)
0088         {
0089           tpg_hard_map[index_hard[j]] = tpg_energy_hard[j];
0090         }
0091       if (index_soft[j] != 0)
0092         {
0093           tpg_soft_map[index_soft[j]] = tpg_energy_soft[j];
0094         }
0095       //initialize
0096       if (i==0)
0097         {
0098           sprintf(name,"difference_%d",true_index[j]);
0099           sprintf(title,"difference:%d",true_index[j]);
0100           hard_minus_soft[j] = new TH1F(name,title,160,-80,80);
0101           sprintf(name,"hard_v_soft_%d",true_index[j]);
0102           sprintf(title,"hard_v_soft:%d",true_index[j]);
0103           hard_v_soft[j] = new TH2F(name,title,200,0,200,200,0,200);
0104         }
0105     }
0106
0107       for(int j=0; j < 4176; j++)
0108     {
0109       int tpg_hard, tpg_soft;
0110       if (tpg_hard_map.count(true_index[j])==0)
0111         {
0112           tpg_hard = -1;
0113         }
0114       else
0115         {
0116           tpg_hard = (int)tpg_hard_map.find(true_index[j])->second;
0117         }
0118           if (tpg_soft_map.count(true_index[j])==0)
0119             {
0120               tpg_soft = -1;
0121             }
0122           else
0123             {
0124               tpg_soft = (int)tpg_soft_map.find(true_index[j])->second;
0125             }
0126       //fill histos
0127       if(tpg_hard != -1 && tpg_soft != -1)
0128         {
0129           ((TH2F*)hard_v_soft[j])->Fill(tpg_soft,tpg_hard);
0130           hardware->Fill(tpg_hard);
0131           software->Fill(tpg_soft);
0132           comparison->Fill(tpg_hard-tpg_soft);
0133           ((TH1F*)hard_minus_soft[j])->Fill(tpg_hard-tpg_soft);
0134         }
0135     }
0136
0137       int hard_count = 0;
0138       int soft_count = 0;
0139       //cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n";
0140       for(int j=0; j < 4176; j++)
0141     {
0142       if (tpg_energy_hard[j] != 0 || (tpg_energy_hard[j] == 0 && index_hard[j] !=0))
0143         {
0144           hard_count++;
0145           cout << "Found hardware channel: " << index_hard[j] << " tpg counts = " << tpg_energy_hard[j] << "\n";
0146         }
0147       if (tpg_energy_soft[j] !=0  || (tpg_energy_soft[j] == 0 && index_soft[j] !=0))
0148         {
0149           soft_count++;
0150           cout << "Found software channel: " << index_soft[j] << " tpg counts = " << tpg_energy_soft[j] << "\n";
0151         }
0152     }
0153       cout << "N hardware channels = " << hard_count << " N software channels = " << soft_count << "\n";
0154
0155     }
0156
0157   double mean;
0158   double sigma;
0159   TF1 *line = new TF1("line","x",0,200);
0160   line->SetLineColor(4);
0161   int eta;
0162   int phi;
0163   for (int j=0; j<4176; ++j)
0164     {
0165       mean = ((TH1F*)hard_minus_soft[j])->GetMean(1);
0166       sigma = ((TH1F*)hard_minus_soft[j])->GetRMS(1);
0167       eta = true_index[j]/100;
0168       phi = TMath::Abs(true_index[j]%100);
0169       mean_2d->Fill(eta,phi,mean);
0170       rms_2d->Fill(eta,phi,sigma);
0171       ((TH2F*)hard_v_soft[j])->Draw();
0172       line->Draw("same");
0173     }
0174
0175   c1->SetLogy();
0176   comparison->Draw();
0177   SetupTitle(comparison,"Hardware TP - Software TP (counts)", "Ntowers");
0178   if (comparison->GetMean(1) == 0 && comparison->GetRMS(1) == 0)
0179     {
0180       SetStatus(comparison,"GOOD");
0181     }
0182   else
0183     {
0184       SetStatus(comparison,"BAD");
0185     }
0186   c1->Print("compare_output.ps(");
0187
0188   TCanvas *c2 = new TCanvas("c2","Mean");
0189   c2->cd();
0190   //  mean_2d->SetStats(kFALSE);
0191   mean_2d->Draw("COLZ");
0192   SetupTowerDisplay(mean_2d);
0193   c2->Print("compare_output.ps");
0194
0195   TCanvas *c3 = new TCanvas("c3","RMS");
0196   c3->cd();
0197   //rms_2d->SetStats(kFALSE);
0198   rms_2d->Draw("COLZ");
0199   SetupTowerDisplay(rms_2d);
0200   c3->Print("compare_output.ps");
0201
0202   TCanvas *c4 = new TCanvas("c4","hard and soft");
0203   c4->Divide(1,2);
0204   c4->cd(1);
0205   gPad->SetLogy();
0206   hardware->Draw();
0207   c4->cd(2);
0208   gPad->SetLogy();
0209   software->Draw();
0210   c4->Print("compare_output.ps)");
0211
0212   TFile *f = new TFile("output.root","recreate");
0213   hard_v_soft.Write();
0214   hard_minus_soft.Write();
0215   f->Close();
0216
0217   return;
0218 }