tpc/cm/CMDistortionRecoCart.C

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