tpc/generator/writeTimeOrderedDistortions.C

0001 //basic framework to read and compare events in a time ordered distortion file
0002 //these files contain a TTree with xingnum, and three distortion map branches
0003
0004 #include <TCanvas.h>
0005 #include <TFile.h>
0006 #include <TFileCollection.h>
0007 #include <TFileInfo.h>
0008 #include <TH1.h>
0009 #include <TH2.h>
0010 #include <TH3.h>
0011 #include <THashList.h>
0012 #include <TTree.h>
0013
0014 #include <format>
0015 #include <iostream>
0016
0017 void writeTimeOrderedDistortions(bool subtractFirst=false, const std::string &filename="/sphenix/user/rcorliss/distortion_maps/2022.07/TimeOrderedDistortions.root", const std::string &inputpattern="/sphenix/user/rcorliss/distortion_maps/2022.07/*.distortion_map.hist.root", bool debug=false){
0018
0019   TFile *treefile=TFile::Open(filename.c_str(),"RECREATE");
0020   //TTree *tree=(TTree*)(file->Get("TimeDists"));
0021
0022   const int nhists=8;
0023   TH3F *basehist[nhists];
0024   TH3F *temphist[nhists];
0025   int xingnum=0;
0026   //note that we match the histogram RPhi to the branch P, and the histogram P to the branch Phi, to maintain backwards compatibility with older sets... for now.
0027   std::string branchname[]={"hIntDistortionP_negz",
0028                 "hIntDistortionPhi_negz",
0029                 "hIntDistortionR_negz",
0030                 "hIntDistortionZ_negz",
0031                 "hIntDistortionP_posz",
0032                 "hIntDistortionPhi_posz",
0033                 "hIntDistortionR_posz",
0034                 "hIntDistortionZ_posz"};
0035   std::string histname[]={"hIntDistortionRPhi_negz",
0036               "hIntDistortionP_negz",
0037               "hIntDistortionR_negz",
0038               "hIntDistortionZ_negz",
0039               "hIntDistortionRPhi_posz",
0040               "hIntDistortionP_posz",
0041               "hIntDistortionR_posz",
0042               "hIntDistortionZ_posz"};
0043   TTree *tree=new TTree("TimeDists", "TimeDists");
0044   tree->Branch("xingnum",&xingnum);
0045   for (int i=0;i<nhists;i++){
0046     temphist[i]=new TH3F(std::format("temphist{}",i).c_str(),std::format("temphist{}",i).c_str(),10,0,10,20,0,20,30,0,30);
0047     basehist[i]=new TH3F(std::format("basehist{}",i).c_str(),std::format("basehist{}",i).c_str(),10,0,10,20,0,20,30,0,30);
0048     tree->Branch(branchname[i].c_str(),&(temphist[i]));
0049   }
0050   if (debug) {  std::cout << "histograms built and branched." << std::endl;
0051 }
0052
0053   //find all files that match the input string
0054   TFileCollection *filelist=new TFileCollection();
0055   filelist->Add(inputpattern.c_str());
0056   filelist->Print();
0057   std::cout << "found " << filelist->GetNFiles() << " files like: " << ((TFileInfo*)(filelist->GetList()->At(0)))->GetCurrentUrl()->GetFile() << std::endl;
0058
0059   //return;
0060   TFile *infile;
0061   bool fileIsValid=true;
0062   bool isFirst=true;
0063   int nMaps=0;
0064
0065
0066   for (int i=0;i<filelist->GetNFiles() && i<100;i++){
0067     //for each file, find all histograms in that file.
0068     infile=TFile::Open(((TFileInfo*)(filelist->GetList()->At(i)))->GetCurrentUrl()->GetUrl(),"READ");//gross.
0069     fileIsValid=true;
0070     if (debug)
0071     {
0072       std::cout << "=====> Trying File " << i << std::endl;
0073     }
0074     if (!infile->IsOpen())
0075     {
0076       std::cout << "=====> File " << i << " is NOT openable <=======" << std::endl;
0077       continue; //file didn't open right.  move on to the next one.
0078     }
0079     //TList *keys=infile->GetListOfKeys();
0080     for (int j=0;j<nhists;j++){
0081       temphist[j]=nullptr;
0082       temphist[j]=infile->Get<TH3F>(histname[j].c_str());
0083       if (!temphist[j]){
0084     fileIsValid=false; //histogram doesn't exist.  don't bother loading the other hists.
0085     break;
0086       }
0087       int nbins=temphist[j]->GetNcells();
0088         if (debug)
0089     {
0090       std::cout << "=======> \"" << histname[j] << "\" has " << nbins << " cells" << std::endl;
0091     }
0092
0093     }
0094     if (!fileIsValid)
0095     {
0096       infile->Close();
0097       std::cout << "=====> File " << i << " is NOT valid <=======" << std::endl;
0098
0099       continue; //didn't get all our hists.  move on to the next file.
0100     }
0101     if (debug)
0102     {
0103       std::cout << "=====> File " << i << " is valid" << std::endl;
0104     }
0105     xingnum=i;//temporary fix to paste something in there.
0106     if(subtractFirst){
0107       if (isFirst){
0108     for (int j=0;j<nhists;j++){
0109       treefile->cd();
0110       basehist[j]=(TH3F*)(temphist[j]->Clone());
0111       //temphist[j]->Copy(*(basehist[j]));
0112     }
0113     isFirst=false;
0114       }
0115       for (int j=0;j<nhists;j++){
0116     std::cout << "ptr j=" << j << ":  b:" << basehist[j] << "\tt:" << temphist[j] << std::endl;
0117     int nbins=temphist[j]->GetNcells();
0118     for (int k=0;k<nbins;k++){
0119       double b=basehist[j]->GetBinContent(k);
0120       double t=temphist[j]->GetBinContent(k);
0121       double diff=t-b;
0122       temphist[j]->SetBinContent(k,diff);
0123       //temphist[j]->Add(basehist[j],-1);
0124     }
0125     if (debug)
0126     {
0127       std::cout << "=======> \"" << histname[j] << "\" has " << nbins << " cells when writing diff" << std::endl;
0128     }
0129
0130       }
0131     }
0132
0133     tree->Fill();
0134     nMaps++;
0135     infile->Close();
0136   }
0137   std::cout << "finished tree " << filename << " has nMaps=" << nMaps << std::endl;
0138
0139   treefile->cd();
0140   tree->Write();
0141   treefile->Close();
0142   return;
0143 }