/*
$Log$
+Revision 1.11 2002/03/08 17:25:36 cussonno
+Update absorber energy loss and Branson corrections : simplified functions
+BransonCorrection and TotalMomentumEnergyLoss.
+
Revision 1.10 2001/04/25 14:50:42 gosset
Corrections to violations of coding conventions
Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
// Extrapolates track parameters upstream to the "Z" end of the front absorber
- ExtrapToZ(zAbsorber);
+ ExtrapToZ(zAbsorber); // !!!
// Makes Branson correction (multiple scattering + energy loss)
BransonCorrection();
+ // Makes a simple magnetic field correction through the absorber
+ FieldCorrection(zAbsorber);
}
Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
Int_t sign;
static Bool_t first = kTRUE;
- static Double_t zBP1, zBP2, rLimit, zEndAbsorber;
+ static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber;
// zBP1 for outer part and zBP2 for inner part (only at the first call)
if (first) {
first = kFALSE;
zEndAbsorber = 503;
- rLimit = zEndAbsorber * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
- zBP1 = 450;
+ thetaLimit = 3.0 * (TMath::Pi()) / 180.;
+ rLimit = zEndAbsorber * TMath::Tan(thetaLimit);
+ zBP1 = 450; // values close to those calculated with EvalAbso.C
zBP2 = 480;
}
yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
// new parameters after Branson and energy loss corrections
- pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
- pX = pZ * xBP / zBP;
- pY = pZ * yBP / zBP;
+// Float_t zSmear = zBP - gRandom->Gaus(0.,2.); // !!! possible smearing of Z vertex position
+ Float_t zSmear = zBP;
+
+ pZ = pTotal * zSmear / TMath::Sqrt(xBP * xBP + yBP * yBP + zSmear * zSmear);
+ pX = pZ * xBP / zSmear;
+ pY = pZ * yBP / zSmear;
fBendingSlope = pY / pZ;
fNonBendingSlope = pX / pZ;
pT = TMath::Sqrt(pX * pX + pY * pY);
theta = TMath::ATan2(pT, pZ);
- pTotal =
- TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
+ pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta);
fInverseBendingMomentum = (sign / pTotal) *
TMath::Sqrt(1.0 +
fNonBendingCoor = 0;
fZ= 0;
}
+
//__________________________________________________________________________
-Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
+Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta)
{
// Returns the total momentum corrected from energy loss in the front absorber
// One can use the macros MUONTestAbso.C and DrawTestAbso.C
// to test this correction.
+ // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002)
Double_t deltaP, pTotalCorrected;
- Double_t radiusEndAbsorber2 =
- xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
- // Parametrization to be redone according to change of absorber material ????
+ // Parametrization to be redone according to change of absorber material ????
// See remark in function BransonCorrection !!!!
// The name is not so good, and there are many arguments !!!!
- if (radiusEndAbsorber2 < rLimit * rLimit) {
- if (pTotal < 15) {
- deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
+ if (theta < thetaLimit ) {
+ if (pTotal < 20) {
+ deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal;
} else {
- deltaP = 3.0643 + 0.01346 *pTotal;
+ deltaP = 3.0714 + 0.011767 *pTotal;
}
- deltaP = 0.63 * deltaP; // !!!! changes in the absorber composition ????
} else {
- if (pTotal < 15) {
- deltaP = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
+ if (pTotal < 20) {
+ deltaP = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal;
} else {
- deltaP = 2.407 + 0.00702 * pTotal;
+ deltaP = 2.6069 + 0.0051705 * pTotal;
}
- deltaP = 0.67 * deltaP; // !!!! changes in the absorber composition ????
}
pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
return pTotalCorrected;
}
+ //__________________________________________________________________________
+void AliMUONTrackParam::FieldCorrection(Double_t Z)
+{
+ //
+ // Correction of the effect of the magnetic field in the absorber
+ // Assume a constant field along Z axis.
+
+ Float_t b[3],x[3];
+ Double_t bZ;
+ Double_t pYZ,pX,pY,pZ,pT;
+ Double_t pXNew,pYNew;
+ Double_t c;
+
+ pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+ c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge
+
+ pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
+ pX = pZ * fNonBendingSlope;
+ pY = pZ * fBendingSlope;
+ pT = TMath::Sqrt(pX*pX+pY*pY);
+
+ if (pZ <= 0) return;
+ x[2] = Z/2;
+ x[0] = x[2]*fNonBendingSlope;
+ x[1] = x[2]*fBendingSlope;
+
+ // Take magn. field value at position x.
+ gAlice->Field()->Field(x, b);
+ bZ = b[2];
+
+ // Transverse momentum rotation
+ // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ.
+ Double_t phiShift = c*0.436*0.0003*bZ*Z/pZ;
+
+ // Rotate momentum around Z axis.
+ pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift);
+ pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift);
+
+ fBendingSlope = pYNew / pZ;
+ fNonBendingSlope = pXNew / pZ;
+
+ fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ);
+
+}
-void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
+void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0,Int_t testSingle=0)
{
// Useful to test the absorber correction in the reconstruction procedure
// (mean energy loss + Branson correction)
// Create some histograms
- TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 120, 8., 11.);
- TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons r < rLimit", 100, -10., 10.);
- TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons r < rLimit", 100, 0., 0.5);
+ TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 1200, 0., 12.);
+// TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 100, 9., 10.5);
+ TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons theta < 3 deg", 100, -10., 10.);
+ TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons theta < 3 deg", 100, 0., 0.5);
//
- TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons r > rLimit", 100, -10., 10.);
- TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons r > rLimit", 100, 0., 0.5);
+ TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons theta > 3 deg", 100, -10., 10.);
+ TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons theta > 3 deg", 100, 0., 0.5);
+
+ //
+ TH1F *hRap = new TH1F("hRap"," Muon Rapidity gen.",20,2.5,4);
+ //
+
+ Int_t nBinProf=19;
+ Float_t xProfMin=5;
+ Float_t xProfMax=290;
+ char titleHist[50];
+ char numberHist[10];
+ TH1F *DeltaP1[100],*DeltaP2[100];
+ TH1F *DeltaAngleX1[100],*DeltaAngleX2[100];
+ TH1F *DeltaAngleY1[100],*DeltaAngleY2[100];
+ TH1F *hP1 = new TH1F("hP1"," P theta < 3 deg ",nBinProf,xProfMin,xProfMax);
+ TProfile *hProfDeltaPvsP1 = new TProfile("hProfDeltaPvsP1"," delta P vs P theta < 3 deg ",nBinProf,xProfMin,xProfMax,-50,50,"s");
+ TH2F *h2DeltaPvsP1 = new TH2F("h2DeltaPvsP1"," delta P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
+ TH1F *hSigmaPvsP1 = new TH1F("hSigmaPvsP1"," deltaP/P vs P theta < 3 deg",nBinProf,xProfMin,xProfMax);
+ for (Int_t iHist=0; iHist < nBinProf; iHist++){
+ sprintf(titleHist,"%s","hDelta1P");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaP1[iHist] = new TH1F(titleHist," deltaP theta < 3 deg ",400,-50,50);
+
+ sprintf(titleHist,"%s","hDelta1AngleX");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaAngleX1[iHist] = new TH1F(titleHist," delta Angle X theta < 3 deg ",400,-0.05,0.05);
+
+ sprintf(titleHist,"%s","hDelta1AngleY");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaAngleY1[iHist] = new TH1F(titleHist," delta Angle Y theta < 3 deg ",400,-0.05,0.05);
+ }
+ TH1F *hSigGausPvsP1 = new TH1F("hSigGausPvsP1"," deltaP/P vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
+ TH1F *hSigGausAngleXvsP1 = new TH1F("hSigGausAngleXvsP1"," delta theta X vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
+ TH1F *hSigGausAngleYvsP1 = new TH1F("hSigGausAngleYvsP1"," delta theta Y vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
+
+ TH1F *hP2 = new TH1F("hP2"," P theta > 3 deg ",nBinProf,xProfMin,xProfMax);
+ TProfile *hProfDeltaPvsP2 = new TProfile("hProfDeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,-50,50,"s");
+ TH2F *h2DeltaPvsP2 = new TH2F("h2DeltaPvsP2"," delta P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
+ TH1F *hSigmaPvsP2 = new TH1F("hSigmaPvsP2"," deltaP/P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax);
+ for (Int_t iHist=0; iHist < nBinProf; iHist++){
+ sprintf(titleHist,"%s","hDelta2P");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaP2[iHist] = new TH1F(titleHist," deltaP theta > 3 deg ",400,-50,50);
+
+ sprintf(titleHist,"%s","hDelta2AngleX");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaAngleX2[iHist] = new TH1F(titleHist," delta Angle X theta > 3 deg ",400,-0.05,0.05);
+
+ sprintf(titleHist,"%s","hDelta2AngleY");
+ sprintf(numberHist,"%d",iHist+1);
+ strcat(titleHist,numberHist);
+ DeltaAngleY2[iHist] = new TH1F(titleHist," delta Angle Y theta > 3 deg ",400,-0.05,0.05);
+ }
+ TH1F *hSigGausPvsP2 = new TH1F("hSigGausPvsP2"," deltaP/P vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
+
+ TH1F *hSigGausAngleXvsP2 = new TH1F("hSigGausAngleXvsP2"," delta theta X vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
+ TH1F *hSigGausAngleYvsP2 = new TH1F("hSigGausAngleYvsP2"," delta theta Y vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
+
+ TProfile *hDeltaPxvsPhi = new TProfile("hDeltaPxvsPhi"," delta Px vs Phi",45,-180,180,-1,1,"s");
+ TProfile *hDeltaPyvsPhi = new TProfile("hDeltaPyvsPhi"," delta Py vs Phi",45,-180,180,-1,1,"s");
+ TProfile *hDeltaPhivsPz = new TProfile("hDeltaPhivsPz"," delta phi vs pZ",25,0,100,-4,4);
+
// Define constant parameters
Float_t massMuon = 0.105658389; // muon mass
Float_t massUpsilon = 9.46037; // Upsilon mass
- Double_t PI = 3.141592654;
Double_t zEndAbso = 503; // z position of the end of the absorber
- Double_t rLimit = zEndAbso * TMath::Tan(3*PI/180); // 3 deg. limit (different absorber composition)
+ Double_t rLimit = zEndAbso * TMath::Tan(3*TMath::Pi()/180); // 3 deg. limit (different absorber composition)
Double_t printLevel = 0;
if (printLevel >= 1)
Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex
- for (Int_t nev=0; nev<= evNumber2; nev++) { // loop over events
- if (printLevel >= 0) {
- cout << " "<<endl;
+ for (Int_t nev=evNumber1; nev<= evNumber2; nev++) { // loop over events
+ if (printLevel >= 0 && nev%100 == 0) {
cout << " **** Event # " << nev <<endl;
}
if (nev < evNumber1) continue;
if (nparticles <= 0) return;
- Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2]; // parameters of generated particles at the vertex
-
+ Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2],RapGen[2]; // parameters of generated particles at the vertex
+ ThetaGen[0] = ThetaGen[1] = 0;
+ RapGen[0] = RapGen[1] = 0;
+
for (int iPart = 0; iPart < nparticles ; iPart++) { // loop over particles
TParticle *particle;
PyGen[0] = particle->Py();
PzGen[0] = particle->Pz();
PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]);
- ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0])*180/PI;
+ ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0]);
+ RapGen[0] = rapParticle(PxGen[0],PyGen[0],PzGen[0],massMuon);
if (printLevel >= 1) {
cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
}
PyGen[1] = particle->Py();
PzGen[1] = particle->Pz();
PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]);
- ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1])*180/PI;
+ ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1]);
+ RapGen[1] = rapParticle(PxGen[1],PyGen[1],PzGen[1],massMuon);
if (printLevel >= 1)
cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
}
- }
+ } // end loop over particles
AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
TTree *TH = gAlice->TreeH();
Int_t ntracks =(Int_t) TH->GetEntries();
-
+
+ Bool_t hitChamberMUPlus[10],hitChamberMUMinus[10];
+ Bool_t hitStationMUPlus[5],hitStationMUMinus[5];
+ for (Int_t iCh=0; iCh < 10; iCh++){
+ hitChamberMUPlus[iCh] = kFALSE;
+ hitChamberMUMinus[iCh] = kFALSE;
+ }
+ for (Int_t iSt=0; iSt < 5; iSt++){
+ hitStationMUPlus[iSt] = kFALSE;
+ hitStationMUMinus[iSt] = kFALSE;
+ }
+
+ Int_t nHitChamber=0;
+
for (Int_t i = 0; i < 2; i++){
pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
}
{
Int_t nch = mHit->Chamber(); // chamber number
- if (nch != 1) continue;
Int_t ipart = mHit->Particle();
Double_t x = mHit->X(); // x-pos of hit in chamber #1
Double_t y = mHit->Y(); // y-pos of hit in chamber #1
Double_t z = mHit->Z(); // z-pos of hit in chamber #1
+
+ Int_t iSt = (nch-1)/2;
+ if ((nch-1) < 10) {
+ if (ipart == kMuonPlus){
+ hitChamberMUPlus[nch-1] = kTRUE;
+ hitStationMUPlus[iSt] = kTRUE;
+ }
+ if (ipart == kMuonMinus){
+ hitChamberMUMinus[nch-1] = kTRUE;
+ hitStationMUMinus[iSt] = kTRUE;
+ }
+ nHitChamber++;
+ }
+ if (nch != 1) continue;
if (ipart == kMuonPlus || ipart == kMuonMinus) {
Double_t px0=mHit->Px();
// create an AliMUONTrackParam object in chamber #1
AliMUONTrackParam *trackParam = new AliMUONTrackParam();
- Double_t bendingMometum = TMath::Sqrt(pz0*pz0+px0*px0);
+ Int_t sign = -TMath::Sign(1.0,ipart);
+ Double_t bendingMometum = sign* TMath::Sqrt(pz0*pz0+py0*py0); // Bug px -> py !!!
Double_t inverseBendingMomentum = 1/bendingMometum;
trackParam->SetBendingCoor(y);
cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
trackParam->Dump();
}
+ if (!testSingle)
trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
if (printLevel >= 1) {
cout <<" **** after extrap " <<endl;
pX[iMuon] = pZ[iMuon] * nonBendingSlope;
pY[iMuon] = pZ[iMuon] * bendingSlope;
pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]);
- theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon])*180/PI;
+ theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon]);
radius[iMuon] = TMath::Sqrt(x*x+y*y);
if (printLevel >= 1)
} // if MuonPlus or MuonMinus
} // loop over MUON hits
} // if MUON
- } // loop over tracks
+ } // loop over tracks
+
+
+ Bool_t goodMUPlus = kTRUE;
+ Bool_t goodMUMinus = kTRUE;
+ for (Int_t iCh=0; iCh < 1; iCh++) { // !!!! 10 -> 1
+ if (!hitChamberMUPlus[iCh]) {
+ goodMUPlus = kFALSE;
+ printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
+ }
+ if (!hitChamberMUMinus[iCh]) {
+ goodMUMinus = kFALSE;
+ printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
+ }
+ }
+
+// for (Int_t iSt=0; iSt < 5; iSt++) {
+// if (!hitStationMUPlus[iSt]) {
+// goodMUPlus = kFALSE;
+// // printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
+// }
+// if (!hitStationMUMinus[iSt]) {
+// goodMUMinus = kFALSE;
+// // printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
+// }
+// }
if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist
Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
hInvMassRes->Fill(massResonance);
- for (Int_t i = 0; i<2; i++){
- Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
- cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
- Double_t deltaAngle = TMath::ACos(cosAngle)*180/PI;
- if (radius[i] < rLimit) {
- hDeltaP1->Fill(pTot[i]-PTotGen[i]);
- hDeltaAngle1->Fill(deltaAngle);
- }else{
- hDeltaP2->Fill(pTot[i]-PTotGen[i]);
- hDeltaAngle2->Fill(deltaAngle);
+ }
+
+
+ for (Int_t i = 0; i<2; i++){ // !!! 1 -> 2
+ if (i == 0 && !goodMUPlus) continue;
+ if (i == 1 && !goodMUMinus) continue;
+ if (pX[i] == 0) continue;
+
+ hRap->Fill(RapGen[i]);
+
+ Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
+ cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
+ Double_t deltaAngle = TMath::ACos(cosAngle)*180/TMath::Pi();
+ Double_t deltaPNorm = 0;
+ if (testSingle) {
+ deltaPNorm = (PTotGen[i]-pTot[i])* TMath::Cos(ThetaGen[i]);
+ } else {
+ deltaPNorm = PTotGen[i]-pTot[i];
+ }
+
+ Float_t thetaGX = TMath::ATan2(PxGen[i],PzGen[i]);
+ Float_t thetaRX = TMath::ATan2(pX[i],pZ[i]);
+ Float_t deltaAngleX = thetaRX - thetaGX;
+
+ Float_t thetaGY = TMath::ATan2(PyGen[i],PzGen[i]);
+ Float_t thetaRY = TMath::ATan2(pY[i],pZ[i]);
+ Float_t deltaAngleY = thetaRY - thetaGY;
+
+ Float_t phiG = TMath::ATan2(PyGen[i],PxGen[i])*180/TMath::Pi();
+ Float_t phiR = TMath::ATan2(pY[i],pX[i])*180/TMath::Pi();
+
+ hDeltaPxvsPhi->Fill(phiG,PxGen[i]-pX[i]);
+ hDeltaPyvsPhi->Fill(phiG,PyGen[i]-pY[i]);
+ hDeltaPhivsPz->Fill(PzGen[i],phiR-phiG);
+
+ Int_t iHist=0;
+ if (ThetaGen[i] < (3*TMath::Pi()/180)) {
+ hDeltaP1->Fill(pTot[i]-PTotGen[i]);
+ hDeltaAngle1->Fill(deltaAngle);
+ hP1->Fill(PTotGen[i]);
+ hProfDeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
+ h2DeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
+ iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
+ if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) {
+ DeltaP1[iHist]->Fill(deltaPNorm);
+ DeltaAngleX1[iHist]->Fill(deltaAngleX);
+ DeltaAngleY1[iHist]->Fill(deltaAngleY);
+ }
+ } else {
+ hDeltaP2->Fill(pTot[i]-PTotGen[i]);
+ hDeltaAngle2->Fill(deltaAngle);
+ hP2->Fill(PTotGen[i]);
+ hProfDeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
+ h2DeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
+ iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
+ if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) {
+ DeltaP2[iHist]->Fill(deltaPNorm);
+ DeltaAngleX2[iHist]->Fill(deltaAngleX);
+ DeltaAngleY2[iHist]->Fill(deltaAngleY);
}
}
}
} // loop over event
+
+
+ Float_t sigmaP = 0;
+ Float_t xBin = 0;
+ for (Int_t iBin = 1; iBin <= nBinProf; iBin++){
+ sigmaP = hProfDeltaPvsP1->GetBinError(iBin);
+ xBin = hProfDeltaPvsP1->GetBinCenter(iBin);
+ hSigmaPvsP1->SetBinContent(iBin,sigmaP/xBin);
+ sigmaP = hProfDeltaPvsP2->GetBinError(iBin);
+ xBin = hProfDeltaPvsP2->GetBinCenter(iBin);
+ hSigmaPvsP2->SetBinContent(iBin,sigmaP/xBin);
+ }
+
+ TF1 *g1= new TF1("g1","gaus",-25,25) ;
+ TF1 *g2= new TF1("g2","gaus",-25,25) ;
+ Float_t sigmaPGaus;
+ Float_t xBinGaus;
+ for (Int_t iHist=0; iHist < nBinProf; iHist++){
+ DeltaP1[iHist]->Fit("g1","RQ");
+ sigmaPGaus = g1->GetParameter(2);
+ printf(" ** 1 ** iHist= %d , sigmaPGaus = %f \n",iHist,sigmaPGaus);
+ xBinGaus = hSigGausPvsP1->GetBinCenter(iHist+1);
+ hSigGausPvsP1->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
+
+ DeltaP2[iHist]->Fit("g2","RQ");
+ sigmaPGaus = g2->GetParameter(2);
+ printf(" ** 2 ** iHist= %d , sigmaGaus = %f \n",iHist,sigmaPGaus);
+ xBinGaus = hSigGausPvsP2->GetBinCenter(iHist+1);
+ hSigGausPvsP2->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
+ }
+
+ // Angles
+ TF1 *g3= new TF1("g3","gaus") ;
+ TF1 *g4= new TF1("g4","gaus") ;
+ TF1 *g5= new TF1("g5","gaus") ;
+ TF1 *g6= new TF1("g6","gaus") ;
+ Float_t sigmaAngleGaus,limitGaus,errorSigma;
+ Float_t nSigFit = 3;
+ for (Int_t iHist=0; iHist < nBinProf; iHist++){
+ limitGaus = nSigFit * (DeltaAngleX1[iHist]->GetRMS());
+ g3->SetRange(-limitGaus,limitGaus);
+ DeltaAngleX1[iHist]->Fit("g3","RQ");
+ sigmaAngleGaus = g3->GetParameter(2);
+ errorSigma = g3->GetParError(2);
+ printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
+ hSigGausAngleXvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
+ hSigGausAngleXvsP1->SetBinError(iHist+1,errorSigma);
+
+ limitGaus = nSigFit * (DeltaAngleY1[iHist]->GetRMS());
+ g4->SetRange(-limitGaus,limitGaus);
+ DeltaAngleY1[iHist]->Fit("g4","RQ");
+ sigmaAngleGaus = g4->GetParameter(2);
+ errorSigma = g4->GetParError(2);
+ printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
+ hSigGausAngleYvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
+ hSigGausAngleYvsP1->SetBinError(iHist+1,errorSigma);
+
+ limitGaus = nSigFit * (DeltaAngleX2[iHist]->GetRMS());
+ g5->SetRange(-limitGaus,limitGaus);
+ DeltaAngleX2[iHist]->Fit("g5","RQ");
+ sigmaAngleGaus = g5->GetParameter(2);
+ errorSigma = g5->GetParError(2);
+ printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
+ hSigGausAngleXvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
+ hSigGausAngleXvsP2->SetBinError(iHist+1,errorSigma);
+
+ limitGaus = nSigFit * (DeltaAngleY2[iHist]->GetRMS());
+ g6->SetRange(-limitGaus,limitGaus);
+ DeltaAngleY2[iHist]->Fit("g6","RQ");
+ sigmaAngleGaus = g6->GetParameter(2);
+ errorSigma = g6->GetParError(2);
+ printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
+ hSigGausAngleYvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
+ hSigGausAngleYvsP2->SetBinError(iHist+1,errorSigma);
+
+ }
// save histograms in MUONTestAbso.root
TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
+
hInvMassRes->Write();
+ hRap->Write();
+
hDeltaP1->Write();
hDeltaAngle1->Write();
+ hP1->Write();
+ hProfDeltaPvsP1->Write();
+ h2DeltaPvsP1->Write();
+ hSigmaPvsP1->Write();
+ hSigGausPvsP1->Write();
+ hSigGausAngleXvsP1->Write();
+ hSigGausAngleYvsP1->Write();
+
hDeltaP2->Write();
hDeltaAngle2->Write();
+ hP2->Write();
+ hProfDeltaPvsP2->Write();
+ h2DeltaPvsP2->Write();
+ hSigmaPvsP2->Write();
+ hSigGausPvsP2->Write();
+ hSigGausAngleXvsP2->Write();
+ hSigGausAngleYvsP2->Write();
+
+ for (Int_t iHist=0; iHist < nBinProf; iHist++){
+ DeltaP1[iHist]->Write();
+ DeltaP2[iHist]->Write();
+ DeltaAngleX1[iHist]->Write();
+ DeltaAngleY1[iHist]->Write();
+ DeltaAngleX2[iHist]->Write();
+ DeltaAngleY2[iHist]->Write();
+ }
+
+ hDeltaPxvsPhi->Write();
+ hDeltaPyvsPhi->Write();
+ hDeltaPhivsPz->Write();
+
histoFile->Close();
}
Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
}
+
+Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
+{
+ // return rapidity for particle
+
+ // Particle energy
+ Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
+ // Rapidity
+ Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
+ return Rapidity;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff