Modified absorber correction. Added function FieldCorrection() to account
authorcussonno <cussonno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Sep 2002 10:14:00 +0000 (10:14 +0000)
committercussonno <cussonno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Sep 2002 10:14:00 +0000 (10:14 +0000)
for the effect of magnetic field in absorber.

MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h
MUON/MUONTestAbso.C

index ddf51f3..4569ce9 100644 (file)
 
 /*
 $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
 
@@ -267,9 +271,11 @@ void AliMUONTrackParam::ExtrapToVertex()
   
   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);
 }
 
 
@@ -394,14 +400,15 @@ void AliMUONTrackParam::BransonCorrection()
   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;
   }
 
@@ -426,16 +433,18 @@ void AliMUONTrackParam::BransonCorrection()
   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 +
@@ -449,35 +458,77 @@ void AliMUONTrackParam::BransonCorrection()
   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);
+}
index 6d63aab..5ba7d9a 100644 (file)
@@ -38,7 +38,8 @@ class AliMUONTrackParam : public TObject {
   void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
   void ExtrapToVertex();  // extrapolation to vertex through the absorber
   void BransonCorrection(); // makes Branson correction
-  Double_t TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber); // returns total momentum after energy loss correction in the absorber
+  Double_t TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta); // returns total momentum after energy loss correction in the absorber
+  void FieldCorrection(Double_t Z); // makes simple magnetic field correction through the absorber 
 
  protected:
  private:
index 1390c23..1eecc34 100644 (file)
@@ -1,4 +1,4 @@
-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)
@@ -39,19 +39,85 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
     
   //  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)
@@ -59,9 +125,8 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
 
   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;
     }
     
@@ -72,8 +137,10 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
     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;
@@ -84,7 +151,8 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
        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;
        }
@@ -94,11 +162,12 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
        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");
@@ -106,7 +175,20 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
      
     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;
     }
@@ -121,13 +203,26 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=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();
@@ -138,7 +233,8 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
 
              // 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);
@@ -153,6 +249,7 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
                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;
@@ -171,7 +268,7 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
              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)
@@ -181,33 +278,215 @@ void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0)
            } // 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();
   
 }
@@ -221,3 +500,14 @@ Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Fl
   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;
+}