tpc/cm/CMDistortionRecoCart.C

0001 // step 2
0002 #include <TCanvas.h>
0003 #include <TFile.h>
0004 #include <TFileCollection.h>
0005 #include <TFileInfo.h>
0006 #include <TH1.h>
0007 #include <TH2.h>
0008 #include <TH3.h>
0009 #include <THashList.h>
0010 #include <TMath.h>
0011 #include <TSystem.h>
0012 #include <TTree.h>
0013 #include <TTime.h>
0014 #include <TVector3.h>
0015
0016 #include <cmath>
0017 #include <format>
0018 #include <iostream>
0019 #include <vector>
0020
0021 int CMDistortionRecoCart(int nMaxEvents = -1) {
0022   int nbins = 35;
0023   double low = -80.0;
0024   double high = 80.0;
0025 //  double deltaX, deltaY, deltaZ, deltaR, deltaPhi;
0026   int nEvents;
0027
0028   //take in events
0029   const char * inputpattern="/sphenix/u/skurdi/CMCalibration/cmDistHitsTree_Event*.root";
0030
0031   //find all files that match the input string (includes wildcards)
0032   TFileCollection *filelist=new TFileCollection();
0033   filelist->Add(inputpattern);
0034   TString sourcefilename;
0035
0036   //how many events
0037   if(nMaxEvents >= 0 && nMaxEvents<filelist->GetNFiles())
0038   {
0039     nEvents=nMaxEvents;
0040   }
0041   else
0042   {
0043     nEvents= filelist->GetNFiles();
0044   }
0045
0046   //canvas for checking data
0047   //TCanvas *canvas1=new TCanvas("canvas1","CMDistortionRecoCart1",1200,800);
0048   //canvas1->Divide(3,2);
0049
0050   //canvas for time plot
0051   TCanvas *canvas=new TCanvas("canvas","CMDistortionRecoCart2",400,400);
0052
0053   TVector3 *position;
0054   TVector3 *newposition;
0055   position = new TVector3(1.,1.,1.);
0056   newposition = new TVector3(1.,1.,1.);
0057
0058   //histogram to compare times
0059   TH1F *hTimePerEvent = new TH1F("hTimePerEvent","Time Per Event; time (ms)",20,0,10000);
0060
0061   for (int ifile=0;ifile < nEvents;ifile++){
0062     //call to TTime before opening ttree
0063     TTime now;
0064     now=gSystem->Now();
0065     unsigned long before = now;
0066
0067     //get data from ttree
0068     sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
0069
0070 //    char const *treename="cmDistHitsTree";
0071     TFile *input=TFile::Open(sourcefilename, "READ");
0072     TTree *inTree=(TTree*)input->Get("tree");
0073
0074     inTree->SetBranchAddress("position",&position);
0075     inTree->SetBranchAddress("newposition",&newposition);
0076
0077     //for forward only
0078
0079     TH2F *hStripesPerBin = new TH2F("hStripesPerBin","CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0080
0081     TH2F *hCartesianForward[3];
0082     hCartesianForward[0] = new TH2F("hForwardX","X Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0083     hCartesianForward[1] = new TH2F("hForwardY","Y Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0084     hCartesianForward[2] = new TH2F("hForwardZ","Z Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0085
0086     for (int i=0;i<inTree->GetEntries();i++){
0087       inTree->GetEntry(i);
0088
0089       hStripesPerBin->Fill(position->X(),position->Y(),1);
0090
0091       double deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
0092       double deltaY = (newposition->Y() - position->Y())*(1e4);
0093       double deltaZ = (newposition->Z() - position->Z())*(1e4);
0094
0095       // double deltaR = (newposition->Perp() - position->Perp())*(1e4);
0096       // double deltaPhi = newposition->DeltaPhi(*position);
0097
0098       hCartesianForward[0]->Fill(position->X(),position->Y(),deltaX);
0099       hCartesianForward[1]->Fill(position->X(),position->Y(),deltaY);
0100       hCartesianForward[2]->Fill(position->X(),position->Y(),deltaZ);
0101     }
0102
0103     TH2F *hCartesianAveShift[3];
0104     hCartesianAveShift[0] = new TH2F("AveShiftX","Average of CM Model X over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0105     hCartesianAveShift[1] = new TH2F("AveShiftY","Average of CM Model Y over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0106     hCartesianAveShift[2] = new TH2F("AveShiftZ","Average of CM Model Z over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0107
0108     TH2F *hCylindricalAveShift[2];
0109      hCylindricalAveShift[0] = new TH2F("AveShiftRCart","Average of CM Model R over Stripes per Bin from Cartesian (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0110     hCylindricalAveShift[1] = new TH2F("AveShiftPhiCart","Average of CM Model Phi over Stripes per Bin from Cartesian (rad); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0111
0112     for (int i = 0; i < 3; i ++){
0113       hCartesianAveShift[i]->Divide(hCartesianForward[i],hStripesPerBin);
0114     }
0115
0116     //cyl models from cart coords
0117     for(int i = 0; i < nbins; i++){
0118       double x = low + ((high - low)/(1.0*nbins))*(i+0.5); //center of bin
0119
0120       for(int j = 0; j < nbins; j++){
0121     double y = low + ((high - low)/(1.0*nbins))*(j+0.5); //center of bin
0122
0123     int xbin = hCartesianAveShift[0]->FindBin(x,y);
0124     int ybin = hCartesianAveShift[1]->FindBin(x,y);
0125     double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1e-4); // converts  microns to cm
0126     double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1e-4);
0127
0128     TVector3 shifted;
0129     TVector3 original;
0130     original.SetX(x);
0131     original.SetY(y);
0132     shifted.SetX(x+xaveshift);
0133     shifted.SetY(y+yaveshift);
0134     //x n y above for orig
0135     //shifted is orig + ave shift
0136
0137     double raveshift = (shifted.Perp() - original.Perp())*(1e4);
0138     double phiaveshift = shifted.DeltaPhi(original);
0139
0140     //fill with r from x n y
0141     hCylindricalAveShift[0]->Fill(x,y,raveshift);
0142     hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
0143       }
0144     }
0145
0146     //same range and bins for each coordinate, binned in cm
0147     //hardcoded numbers from average distortion file's hIntDistortionPosX
0148     int nphi = 82;
0149     int nr = 54;
0150     int nz = 82;
0151
0152     double minphi = -0.078539819;
0153     double minr = 18.884615;
0154     //double minz = 5.0;
0155     double minz = -1.3187500;
0156
0157     double maxphi = 6.3617253;
0158     double maxr = 79.115387;
0159     double maxz = 106.81875;
0160
0161     TH3F *hCartesianCMModel[3];
0162     hCartesianCMModel[0]=new TH3F("hCMModelX", "CM Model: X Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz); //rad, cm, cm
0163     hCartesianCMModel[1]=new TH3F("hCMModelY", "CM Model: Y Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0164     hCartesianCMModel[2]=new TH3F("hCMModelZ", "CM Model: Z Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0165
0166     TH3F *hCylindricalCMModel[2];
0167     hCylindricalCMModel[0]=new TH3F("hCMModelRCart", "CM Model: Radial Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0168     hCylindricalCMModel[1]=new TH3F("hCMModelPhiCart", "CM Model: Phi Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0169
0170     double xshift;
0171     double yshift;
0172     double zshift;
0173     double rshiftcart;
0174     double phishiftcart;
0175
0176     for(int i = 0; i < nphi; i++){
0177       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0178
0179       for(int j = 0; j < nr; j++){
0180     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0181
0182     double x = r*cos(phi); //cm
0183     double y = r*sin(phi);
0184
0185     for(int k = 0; k < nz; k++){
0186       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0187
0188       xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1e-4);//coordinate of your stripe
0189       yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1e-4);//convert micron to cm
0190       zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1e-4);
0191
0192       rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1e-4);
0193       phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
0194
0195       hCartesianCMModel[0]->Fill(phi,r,z,xshift*(1-z/105.5));
0196       hCartesianCMModel[1]->Fill(phi,r,z,yshift*(1-z/105.5));
0197       hCartesianCMModel[2]->Fill(phi,r,z,zshift*(1-z/105.5));
0198
0199       hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
0200       hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
0201     }
0202       }
0203     }
0204
0205     TFile *plots;
0206
0207     plots=TFile::Open(std::format("CMModelsCart_Event{}.root",ifile).c_str(),"RECREATE");
0208     hStripesPerBin->Write();
0209
0210     for(int i = 0; i < 3; i++){
0211       hCartesianForward[i]->Write();
0212       hCartesianAveShift[i]->Write();
0213       hCartesianCMModel[i]->Write();
0214     }
0215
0216     for(int i = 0; i < 2; i++){
0217       hCylindricalAveShift[i]->Write();
0218       hCylindricalCMModel[i]->Write();
0219     }
0220
0221     plots->Close();
0222
0223
0224     //call to TTime after outputting TH3Fs
0225     now=gSystem->Now();
0226     unsigned long after = now;
0227
0228     hTimePerEvent->Fill(after-before);
0229
0230     /*
0231     //to check histograms
0232     for (int i = 0; i < 3; i++){
0233       hCartesianForward[i]->SetStats(0);
0234     }
0235
0236     hStripesPerBin->SetStats(0);
0237     canvas1->cd(1);
0238     hCartesianForward[0]->Draw("colz");
0239     canvas1->cd(2);
0240     hCartesianForward[1]->Draw("colz");
0241     canvas1->cd(3);
0242     hCartesianForward[2]->Draw("colz");
0243     canvas1->cd(4)->Clear();
0244     canvas1->cd(5)->Clear();
0245     canvas1->cd(6)->Clear();
0246
0247     if(ifile == 0){
0248       canvas1->Print("CMDistortionRecoCart1.pdf(","pdf");
0249     } else if (ifile == nEvents - 1){
0250       canvas1->Print("CMDistortionRecoCart1.pdf)","pdf");
0251     } else {
0252       canvas1->Print("CMDistortionRecoCart1.pdf","pdf");
0253     }*/
0254
0255
0256   }
0257
0258   canvas->cd();
0259   hTimePerEvent->Draw();
0260   canvas->Print("CMDistortionRecoCart2.pdf","pdf");
0261
0262   return 0;
0263 }