00001 #include "subtract.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 const char* infileA="links/DEV.nofield.1244036.noExB.refitOS.slice/dip5.slice.root";
00016 const char* infileB="links/P00hk.nofield.refitOS.0917.slice/dip5.typec.slice.new.root";
00017 const char* outfile="links/DEV.nofield.1244036.noExB.refitOS.slice/dip5.typec.DEVminusP00hk.slice.root";
00018
00019
00020 int main(int argc,char** argv)
00021 {
00022 char* argvZ[] = {"-b"};
00023 Int_t argcZ = 1;
00024 TApplication r00t("r00t", &argcZ, argvZ);
00025
00026 extern char* optarg;
00027 const char* options = "a:b:o:";
00028 Int_t chr, a=0,b=0,o=0;
00029 char inFileA[300],inFileB[300],outFile[300];
00030 while ((chr = getopt(argc, argv, options)) >= 0){
00031 switch(chr){
00032 case 'a': strcpy(inFileA,optarg); a=1; break;
00033 case 'b': strcpy(inFileB,optarg); b=1; break;
00034 case 'o': strcpy(outFile,optarg); o=1; break;
00035 }
00036 }
00037 if(!(a*b*o)){
00038 cout << "-a infileA -b infileB -o outFile" << endl;
00039 exit(-1);
00040 }
00041
00042
00043
00044
00045
00046 cout << "infileA=" << inFileA << endl
00047 << "infileB=" << inFileB << endl
00048 << "outFile=" << outFile << endl;
00049
00050 TFile inRootA(inFileA); if(!inRootA.IsOpen()) return -1;
00051 TFile inRootB(inFileB); if(!inRootB.IsOpen()) return -1;
00052 TFile outRoot(outFile,"RECREATE"); if(!outRoot.IsOpen()) return -1;
00053
00054 TIterator* iterator = inRootA.GetListOfKeys()->MakeIterator();
00055 TKey* key;
00056
00057 while( (key=dynamic_cast<TKey*>(iterator->Next())) != 0){
00058
00059 TH1 *hA=0,*hB=0,*hOut=0;
00060 hA = (TH1*)inRootA.Get(key->GetName());
00061 if(hA->GetDimension()!=1) continue;
00062
00063 hB = (TH1*)inRootB.Get(key->GetName());
00064 if(!hB) {
00065 cout << "Cannot find " << key->GetName()
00066 << " in file " << inFileB << endl;
00067 continue;
00068 }
00069 if(hA->GetNbinsX()!=hB->GetNbinsX() ||
00070 hA->GetXaxis()->GetBinLowEdge(1) != hB->GetXaxis()->GetBinLowEdge(1)){
00071 cout << "Different bins " << key->GetName() << endl;
00072 continue;
00073 }
00074 if(hA->GetXaxis()->GetXbins()){
00075 hOut = new TH1D(hA->GetName(),hA->GetTitle(),
00076 hA->GetNbinsX(),hA->GetXaxis()->GetXbins()->GetArray());
00077 }
00078 else{
00079 hOut = new TH1D(hA->GetName(),hA->GetTitle(),
00080 hA->GetNbinsX(),
00081 hA->GetXaxis()->GetBinLowEdge(1),
00082 hA->GetXaxis()->GetBinUpEdge(hA->GetNbinsX()));
00083 }
00084
00085 for(int i=1; i<=hA->GetNbinsX(); i++){
00086 double val=hA->GetBinContent(i)-hB->GetBinContent(i);
00087 double err=sqrt(hA->GetBinError(i)*hA->GetBinError(i)+
00088 hB->GetBinError(i)*hB->GetBinError(i));
00089 hOut->SetBinContent(i,val);
00090 hOut->SetBinError(i,err);
00091 }
00092 hOut->Write();
00093
00094
00095 }
00096 outRoot.Close();
00097
00098 cout << "done" << endl;
00099
00100
00101 exit(0);
00102 }