SteppingHelixPropagator/test/MuProp.C

0001 #define MuProp_cxx
0002 #include "MuProp.h"
0003 #include <TStyle.h>
0004 #include <TCanvas.h>
0005
0006 #include "TVector3.h"
0007 #include "TLorentzVector.h"
0008 #include <iostream>
0009
0010 Double_t MUON_MASS = 0.1056584;
0011
0012 void MuProp::Loop(Double_t maxEloss)
0013 {
0014    if (fChain == 0) return;
0015
0016    Long64_t nentries = fChain->GetEntriesFast();
0017
0018    Long64_t nbytes = 0, nb = 0;
0019    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0020       Long64_t ientry = LoadTree(jentry);
0021       if (ientry < 0) break;
0022       nb = fChain->GetEntry(jentry);   nbytes += nb;
0023
0024       //      Float_t p0 = sqrt(p3R[0][0]**2+p3R[0][1]**2+p3R[0][2]**2);
0025       TVector3 p3_0; p3_0.SetXYZ(p3R[0][0], p3R[0][1], p3R[0][2]);
0026
0027       for (Int_t iPoint = 1; iPoint < nPoints; iPoint++){
0028     TVector3 r3Exp; r3Exp.SetXYZ(r3R[iPoint][0], r3R[iPoint][1], r3R[iPoint][2]);
0029     TVector3 r3Sim; r3Sim.SetXYZ( r3[iPoint][0],  r3[iPoint][1],  r3[iPoint][2]);
0030
0031     TLorentzVector p3Exp; p3Exp.SetXYZM(p3R[iPoint][0], p3R[iPoint][1], p3R[iPoint][2], MUON_MASS);
0032     TLorentzVector p3Sim; p3Sim.SetXYZM( p3[iPoint][0],  p3[iPoint][1],  p3[iPoint][2], MUON_MASS);
0033     if (fabs(p3Exp.E() - p3Sim.E())>maxEloss) continue;
0034
0035     Float_t rExp = r3Exp.Perp();
0036     Float_t dX = rExp*r3Exp.DeltaPhi(r3Sim);
0037
0038     Float_t ddXdX = (r3Exp.X()*dX + rExp*r3Exp.Y())/rExp/rExp;
0039     Float_t ddXdY = (r3Exp.Y()*dX - rExp*r3Exp.X())/rExp/rExp;
0040     Float_t dXSigma = ddXdX*ddXdX*covFlat[iPoint][0]
0041       + ddXdY*ddXdY*covFlat[iPoint][2]
0042       + 2.*ddXdX*ddXdY*covFlat[iPoint][1];
0043     dXSigma = sqrt(dXSigma);
0044
0045     Float_t dXPull = dX/dXSigma;
0046
0047     Int_t idF[2];
0048     Int_t dSubDMask = (0xF<<28) | (0x7<<25);
0049     //wheel*station
0050     Int_t dtWSMask = dSubDMask | (0x3F<<19);
0051     //endcap*station*ring
0052     Int_t cscESRMask = dSubDMask | (0xFF<<10);
0053     idF[0] = id[iPoint] & dSubDMask;
0054     Int_t det =  (idF[0]>>28) & 0xF;
0055     Int_t subDet =  (idF[0]>>25) & 0x7;
0056     //muon detector only
0057     if ( det != 2) continue;
0058     //DT and CSC only
0059     if ( subDet != 1 && subDet != 2 ) continue;
0060
0061     if ( subDet == 1 ){
0062       idF[1] = id[iPoint] & dtWSMask;
0063     } else if ( subDet == 2){
0064       idF[1] = id[iPoint] & cscESRMask;
0065     }
0066
0067     for (Int_t iF = 0; iF < 2; iF++){
0068       if (dX_mh1[idF[iF]] == 0){
0069         std::string hName = Form("dX_%X", idF[iF]);
0070         dX_mh1[idF[iF]] = new TH1F(hName.c_str(), hName.c_str(), 100, -50, 50);
0071       }
0072       dX_mh1[idF[iF]]->Fill(dX);
0073
0074
0075       if (dXPull_mh1[idF[iF]] == 0){
0076         std::string hName = Form("dXPull_%X", idF[iF]);
0077         dXPull_mh1[idF[iF]] = new TH1F(hName.c_str(), hName.c_str(), 100, -5, 5);
0078       }
0079       dXPull_mh1[idF[iF]]->Fill(dXPull);
0080
0081     }
0082       }
0083
0084    }
0085 }