tpc/cm/CMDistortionReco.C

0001 // step 2 with phi,r coords
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 CMDistortionReco(int nMaxEvents = -1) {
0022   int nbins = 35;
0023   double low = -80.0;
0024   double high = 80.0;
0025   double deltaX;
0026   double deltaY;
0027   double deltaZ;
0028   double deltaR;
0029   double deltaPhi;
0030   int nEvents;
0031
0032   //take in events
0033   const char * inputpattern="/sphenix/u/skurdi/CMCalibration/cmDistHitsTree_Event*.root";
0034
0035   //find all files that match the input string (includes wildcards)
0036   TFileCollection *filelist=new TFileCollection();
0037   filelist->Add(inputpattern);
0038   TString sourcefilename;
0039
0040   //how many events
0041   if(nMaxEvents >= 0 && nMaxEvents<filelist->GetNFiles())
0042   {
0043     nEvents=nMaxEvents;
0044   }
0045   else
0046   {
0047     nEvents= filelist->GetNFiles();
0048   }
0049   //  nEvents = 2;
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 *hStripesPerBin = new TH2F("hStripesPerBin","CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0092
0093     //phi,r binning
0094     TH2F *hStripesPerBinPhiR = new TH2F("hStripesPerBinPhiR","CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
0095
0096     TH2F *hCartesianForward[3];
0097     hCartesianForward[0] = new TH2F("hForwardX","X Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0098     hCartesianForward[1] = new TH2F("hForwardY","Y Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0099     hCartesianForward[2] = new TH2F("hForwardZ","Z Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0100
0101     //phi,r binning
0102     TH2F *hCartesianForwardPhiR[3];
0103     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);
0104     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);
0105     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);
0106
0107      TH2F *hCylindricalForwardPhiR[2];
0108     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);
0109     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);
0110
0111     for (int i=0;i<inTree->GetEntries();i++){
0112       inTree->GetEntry(i);
0113
0114       hStripesPerBin->Fill(position->X(),position->Y(),1);
0115
0116       //for phi,r binning
0117       double r = position->Perp();
0118       double phi = position->Phi();
0119
0120       if(position->Phi() < 0.0){
0121     phi = position->Phi() + 2.0*TMath::Pi();
0122       }
0123
0124       hStripesPerBinPhiR->Fill(phi,r,1);
0125
0126       deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
0127       deltaY = (newposition->Y() - position->Y())*(1e4);
0128       deltaZ = (newposition->Z() - position->Z())*(1e4);
0129
0130       deltaR = (newposition->Perp() - position->Perp())*(1e4);
0131       deltaPhi = newposition->DeltaPhi(*position);
0132
0133       hCartesianForward[0]->Fill(position->X(),position->Y(),deltaX);
0134       hCartesianForward[1]->Fill(position->X(),position->Y(),deltaY);
0135       hCartesianForward[2]->Fill(position->X(),position->Y(),deltaZ);
0136
0137       // phi,r binning
0138       hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
0139       hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
0140       hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
0141
0142       hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
0143       hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
0144     }
0145
0146     TH2F *hCartesianAveShift[3];
0147     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);
0148     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);
0149     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);
0150
0151     TH2F *hCylindricalAveShift[2];
0152      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);
0153     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);
0154
0155     // phi,r binning
0156      TH2F *hCartesianAveShiftPhiR[3];
0157     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);
0158     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);
0159     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);
0160
0161     TH2F *hCylindricalAveShiftPhiR[2];
0162      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);
0163     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);
0164
0165     for (int i = 0; i < 3; i ++){
0166       hCartesianAveShift[i]->Divide(hCartesianForward[i],hStripesPerBin);
0167       hCartesianAveShiftPhiR[i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
0168     }
0169
0170     hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
0171     hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
0172
0173     //cyl models from cart coords
0174     for(int i = 0; i < nbins; i++){
0175       double x = low + ((high - low)/(1.0*nbins))*(i+0.5); //center of bin
0176
0177       for(int j = 0; j < nbins; j++){
0178     double y = low + ((high - low)/(1.0*nbins))*(j+0.5); //center of bin
0179
0180     int xbin = hCartesianAveShift[0]->FindBin(x,y);
0181     int ybin = hCartesianAveShift[1]->FindBin(x,y);
0182     double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1e-4); // converts  microns to cm
0183     double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1e-4);
0184
0185     TVector3 shifted;
0186     TVector3 original;
0187     original.SetX(x);
0188     original.SetY(y);
0189     shifted.SetX(x+xaveshift);
0190     shifted.SetY(y+yaveshift);
0191     //x n y above for orig
0192     //shifted is orig + ave shift
0193
0194     double raveshift = (shifted.Perp() - original.Perp())*(1e4);
0195     double phiaveshift = shifted.DeltaPhi(original);
0196
0197     //fill with r from x n y
0198     hCylindricalAveShift[0]->Fill(x,y,raveshift);
0199     hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
0200       }
0201     }
0202
0203     //same range and bins for each coordinate, binned in cm
0204     //hardcoded numbers from average distortion file's hIntDistortionPosX
0205     int nphi = 82;
0206     int nr = 54;
0207     int nz = 82;
0208
0209     double minphi = -0.078539819;
0210     double minr = 18.884615;
0211     //double minz = 5.0;
0212     double minz = -1.3187500;
0213
0214     double maxphi = 6.3617253;
0215     double maxr = 79.115387;
0216     double maxz = 106.81875;
0217
0218     TH3F *hCartesianCMModel[3];
0219     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
0220     hCartesianCMModel[1]=new TH3F("hCMModelY", "CM Model: Y Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0221     hCartesianCMModel[2]=new TH3F("hCMModelZ", "CM Model: Z Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0222
0223     TH3F *hCylindricalCMModel[2];
0224     hCylindricalCMModel[0]=new TH3F("hCMModelRCart", "CM Model: Radial Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0225     hCylindricalCMModel[1]=new TH3F("hCMModelPhiCart", "CM Model: Phi Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
0226
0227     //phi,r binning
0228     TH3F *hCartesianCMModelPhiR[3];
0229     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
0230     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);
0231     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);
0232
0233     TH3F *hCylindricalCMModelPhiR[2];
0234     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);
0235     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);
0236
0237     double xshift;
0238     double yshift;
0239     double zshift;
0240     double rshiftcart;
0241     double phishiftcart;
0242     double xshiftPhiR;
0243     double yshiftPhiR;
0244     double zshiftPhiR;
0245     double rshiftcartPhiR;
0246     double phishiftcartPhiR;
0247
0248     for(int i = 0; i < nphi; i++){
0249       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0250
0251       for(int j = 0; j < nr; j++){
0252     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0253
0254     double x = r*cos(phi); //cm
0255     double y = r*sin(phi);
0256
0257     for(int k = 0; k < nz; k++){
0258       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0259
0260       xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1e-4);//coordinate of your stripe
0261       yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1e-4);//convert micron to cm
0262       zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1e-4);
0263
0264       rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1e-4);
0265       phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
0266
0267       hCartesianCMModel[0]->Fill(phi,r,z,xshift*(1-z/105.5));
0268       hCartesianCMModel[1]->Fill(phi,r,z,yshift*(1-z/105.5));
0269       hCartesianCMModel[2]->Fill(phi,r,z,zshift*(1-z/105.5));
0270
0271       hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
0272       hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
0273
0274       //phi,r binning
0275       xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);//coordinate of your stripe
0276       yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1e-4);//convert micron to cm
0277       zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1e-4);
0278
0279       rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);
0280       phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
0281
0282       hCartesianCMModelPhiR[0]->Fill(phi,r,z,xshiftPhiR*(1-z/105.5));
0283       hCartesianCMModelPhiR[1]->Fill(phi,r,z,yshiftPhiR*(1-z/105.5));
0284       hCartesianCMModelPhiR[2]->Fill(phi,r,z,zshiftPhiR*(1-z/105.5));
0285
0286       hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
0287       hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
0288     }
0289       }
0290     }
0291
0292     TFile *plots;
0293
0294     plots=TFile::Open(std::format("CMModels_Event{}.root",ifile).c_str(),"RECREATE");
0295     hStripesPerBin->Write();
0296
0297     for(int i = 0; i < 3; i++){
0298       hCartesianForward[i]->Write();
0299       hCartesianAveShift[i]->Write();
0300       hCartesianCMModel[i]->Write();
0301
0302       hCartesianForwardPhiR[i]->Write();
0303       hCartesianAveShiftPhiR[i]->Write();
0304       hCartesianCMModelPhiR[i]->Write();
0305     }
0306
0307     for(int i = 0; i < 2; i++){
0308       hCylindricalAveShift[i]->Write();
0309       hCylindricalCMModel[i]->Write();
0310
0311       hCylindricalForwardPhiR[i]->Write();
0312       hCylindricalAveShiftPhiR[i]->Write();
0313       hCylindricalCMModelPhiR[i]->Write();
0314     }
0315
0316     plots->Close();
0317
0318
0319     //call to TTime after outputting TH3Fs
0320     now=gSystem->Now();
0321     unsigned long after = now;
0322
0323     hTimePerEvent->Fill(after-before);
0324
0325     //to check histograms
0326     for (int i = 0; i < 3; i++){
0327       hCartesianForward[i]->SetStats(false);
0328       hCartesianForwardPhiR[i]->SetStats(false);
0329     }
0330
0331     hStripesPerBin->SetStats(false);
0332     canvas1->cd(1);
0333     hStripesPerBinPhiR->Draw("colz");
0334     canvas1->cd(2);
0335     hCartesianForward[1]->Draw("colz");
0336     canvas1->cd(3);
0337     hCartesianForward[2]->Draw("colz");
0338     canvas1->cd(4);
0339     hCartesianForwardPhiR[0]->Draw("colz");
0340     canvas1->cd(5);
0341     hCartesianForwardPhiR[1]->Draw("colz");
0342     canvas1->cd(6);
0343     hCartesianForwardPhiR[2]->Draw("colz");
0344
0345     if(ifile == 0){
0346       canvas1->Print("CMDistortionReco1.pdf(","pdf");
0347     } else if (ifile == nEvents - 1){
0348       canvas1->Print("CMDistortionReco1.pdf)","pdf");
0349     } else {
0350       canvas1->Print("CMDistortionReco1.pdf","pdf");
0351     }
0352
0353   }
0354
0355   canvas->cd();
0356   hTimePerEvent->Draw();
0357   canvas->Print("CMDistortionReco2.pdf","pdf");
0358
0359   return 0;
0360 }