tpc/cm/CMDistortionReco.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 CMDistortionReco(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   //  nEvents = 2;
0036
0037   TCanvas *canvas1=new TCanvas("canvas1","CMDistortionReco1",1200,800);
0038   canvas1->Divide(3,2);
0039
0040   //canvas for time plot
0041   TCanvas *canvas=new TCanvas("canvas","CMDistortionReco2",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     //hardcoded numbers from average distortion file's hIntDistortionPosX
0067     int nbinsphi = 30; //when using 35, blank spots at around r = 22 cm, phi just above n below pi
0068     double lowphi = -0.078539819;
0069     double highphi = 6.3617253;
0070     int nbinsr = 30; // when using 35, blank stripe around r = 58 cm
0071     double lowr = 0.0;
0072     double highr = 90.0;
0073
0074     //for forward only
0075
0076     TH2F *hStripesPerBin = new TH2F("hStripesPerBin","CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0077
0078     //phi,r binning
0079     TH2F *hStripesPerBinPhiR = new TH2F("hStripesPerBinPhiR","CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
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     //phi,r binning
0087     TH2F *hCartesianForwardPhiR[3];
0088     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);
0089     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);
0090     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);
0091
0092      TH2F *hCylindricalForwardPhiR[2];
0093     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);
0094     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);
0095
0096     for (int i=0;i<inTree->GetEntries();i++){
0097       inTree->GetEntry(i);
0098
0099       hStripesPerBin->Fill(position->X(),position->Y(),1);
0100
0101       //for phi,r binning
0102       double r = position->Perp();
0103       double phi = position->Phi();
0104
0105       if(position->Phi() < 0.0){
0106     phi = position->Phi() + 2.0*TMath::Pi();
0107       }
0108
0109       hStripesPerBinPhiR->Fill(phi,r,1);
0110
0111       deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
0112       deltaY = (newposition->Y() - position->Y())*(1e4);
0113       deltaZ = (newposition->Z() - position->Z())*(1e4);
0114
0115       deltaR = (newposition->Perp() - position->Perp())*(1e4);
0116       deltaPhi = newposition->DeltaPhi(*position);
0117
0118       hCartesianForward[0]->Fill(position->X(),position->Y(),deltaX);
0119       hCartesianForward[1]->Fill(position->X(),position->Y(),deltaY);
0120       hCartesianForward[2]->Fill(position->X(),position->Y(),deltaZ);
0121
0122       // phi,r binning
0123       hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
0124       hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
0125       hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
0126
0127       hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
0128       hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
0129     }
0130
0131     TH2F *hCartesianAveShift[3];
0132     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);
0133     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);
0134     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);
0135
0136     TH2F *hCylindricalAveShift[2];
0137      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);
0138     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);
0139
0140     // phi,r binning
0141      TH2F *hCartesianAveShiftPhiR[3];
0142     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);
0143     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);
0144     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);
0145
0146     TH2F *hCylindricalAveShiftPhiR[2];
0147      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);
0148     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);
0149
0150     for (int i = 0; i < 3; i ++){
0151       hCartesianAveShift[i]->Divide(hCartesianForward[i],hStripesPerBin);
0152       hCartesianAveShiftPhiR[i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
0153     }
0154
0155     hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
0156     hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
0157
0158     //cyl models from cart coords
0159     for(int i = 0; i < nbins; i++){
0160       double x = low + ((high - low)/(1.0*nbins))*(i+0.5); //center of bin
0161
0162       for(int j = 0; j < nbins; j++){
0163     double y = low + ((high - low)/(1.0*nbins))*(j+0.5); //center of bin
0164
0165     int xbin = hCartesianAveShift[0]->FindBin(x,y);
0166     int ybin = hCartesianAveShift[1]->FindBin(x,y);
0167     double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1e-4); // converts  microns to cm
0168     double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1e-4);
0169
0170     TVector3 shifted, original;
0171     original.SetX(x);
0172     original.SetY(y);
0173     shifted.SetX(x+xaveshift);
0174     shifted.SetY(y+yaveshift);
0175     //x n y above for orig
0176     //shifted is orig + ave shift
0177
0178     double raveshift = (shifted.Perp() - original.Perp())*(1e4);
0179     double phiaveshift = shifted.DeltaPhi(original);
0180
0181     //fill with r from x n y
0182     hCylindricalAveShift[0]->Fill(x,y,raveshift);
0183     hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
0184       }
0185     }
0186
0187     //same range and bins for each coordinate, binned in cm
0188     //hardcoded numbers from average distortion file's hIntDistortionPosX
0189     int nphi = 82;
0190     int nr = 54;
0191     int nz = 82;
0192
0193     double minphi = -0.078539819;
0194     double minr = 18.884615;
0195     //double minz = 5.0;
0196     double minz = -1.3187500;
0197
0198     double maxphi = 6.3617253;
0199     double maxr = 79.115387;
0200     double maxz = 106.81875;
0201
0202     TH3F *hCartesianCMModel[3];
0203     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
0204     hCartesianCMModel[1]=new TH3F("hCMModelY", "CM Model: Y Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0205     hCartesianCMModel[2]=new TH3F("hCMModelZ", "CM Model: Z Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0206
0207     TH3F *hCylindricalCMModel[2];
0208     hCylindricalCMModel[0]=new TH3F("hCMModelRCart", "CM Model: Radial Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0209     hCylindricalCMModel[1]=new TH3F("hCMModelPhiCart", "CM Model: Phi Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0210
0211     //phi,r binning
0212     TH3F *hCartesianCMModelPhiR[3];
0213     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
0214     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);
0215     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);
0216
0217     TH3F *hCylindricalCMModelPhiR[2];
0218     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);
0219     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);
0220
0221     double xshift, yshift, zshift, rshiftcart, phishiftcart;
0222     double xshiftPhiR, yshiftPhiR, zshiftPhiR, rshiftcartPhiR, phishiftcartPhiR;
0223
0224     for(int i = 0; i < nphi; i++){
0225       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0226
0227       for(int j = 0; j < nr; j++){
0228     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0229
0230     double x = r*cos(phi); //cm
0231     double y = r*sin(phi);
0232
0233     for(int k = 0; k < nz; k++){
0234       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0235
0236       xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1e-4);//coordinate of your stripe
0237       yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1e-4);//convert micron to cm
0238       zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1e-4);
0239
0240       rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1e-4);
0241       phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
0242
0243       hCartesianCMModel[0]->Fill(phi,r,z,xshift*(1-z/105.5));
0244       hCartesianCMModel[1]->Fill(phi,r,z,yshift*(1-z/105.5));
0245       hCartesianCMModel[2]->Fill(phi,r,z,zshift*(1-z/105.5));
0246
0247       hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
0248       hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
0249
0250       //phi,r binning
0251       xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);//coordinate of your stripe
0252       yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1e-4);//convert micron to cm
0253       zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1e-4);
0254
0255       rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);
0256       phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
0257
0258       hCartesianCMModelPhiR[0]->Fill(phi,r,z,xshiftPhiR*(1-z/105.5));
0259       hCartesianCMModelPhiR[1]->Fill(phi,r,z,yshiftPhiR*(1-z/105.5));
0260       hCartesianCMModelPhiR[2]->Fill(phi,r,z,zshiftPhiR*(1-z/105.5));
0261
0262       hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
0263       hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
0264     }
0265       }
0266     }
0267
0268     TFile *plots;
0269
0270     plots=TFile::Open(Form("CMModels_Event%d.root",ifile),"RECREATE");
0271     hStripesPerBin->Write();
0272
0273     for(int i = 0; i < 3; i++){
0274       hCartesianForward[i]->Write();
0275       hCartesianAveShift[i]->Write();
0276       hCartesianCMModel[i]->Write();
0277
0278       hCartesianForwardPhiR[i]->Write();
0279       hCartesianAveShiftPhiR[i]->Write();
0280       hCartesianCMModelPhiR[i]->Write();
0281     }
0282
0283     for(int i = 0; i < 2; i++){
0284       hCylindricalAveShift[i]->Write();
0285       hCylindricalCMModel[i]->Write();
0286
0287       hCylindricalForwardPhiR[i]->Write();
0288       hCylindricalAveShiftPhiR[i]->Write();
0289       hCylindricalCMModelPhiR[i]->Write();
0290     }
0291
0292     plots->Close();
0293
0294
0295     //call to TTime after outputting TH3Fs
0296     now=gSystem->Now();
0297     unsigned long after = now;
0298
0299     hTimePerEvent->Fill(after-before);
0300
0301     //to check histograms
0302     for (int i = 0; i < 3; i++){
0303       hCartesianForward[i]->SetStats(0);
0304       hCartesianForwardPhiR[i]->SetStats(0);
0305     }
0306
0307     hStripesPerBin->SetStats(0);
0308     canvas1->cd(1);
0309     hStripesPerBinPhiR->Draw("colz");
0310     canvas1->cd(2);
0311     hCartesianForward[1]->Draw("colz");
0312     canvas1->cd(3);
0313     hCartesianForward[2]->Draw("colz");
0314     canvas1->cd(4);
0315     hCartesianForwardPhiR[0]->Draw("colz");
0316     canvas1->cd(5);
0317     hCartesianForwardPhiR[1]->Draw("colz");
0318     canvas1->cd(6);
0319     hCartesianForwardPhiR[2]->Draw("colz");
0320
0321     if(ifile == 0){
0322       canvas1->Print("CMDistortionReco1.pdf(","pdf");
0323     } else if (ifile == nEvents - 1){
0324       canvas1->Print("CMDistortionReco1.pdf)","pdf");
0325     } else {
0326       canvas1->Print("CMDistortionReco1.pdf","pdf");
0327     }
0328
0329   }
0330
0331   canvas->cd();
0332   hTimePerEvent->Draw();
0333   canvas->Print("CMDistortionReco2.pdf","pdf");
0334
0335   return 0;
0336 }