tpc/cm/CMDistortionRecoPhiR.C

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