tpc/cm/CMDistortionRecoPhiR.C

0001 // step 2 with phi,r coords
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 CMDistortionRecoPhiR(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   TCanvas *canvas1=new TCanvas("canvas1","CMDistortionReco1",1200,800);
0037   canvas1->Divide(3,2);
0038
0039   //canvas for time plot
0040   TCanvas *canvas=new TCanvas("canvas","CMDistortionReco2",400,400);
0041
0042   TVector3 *position, *newposition;
0043   position = new TVector3(1.,1.,1.);
0044   newposition = new TVector3(1.,1.,1.);
0045
0046   //histogram to compare times
0047   TH1F *hTimePerEvent = new TH1F("hTimePerEvent","Time Per Event; time (ms)",20,0,10000);
0048
0049   for (int ifile=0;ifile < nEvents;ifile++){
0050     //call to TTime before opening ttree
0051     TTime now;
0052     now=gSystem->Now();
0053     unsigned long before = now;
0054
0055     //get data from ttree
0056     sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
0057
0058     char const *treename="cmDistHitsTree";
0059     TFile *input=TFile::Open(sourcefilename, "READ");
0060     TTree *inTree=(TTree*)input->Get("tree");
0061
0062     inTree->SetBranchAddress("position",&position);
0063     inTree->SetBranchAddress("newposition",&newposition);
0064
0065     //hardcoded numbers from average distortion file's hIntDistortionPosX
0066     int nbinsphi = 30; //when using 35, blank spots at around r = 22 cm, phi just above n below pi
0067     double lowphi = -0.078539819;
0068     double highphi = 6.3617253;
0069     int nbinsr = 30; // when using 35, blank stripe around r = 58 cm
0070     double lowr = 0.0;
0071     double highr = 90.0;
0072
0073     //for forward only
0074
0075     TH2F *hStripesPerBinPhiR = new TH2F("hStripesPerBinPhiR","CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0076
0077     TH2F *hCartesianForwardPhiR[3];
0078     hCartesianForwardPhiR[0] = new TH2F("hForwardX_PhiR","X Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0079     hCartesianForwardPhiR[1] = new TH2F("hForwardY_PhiR","Y Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0080     hCartesianForwardPhiR[2] = new TH2F("hForwardZ_PhiR","Z Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0081
0082      TH2F *hCylindricalForwardPhiR[2];
0083     hCylindricalForwardPhiR[0] = new TH2F("hForwardR_PhiR","Radial Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0084     hCylindricalForwardPhiR[1] = new TH2F("hForwardPhi_PhiR","Phi Shift Forward of Stripe Centers, Phi,R binning (rad); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0085
0086     for (int i=0;i<inTree->GetEntries();i++){
0087       inTree->GetEntry(i);
0088
0089       double r = position->Perp();
0090       double phi = position->Phi();
0091
0092       if(position->Phi() < 0.0){
0093     phi = position->Phi() + 2.0*TMath::Pi();
0094       }
0095
0096       hStripesPerBinPhiR->Fill(phi,r,1);
0097
0098       deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
0099       deltaY = (newposition->Y() - position->Y())*(1e4);
0100       deltaZ = (newposition->Z() - position->Z())*(1e4);
0101
0102       deltaR = (newposition->Perp() - position->Perp())*(1e4);
0103       deltaPhi = newposition->DeltaPhi(*position);
0104
0105       hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
0106       hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
0107       hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
0108
0109       hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
0110       hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
0111     }
0112
0113     TH2F *hCartesianAveShiftPhiR[3];
0114     hCartesianAveShiftPhiR[0] = new TH2F("AveShiftX_PhiR","Average of CM Model X over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0115     hCartesianAveShiftPhiR[1] = new TH2F("AveShiftY_PhiR","Average of CM Model Y over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0116     hCartesianAveShiftPhiR[2] = new TH2F("AveShiftZ_PhiR","Average of CM Model Z over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0117
0118     TH2F *hCylindricalAveShiftPhiR[2];
0119      hCylindricalAveShiftPhiR[0] = new TH2F("AveShiftR_PhiR","Average of CM Model R over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0120     hCylindricalAveShiftPhiR[1] = new TH2F("AveShiftPhi_PhiR","Average of CM Model Phi over Stripes per Bin, Phi,R binning (rad); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0121
0122     for (int i = 0; i < 3; i ++){
0123       hCartesianAveShiftPhiR[i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
0124     }
0125
0126     hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
0127     hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
0128
0129     //same range and bins for each coordinate, binned in cm
0130     //hardcoded numbers from average distortion file's hIntDistortionPosX
0131     int nphi = 82;
0132     int nr = 54;
0133     int nz = 82;
0134
0135     double minphi = -0.078539819;
0136     double minr = 18.884615;
0137     double minz = -1.3187500;
0138
0139     double maxphi = 6.3617253;
0140     double maxr = 79.115387;
0141     double maxz = 106.81875;
0142
0143     TH3F *hCartesianCMModelPhiR[3];
0144     hCartesianCMModelPhiR[0]=new TH3F("hCMModelX_PhiR", "CM Model: X Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz); //rad, cm, cm
0145     hCartesianCMModelPhiR[1]=new TH3F("hCMModelY_PhiR", "CM Model: Y Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0146     hCartesianCMModelPhiR[2]=new TH3F("hCMModelZ_PhiR", "CM Model: Z Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0147
0148     TH3F *hCylindricalCMModelPhiR[2];
0149     hCylindricalCMModelPhiR[0]=new TH3F("hCMModelR_PhiR", "CM Model: Radial Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0150     hCylindricalCMModelPhiR[1]=new TH3F("hCMModelPhi_PhiR", "CM Model: Phi Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0151
0152     double xshiftPhiR, yshiftPhiR, zshiftPhiR, rshiftcartPhiR, phishiftcartPhiR;
0153
0154     for(int i = 0; i < nphi; i++){
0155       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0156
0157       for(int j = 0; j < nr; j++){
0158     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0159
0160     for(int k = 0; k < nz; k++){
0161       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0162
0163       xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);//coordinate of your stripe
0164       yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1e-4);//convert micron to cm
0165       zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1e-4);
0166
0167       rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);
0168       phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
0169
0170       hCartesianCMModelPhiR[0]->Fill(phi,r,z,xshiftPhiR*(1-z/105.5));
0171       hCartesianCMModelPhiR[1]->Fill(phi,r,z,yshiftPhiR*(1-z/105.5));
0172       hCartesianCMModelPhiR[2]->Fill(phi,r,z,zshiftPhiR*(1-z/105.5));
0173
0174       hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
0175       hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
0176     }
0177       }
0178     }
0179
0180     TFile *plots;
0181
0182     plots=TFile::Open(Form("CMModelsPhiR_Event%d.root",ifile),"RECREATE");
0183
0184     for(int i = 0; i < 3; i++){
0185       hCartesianForwardPhiR[i]->Write();
0186       hCartesianAveShiftPhiR[i]->Write();
0187       hCartesianCMModelPhiR[i]->Write();
0188     }
0189
0190     for(int i = 0; i < 2; i++){
0191       hCylindricalForwardPhiR[i]->Write();
0192       hCylindricalAveShiftPhiR[i]->Write();
0193       hCylindricalCMModelPhiR[i]->Write();
0194     }
0195
0196     plots->Close();
0197
0198
0199     //call to TTime after outputting TH3Fs
0200     now=gSystem->Now();
0201     unsigned long after = now;
0202
0203     hTimePerEvent->Fill(after-before);
0204
0205     //to check histograms
0206     for (int i = 0; i < 3; i++){
0207       hCartesianForwardPhiR[i]->SetStats(0);
0208     }
0209
0210     canvas1->cd(1);
0211     hStripesPerBinPhiR->Draw("colz");
0212     canvas1->cd(2)->Clear();
0213     canvas1->cd(3)->Clear();
0214      canvas1->cd(4);
0215     hCartesianForwardPhiR[0]->Draw("colz");
0216     canvas1->cd(5);
0217     hCartesianForwardPhiR[1]->Draw("colz");
0218     canvas1->cd(6);
0219     hCartesianForwardPhiR[2]->Draw("colz");
0220
0221     if(ifile == 0){
0222       canvas1->Print("CMDistortionReco1.pdf(","pdf");
0223     } else if (ifile == nEvents - 1){
0224       canvas1->Print("CMDistortionReco1.pdf)","pdf");
0225     } else {
0226       canvas1->Print("CMDistortionReco1.pdf","pdf");
0227     }
0228
0229   }
0230
0231   canvas->cd();
0232   hTimePerEvent->Draw();
0233   canvas->Print("CMDistortionReco2.pdf","pdf");
0234
0235   return 0;
0236 }