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 #include "TTree.h" //this prevents a lazy binding issue and/or is a magic spell.
0004 #include "TCanvas.h" //this prevents a lazy binding issue and/or is a magic spell.
0005 #include "TH3F.h"
0006 #include "TH2F.h"
0007 #include "TH1F.h"
0008
0009 void writeTimeOrderedDistortions(bool subtractFirst=false, char *filename="/sphenix/user/rcorliss/distortion_maps/2022.07/TimeOrderedDistortions.root", char *inputpattern="/sphenix/user/rcorliss/distortion_maps/2022.07/*.distortion_map.hist.root", bool debug=false){
0010
0011   TFile *treefile=TFile::Open(filename,"RECREATE");
0012   //TTree *tree=(TTree*)(file->Get("TimeDists"));
0013
0014   int nhists=8;
0015   TH3F *basehist[nhists];
0016   TH3F *temphist[nhists];
0017   int xingnum=0;
0018   //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.
0019   std::string branchname[]={"hIntDistortionP_negz",
0020                 "hIntDistortionPhi_negz",
0021                 "hIntDistortionR_negz",
0022                 "hIntDistortionZ_negz",
0023                 "hIntDistortionP_posz",
0024                 "hIntDistortionPhi_posz",
0025                 "hIntDistortionR_posz",
0026                 "hIntDistortionZ_posz"};
0027   std::string histname[]={"hIntDistortionRPhi_negz",
0028               "hIntDistortionP_negz",
0029               "hIntDistortionR_negz",
0030               "hIntDistortionZ_negz",
0031               "hIntDistortionRPhi_posz",
0032               "hIntDistortionP_posz",
0033               "hIntDistortionR_posz",
0034               "hIntDistortionZ_posz"};
0035   TTree *tree=new TTree("TimeDists", "TimeDists");
0036   tree->Branch("xingnum",&xingnum);
0037   for (int i=0;i<nhists;i++){
0038     temphist[i]=new TH3F(Form("temphist%d",i),Form("temphist%d",i),10,0,10,20,0,20,30,0,30);
0039     basehist[i]=new TH3F(Form("basehist%d",i),Form("basehist%d",i),10,0,10,20,0,20,30,0,30);
0040     tree->Branch(branchname[i].c_str(),&(temphist[i]));
0041   }
0042   if (debug)  printf("histograms built and branched.\n");
0043
0044   //find all files that match the input string
0045   TFileCollection *filelist=new TFileCollection();
0046   filelist->Add(inputpattern);
0047   filelist->Print();
0048   printf("found %d files like: %s\n",filelist->GetNFiles(),((TFileInfo*)(filelist->GetList()->At(0)))->GetCurrentUrl()->GetFile());//Title());//Print();
0049
0050   //return;
0051   TFile *infile;
0052   bool fileIsValid=true;
0053   bool isFirst=true;
0054   int nMaps=0;
0055
0056
0057   for (int i=0;i<filelist->GetNFiles() && i<100;i++){
0058     //for each file, find all histograms in that file.
0059     infile=TFile::Open(((TFileInfo*)(filelist->GetList()->At(i)))->GetCurrentUrl()->GetUrl(),"READ");//gross.
0060     fileIsValid=true;
0061     if (debug)   printf("=====> Trying File %d\n",i);
0062     if (!infile->IsOpen()) {
0063           printf("=====> File %d is NOT openable <=======\n",i);
0064       continue; //file didn't open right.  move on to the next one.
0065     }
0066     //TList *keys=infile->GetListOfKeys();
0067     for (int j=0;j<nhists;j++){
0068       temphist[j]=NULL;
0069       temphist[j]=infile->Get<TH3F>(histname[j].c_str());
0070       if (!temphist[j]){
0071     fileIsValid=false; //histogram doesn't exist.  don't bother loading the other hists.
0072     break;
0073       }
0074       int nbins=temphist[j]->GetNcells();
0075         if (debug) printf("=======> \"%s\" has %d cells\n",histname[j].c_str(),nbins);
0076
0077     }
0078     if (!fileIsValid) {
0079       infile->Close();
0080           printf("=====> File %d is NOT valid <=======\n",i);
0081
0082     continue; //didn't get all our hists.  move on to the next file.
0083     }
0084       if (debug) printf("=====> File %d is valid\n",i);
0085     xingnum=i;//temporary fix to paste something in there.
0086     if(subtractFirst){
0087       if (isFirst){
0088     for (int j=0;j<nhists;j++){
0089       treefile->cd();
0090       basehist[j]=(TH3F*)(temphist[j]->Clone());
0091       //temphist[j]->Copy(*(basehist[j]));
0092     }
0093     isFirst=false;
0094       }
0095       for (int j=0;j<nhists;j++){
0096     printf("ptr j=%d:  b:%p\tt:%p\n",j,basehist[j],temphist[j]);
0097     int nbins=temphist[j]->GetNcells();
0098     for (int k=0;k<nbins;k++){
0099       double b=basehist[j]->GetBinContent(k);
0100       double t=temphist[j]->GetBinContent(k);
0101       double diff=t-b;
0102       temphist[j]->SetBinContent(k,diff);
0103       //temphist[j]->Add(basehist[j],-1);
0104     }
0105       if (debug) printf("=======> \"%s\" has %d cells when writing diff\n",histname[j].c_str(),nbins);
0106
0107       }
0108     }
0109
0110     tree->Fill();
0111     nMaps++;
0112     infile->Close();
0113   }
0114   printf("finished tree %s has nMaps=%d\n",filename, nMaps);
0115
0116   treefile->cd();
0117   tree->Write();
0118   treefile->Close();
0119   return;
0120 }