MUON/MUONTestAbso.C

   1 void MUONTestAbso (Int_t evNumber1=0,Int_t evNumber2=0,Int_t testSingle=0)
   2 {
   3   // Useful to test the absorber correction in the reconstruction procedure
   4   // (mean energy loss + Branson correction)
   5   // The AliMUONTrackParam class from the reconstruction is directly checked
   6   // in this macro using the GEANT particle momentum downstream of the absorber.
   7   // Histograms are saved in file MUONTestAbso.root.
   8   // Use DrawTestAbso.C to display control histograms.
   9
  10
  11   // Dynamically link some shared libs
  12   if (gClassTable->GetID("AliRun") < 0) {
  13     gROOT->LoadMacro("loadlibs.C");
  14     loadlibs();
  15   }
  16
  17   // Connect the Root Galice file containing Geometry, Kine and Hits
  18
  19   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  20
  21
  22   if (!file) {
  23     printf("\n Creating galice.root \n");
  24     file = new TFile("galice.root");
  25   } else {
  26     printf("\n galice.root found in file list");
  27   }
  28
  29 // Get AliRun object from file or create it if not on file
  30   if (!gAlice) {
  31     gAlice = (AliRun*)(file->Get("gAlice"));
  32     if (gAlice) printf("AliRun object found on file\n");
  33     if (!gAlice) {
  34       printf("\n Create new gAlice object");
  35       gAlice = new AliRun("gAlice","Alice test program");
  36     }
  37   }
  38
  39
  40   //  Create some histograms
  41
  42   TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 1200, 0., 12.);
  43 //   TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2)", 100, 9., 10.5);
  44   TH1F *hDeltaP1 = new TH1F("hDeltaP1", " delta P for muons theta < 3 deg", 100, -10., 10.);
  45   TH1F *hDeltaAngle1 = new TH1F("hDeltaAngle1", " delta angle for muons theta < 3 deg", 100, 0., 0.5);
  46   //
  47   TH1F *hDeltaP2 = new TH1F("hDeltaP2", " delta P for muons theta > 3 deg", 100, -10., 10.);
  48   TH1F *hDeltaAngle2 = new TH1F("hDeltaAngle2", " delta angle for muons theta > 3 deg", 100, 0., 0.5);
  49
  50   //
  51   TH1F *hRap = new TH1F("hRap"," Muon Rapidity gen.",20,2.5,4);
  52   //
  53
  54   Int_t nBinProf=19;
  55   Float_t xProfMin=5;
  56   Float_t xProfMax=290;
  57   char titleHist[50];
  58   char numberHist[10];
  59   TH1F *DeltaP1[100],*DeltaP2[100];
  60   TH1F *DeltaAngleX1[100],*DeltaAngleX2[100];
  61   TH1F *DeltaAngleY1[100],*DeltaAngleY2[100];
  62
  63   TH1F *hP1 = new TH1F("hP1"," P  theta < 3 deg ",nBinProf,xProfMin,xProfMax);
  64   TProfile *hProfDeltaPvsP1 = new TProfile("hProfDeltaPvsP1"," delta P  vs  P  theta < 3 deg ",nBinProf,xProfMin,xProfMax,-50,50,"s");
  65   TH2F *h2DeltaPvsP1 = new TH2F("h2DeltaPvsP1"," delta P  vs  P  theta < 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
  66    TH1F *hSigmaPvsP1 = new TH1F("hSigmaPvsP1"," deltaP/P vs P  theta < 3 deg",nBinProf,xProfMin,xProfMax);
  67    for (Int_t iHist=0; iHist < nBinProf; iHist++){
  68       sprintf(titleHist,"%s","hDelta1P");
  69       sprintf(numberHist,"%d",iHist+1);
  70       strcat(titleHist,numberHist);
  71       DeltaP1[iHist] = new TH1F(titleHist," deltaP  theta < 3 deg ",400,-50,50);
  72
  73       sprintf(titleHist,"%s","hDelta1AngleX");
  74       sprintf(numberHist,"%d",iHist+1);
  75       strcat(titleHist,numberHist);
  76       DeltaAngleX1[iHist] = new TH1F(titleHist," delta Angle X theta < 3 deg ",400,-0.05,0.05);
  77
  78       sprintf(titleHist,"%s","hDelta1AngleY");
  79       sprintf(numberHist,"%d",iHist+1);
  80       strcat(titleHist,numberHist);
  81       DeltaAngleY1[iHist] = new TH1F(titleHist," delta Angle Y theta < 3 deg ",400,-0.05,0.05);
  82    }
  83    TH1F *hSigGausPvsP1 = new TH1F("hSigGausPvsP1"," deltaP/P vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
  84    TH1F *hSigGausAngleXvsP1 = new TH1F("hSigGausAngleXvsP1"," delta theta X vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
  85    TH1F *hSigGausAngleYvsP1 = new TH1F("hSigGausAngleYvsP1"," delta theta Y vs P gauss theta < 3 deg",nBinProf,xProfMin,xProfMax);
  86
  87   TH1F *hP2 = new TH1F("hP2"," P  theta > 3 deg ",nBinProf,xProfMin,xProfMax);
  88   TProfile *hProfDeltaPvsP2 = new TProfile("hProfDeltaPvsP2"," delta P  vs  P  theta > 3 deg",nBinProf,xProfMin,xProfMax,-50,50,"s");
  89   TH2F *h2DeltaPvsP2 = new TH2F("h2DeltaPvsP2"," delta P  vs  P  theta > 3 deg",nBinProf,xProfMin,xProfMax,50,-50,50);
  90    TH1F *hSigmaPvsP2 = new TH1F("hSigmaPvsP2"," deltaP/P vs P theta > 3 deg",nBinProf,xProfMin,xProfMax);
  91    for (Int_t iHist=0; iHist < nBinProf; iHist++){
  92       sprintf(titleHist,"%s","hDelta2P");
  93       sprintf(numberHist,"%d",iHist+1);
  94       strcat(titleHist,numberHist);
  95       DeltaP2[iHist] = new TH1F(titleHist," deltaP  theta > 3 deg ",400,-50,50);
  96
  97       sprintf(titleHist,"%s","hDelta2AngleX");
  98       sprintf(numberHist,"%d",iHist+1);
  99       strcat(titleHist,numberHist);
 100       DeltaAngleX2[iHist] = new TH1F(titleHist," delta Angle X theta > 3 deg ",400,-0.05,0.05);
 101
 102       sprintf(titleHist,"%s","hDelta2AngleY");
 103       sprintf(numberHist,"%d",iHist+1);
 104       strcat(titleHist,numberHist);
 105       DeltaAngleY2[iHist] = new TH1F(titleHist," delta Angle Y theta > 3 deg ",400,-0.05,0.05);
 106    }
 107    TH1F *hSigGausPvsP2 = new TH1F("hSigGausPvsP2"," deltaP/P vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
 108
 109   TH1F *hSigGausAngleXvsP2 = new TH1F("hSigGausAngleXvsP2"," delta theta X vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
 110    TH1F *hSigGausAngleYvsP2 = new TH1F("hSigGausAngleYvsP2"," delta theta Y vs P gauss theta > 3 deg",nBinProf,xProfMin,xProfMax);
 111
 112    TProfile *hDeltaPxvsPhi = new TProfile("hDeltaPxvsPhi"," delta Px  vs  Phi",45,-180,180,-1,1,"s");
 113   TProfile *hDeltaPyvsPhi = new TProfile("hDeltaPyvsPhi"," delta Py  vs  Phi",45,-180,180,-1,1,"s");
 114   TProfile *hDeltaPhivsPz = new TProfile("hDeltaPhivsPz"," delta phi  vs pZ",25,0,100,-4,4);
 115
 116   // Define constant parameters
 117   Float_t massMuon = 0.105658389;  // muon mass
 118   Float_t massUpsilon = 9.46037;   // Upsilon mass
 119   Double_t zEndAbso = 503;         // z position of the end of the absorber
 120   Double_t rLimit = zEndAbso * TMath::Tan(3*TMath::Pi()/180); // 3 deg. limit (different absorber composition)
 121   Double_t printLevel = 0;
 122
 123   if (printLevel >= 1)
 124   cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;
 125
 126   Double_t pX[2],pY[2],pZ[2],pTot[2],theta[2],radius[2]; // extrapolated parameters of particles from chamber 1 to vertex
 127
 128   for (Int_t nev=evNumber1; nev<= evNumber2; nev++) {  // loop over events
 129     if (printLevel >= 0 && nev%100 == 0) {
 130       cout << " **** Event # " << nev <<endl;
 131     }
 132
 133     Int_t nparticles = gAlice->GetEvent(nev);
 134     if (printLevel >= 1) {
 135       cout << " Total number of nparticles  " << nparticles <<endl;
 136     }
 137     if (nev < evNumber1) continue;
 138     if (nparticles <= 0) return;
 139
 140     Double_t PxGen[2],PyGen[2],PzGen[2],PTotGen[2],ThetaGen[2],RapGen[2]; // parameters of generated  particles at the vertex
 141     ThetaGen[0] = ThetaGen[1] = 0;
 142     RapGen[0] = RapGen[1] = 0;
 143
 144     for (int iPart = 0; iPart < nparticles  ; iPart++) {  // loop over particles
 145
 146       TParticle *particle;
 147       particle = gAlice->Particle(iPart);
 148
 149       if ( particle->GetPdgCode() == kMuonPlus ) {
 150         PxGen[0] = particle->Px();
 151         PyGen[0] = particle->Py();
 152         PzGen[0] = particle->Pz();
 153         PTotGen[0] = TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]+PzGen[0]*PzGen[0]);
 154         ThetaGen[0] = TMath::ATan2(TMath::Sqrt(PxGen[0]*PxGen[0]+PyGen[0]*PyGen[0]),PzGen[0]);
 155         RapGen[0] = rapParticle(PxGen[0],PyGen[0],PzGen[0],massMuon);
 156         if (printLevel >= 1) {
 157           cout << " Particle id : "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[0]<<" "<<PyGen[0]<<" "<<PzGen[0]<<" "<<ThetaGen[0]<<endl;
 158         }
 159       }
 160       if ( particle->GetPdgCode() == kMuonMinus ) {
 161         PxGen[1] = particle->Px();
 162         PyGen[1] = particle->Py();
 163         PzGen[1] = particle->Pz();
 164         PTotGen[1] = TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]+PzGen[1]*PzGen[1]);
 165         ThetaGen[1] = TMath::ATan2(TMath::Sqrt(PxGen[1]*PxGen[1]+PyGen[1]*PyGen[1]),PzGen[1]);
 166         RapGen[1] = rapParticle(PxGen[1],PyGen[1],PzGen[1],massMuon);
 167         if (printLevel >= 1)
 168           cout << " Particle #: "<<particle->GetPdgCode()<<" px,py,pz,theta : "<<PxGen[1]<<" "<<PyGen[1]<<" "<<PzGen[1]<<" "<<ThetaGen[1]<<endl;
 169       }
 170     } // end loop over particles
 171
 172
 173     AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
 174
 175
 176     TTree *TH = gAlice->TreeH();
 177     Int_t ntracks =(Int_t) TH->GetEntries();
 178
 179     Bool_t hitChamberMUPlus[10],hitChamberMUMinus[10];
 180     Bool_t hitStationMUPlus[5],hitStationMUMinus[5];
 181     for (Int_t iCh=0; iCh < 10; iCh++){
 182        hitChamberMUPlus[iCh] = kFALSE;
 183        hitChamberMUMinus[iCh] = kFALSE;
 184     }
 185     for (Int_t iSt=0; iSt < 5; iSt++){
 186        hitStationMUPlus[iSt] = kFALSE;
 187        hitStationMUMinus[iSt] = kFALSE;
 188     }
 189
 190     Int_t nHitChamber=0;
 191
 192     for (Int_t i = 0; i < 2; i++){
 193       pX[i] = pY[i] = pZ[i] = pTot[i] = theta[i] = radius[i] = 0;
 194     }
 195
 196     for (Int_t track=0; track<ntracks;track++) { // loop over tracks
 197       gAlice->ResetHits();
 198       Int_t nbytes += TH->GetEvent(track);
 199       if (MUON)  {
 200         for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 201             mHit;
 202             mHit=(AliMUONHit*)MUON->NextHit())  // loop over GEANT muon hits
 203           {
 204             Int_t   nch   = mHit->Chamber();  // chamber number
 205
 206
 207             Int_t ipart = mHit->Particle();
 208             Double_t x     = mHit->X();        // x-pos of hit in chamber #1
 209             Double_t y     = mHit->Y();        // y-pos of hit in chamber #1
 210             Double_t z     = mHit->Z();        // z-pos of hit in chamber #1
 211
 212             Int_t iSt = (nch-1)/2;
 213             if ((nch-1) < 10) {
 214               if (ipart == kMuonPlus){
 215                 hitChamberMUPlus[nch-1] = kTRUE;
 216                 hitStationMUPlus[iSt] = kTRUE;
 217               }
 218               if (ipart == kMuonMinus){
 219                 hitChamberMUMinus[nch-1] = kTRUE;
 220                 hitStationMUMinus[iSt] = kTRUE;
 221               }
 222               nHitChamber++;
 223             }
 224
 225             if (nch != 1) continue;
 226
 227             if (ipart == kMuonPlus || ipart == kMuonMinus) {
 228               Double_t px0=mHit->Px();
 229               Double_t py0=mHit->Py();
 230               Double_t pz0=mHit->Pz();
 231               Double_t nonBendingSlope=px0/pz0;
 232               Double_t bendingSlope=py0/pz0;
 233
 234               // create an AliMUONTrackParam object in chamber #1
 235               AliMUONTrackParam *trackParam = new AliMUONTrackParam();
 236               Int_t sign = -TMath::Sign(1.0,ipart);
 237               Double_t bendingMometum = sign* TMath::Sqrt(pz0*pz0+py0*py0); // Bug px -> py !!!
 238               Double_t inverseBendingMomentum = 1/bendingMometum;
 239
 240               trackParam->SetBendingCoor(y);
 241               trackParam->SetNonBendingCoor(x);
 242               trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
 243               trackParam->SetBendingSlope(bendingSlope);
 244               trackParam->SetNonBendingSlope(nonBendingSlope);
 245               trackParam->SetZ(z);
 246
 247               if (printLevel >= 1) {
 248                 cout <<" **** before extrap " <<endl;
 249                 cout << " px, py, pz = " <<px0<<" "<<py0<<" "<<pz0<<endl;
 250                 trackParam->Dump();
 251               }
 252               if (!testSingle)
 253               trackParam->ExtrapToVertex(); // extrapolate track parameters through the absorber
 254               if (printLevel >= 1) {
 255                 cout <<" **** after extrap " <<endl;
 256                 trackParam->Dump();
 257               }
 258
 259               Int_t iMuon;
 260               if (ipart == kMuonPlus) iMuon = 0;
 261               else iMuon = 1;
 262
 263               // calculate track parameters extrapolated to vertex.
 264               bendingSlope = trackParam->GetBendingSlope();
 265               nonBendingSlope = trackParam->GetNonBendingSlope();
 266               Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
 267               pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
 268               pX[iMuon] = pZ[iMuon] * nonBendingSlope;
 269               pY[iMuon] = pZ[iMuon] * bendingSlope;
 270               pTot[iMuon] = TMath::Sqrt(pYZ*pYZ+pX[iMuon]*pX[iMuon]);
 271               theta[iMuon] = TMath::ATan2(TMath::Sqrt(pX[iMuon]*pX[iMuon]+pY[iMuon]*pY[iMuon]),pZ[iMuon]);
 272               radius[iMuon] = TMath::Sqrt(x*x+y*y);
 273
 274               if (printLevel >= 1)
 275                 cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
 276               delete trackParam;
 277
 278             } // if MuonPlus or MuonMinus
 279           } // loop over MUON hits
 280       } // if MUON
 281     } // loop over tracks
 282
 283
 284     Bool_t goodMUPlus = kTRUE;
 285     Bool_t goodMUMinus = kTRUE;
 286     for (Int_t iCh=0; iCh < 1; iCh++) { // !!!!  10 -> 1
 287       if (!hitChamberMUPlus[iCh]) {
 288         goodMUPlus = kFALSE;
 289         printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
 290       }
 291       if (!hitChamberMUMinus[iCh]) {
 292         goodMUMinus = kFALSE;
 293         printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
 294       }
 295     }
 296
 297 //     for (Int_t iSt=0; iSt < 5; iSt++) {
 298 //       if (!hitStationMUPlus[iSt]) {
 299 //      goodMUPlus = kFALSE;
 300 // //   printf(" evt number %d no muon+ in chamber %d \n",nev,iCh+1);
 301 //       }
 302 //       if (!hitStationMUMinus[iSt]) {
 303 //      goodMUMinus = kFALSE;
 304 // //   printf(" evt number %d no muon- in chamber %d \n",nev,iCh+1);
 305 //       }
 306 //     }
 307
 308     if (pX[0] != 0 && pX[1] != 0) { // if track 1 & 2 exist
 309       Float_t massResonance = MuPlusMuMinusMass(pX[0],pY[0],pZ[0],pX[1],pY[1],pZ[1],massMuon);
 310       hInvMassRes->Fill(massResonance);
 311     }
 312
 313
 314     for (Int_t i = 0; i<2; i++){  // !!! 1 -> 2
 315       if (i == 0 && !goodMUPlus) continue;
 316       if (i == 1 && !goodMUMinus) continue;
 317       if (pX[i] == 0) continue;
 318
 319       hRap->Fill(RapGen[i]);
 320
 321       Double_t cosAngle = pX[i]*PxGen[i]+pY[i]*PyGen[i]+pZ[i]*PzGen[i];
 322       cosAngle = cosAngle/(pTot[i]*PTotGen[i]);
 323       Double_t deltaAngle = TMath::ACos(cosAngle)*180/TMath::Pi();
 324       Double_t deltaPNorm = 0;
 325       if (testSingle) {
 326         deltaPNorm = (PTotGen[i]-pTot[i])* TMath::Cos(ThetaGen[i]);
 327       } else {
 328         deltaPNorm = PTotGen[i]-pTot[i];
 329       }
 330
 331       Float_t thetaGX = TMath::ATan2(PxGen[i],PzGen[i]);
 332       Float_t thetaRX = TMath::ATan2(pX[i],pZ[i]);
 333       Float_t deltaAngleX = thetaRX - thetaGX;
 334
 335       Float_t thetaGY = TMath::ATan2(PyGen[i],PzGen[i]);
 336       Float_t thetaRY = TMath::ATan2(pY[i],pZ[i]);
 337       Float_t deltaAngleY = thetaRY - thetaGY;
 338
 339       Float_t phiG = TMath::ATan2(PyGen[i],PxGen[i])*180/TMath::Pi();
 340       Float_t phiR = TMath::ATan2(pY[i],pX[i])*180/TMath::Pi();
 341
 342       hDeltaPxvsPhi->Fill(phiG,PxGen[i]-pX[i]);
 343       hDeltaPyvsPhi->Fill(phiG,PyGen[i]-pY[i]);
 344       hDeltaPhivsPz->Fill(PzGen[i],phiR-phiG);
 345
 346       Int_t iHist=0;
 347       if (ThetaGen[i] < (3*TMath::Pi()/180)) {
 348         hDeltaP1->Fill(pTot[i]-PTotGen[i]);
 349         hDeltaAngle1->Fill(deltaAngle);
 350         hP1->Fill(PTotGen[i]);
 351         hProfDeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
 352         h2DeltaPvsP1->Fill(PTotGen[i],deltaPNorm);
 353         iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
 354         if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax  ) {
 355           DeltaP1[iHist]->Fill(deltaPNorm);
 356           DeltaAngleX1[iHist]->Fill(deltaAngleX);
 357           DeltaAngleY1[iHist]->Fill(deltaAngleY);
 358         }
 359       } else {
 360         hDeltaP2->Fill(pTot[i]-PTotGen[i]);
 361         hDeltaAngle2->Fill(deltaAngle);
 362         hP2->Fill(PTotGen[i]);
 363         hProfDeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
 364         h2DeltaPvsP2->Fill(PTotGen[i],deltaPNorm);
 365         iHist = (PTotGen[i]-xProfMin)*nBinProf/(xProfMax-xProfMin);
 366         if (PTotGen[i] > xProfMin && PTotGen[i] < xProfMax ) {
 367           DeltaP2[iHist]->Fill(deltaPNorm);
 368           DeltaAngleX2[iHist]->Fill(deltaAngleX);
 369           DeltaAngleY2[iHist]->Fill(deltaAngleY);
 370         }
 371       }
 372     }
 373   } // loop over event
 374
 375
 376   Float_t sigmaP = 0;
 377   Float_t xBin = 0;
 378   for (Int_t iBin = 1; iBin <= nBinProf; iBin++){
 379     sigmaP = hProfDeltaPvsP1->GetBinError(iBin);
 380     xBin = hProfDeltaPvsP1->GetBinCenter(iBin);
 381     hSigmaPvsP1->SetBinContent(iBin,sigmaP/xBin);
 382     sigmaP = hProfDeltaPvsP2->GetBinError(iBin);
 383     xBin = hProfDeltaPvsP2->GetBinCenter(iBin);
 384     hSigmaPvsP2->SetBinContent(iBin,sigmaP/xBin);
 385   }
 386
 387   TF1  *g1= new TF1("g1","gaus",-25,25) ;
 388   TF1  *g2= new TF1("g2","gaus",-25,25) ;
 389   Float_t sigmaPGaus;
 390   Float_t xBinGaus;
 391   for (Int_t iHist=0; iHist < nBinProf; iHist++){
 392     DeltaP1[iHist]->Fit("g1","RQ");
 393     sigmaPGaus =  g1->GetParameter(2);
 394     printf(" ** 1 ** iHist= %d , sigmaPGaus = %f \n",iHist,sigmaPGaus);
 395     xBinGaus = hSigGausPvsP1->GetBinCenter(iHist+1);
 396     hSigGausPvsP1->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
 397
 398     DeltaP2[iHist]->Fit("g2","RQ");
 399     sigmaPGaus =  g2->GetParameter(2);
 400     printf(" ** 2 ** iHist= %d , sigmaGaus = %f \n",iHist,sigmaPGaus);
 401     xBinGaus = hSigGausPvsP2->GetBinCenter(iHist+1);
 402     hSigGausPvsP2->SetBinContent(iHist+1,sigmaPGaus/xBinGaus);
 403   }
 404
 405
 406   // Angles
 407   TF1  *g3= new TF1("g3","gaus") ;
 408   TF1  *g4= new TF1("g4","gaus") ;
 409   TF1  *g5= new TF1("g5","gaus") ;
 410   TF1  *g6= new TF1("g6","gaus") ;
 411   Float_t sigmaAngleGaus,limitGaus,errorSigma;
 412   Float_t nSigFit = 3;
 413   for (Int_t iHist=0; iHist < nBinProf; iHist++){
 414     limitGaus = nSigFit * (DeltaAngleX1[iHist]->GetRMS());
 415     g3->SetRange(-limitGaus,limitGaus);
 416     DeltaAngleX1[iHist]->Fit("g3","RQ");
 417     sigmaAngleGaus =  g3->GetParameter(2);
 418     errorSigma = g3->GetParError(2);
 419     printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
 420     hSigGausAngleXvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
 421     hSigGausAngleXvsP1->SetBinError(iHist+1,errorSigma);
 422
 423     limitGaus = nSigFit * (DeltaAngleY1[iHist]->GetRMS());
 424     g4->SetRange(-limitGaus,limitGaus);
 425     DeltaAngleY1[iHist]->Fit("g4","RQ");
 426     sigmaAngleGaus =  g4->GetParameter(2);
 427     errorSigma = g4->GetParError(2);
 428     printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
 429     hSigGausAngleYvsP1->SetBinContent(iHist+1,sigmaAngleGaus);
 430     hSigGausAngleYvsP1->SetBinError(iHist+1,errorSigma);
 431
 432     limitGaus = nSigFit * (DeltaAngleX2[iHist]->GetRMS());
 433     g5->SetRange(-limitGaus,limitGaus);
 434     DeltaAngleX2[iHist]->Fit("g5","RQ");
 435     sigmaAngleGaus =  g5->GetParameter(2);
 436     errorSigma = g5->GetParError(2);
 437     printf(" ** 1 ** iHist= %d , sigmaAngleGaus X = %f \n",iHist,sigmaAngleGaus);
 438     hSigGausAngleXvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
 439     hSigGausAngleXvsP2->SetBinError(iHist+1,errorSigma);
 440
 441     limitGaus = nSigFit * (DeltaAngleY2[iHist]->GetRMS());
 442     g6->SetRange(-limitGaus,limitGaus);
 443     DeltaAngleY2[iHist]->Fit("g6","RQ");
 444     sigmaAngleGaus =  g6->GetParameter(2);
 445     errorSigma = g6->GetParError(2);
 446     printf(" ** 1 ** iHist= %d , sigmaAngleGaus Y = %f \n",iHist,sigmaAngleGaus);
 447     hSigGausAngleYvsP2->SetBinContent(iHist+1,sigmaAngleGaus);
 448     hSigGausAngleYvsP2->SetBinError(iHist+1,errorSigma);
 449
 450   }
 451   // save histograms in MUONTestAbso.root
 452   TFile *histoFile = new TFile("MUONTestAbso.root", "RECREATE");
 453
 454   hInvMassRes->Write();
 455   hRap->Write();
 456
 457   hDeltaP1->Write();
 458   hDeltaAngle1->Write();
 459   hP1->Write();
 460   hProfDeltaPvsP1->Write();
 461   h2DeltaPvsP1->Write();
 462   hSigmaPvsP1->Write();
 463   hSigGausPvsP1->Write();
 464   hSigGausAngleXvsP1->Write();
 465   hSigGausAngleYvsP1->Write();
 466
 467   hDeltaP2->Write();
 468   hDeltaAngle2->Write();
 469   hP2->Write();
 470   hProfDeltaPvsP2->Write();
 471   h2DeltaPvsP2->Write();
 472   hSigmaPvsP2->Write();
 473   hSigGausPvsP2->Write();
 474   hSigGausAngleXvsP2->Write();
 475   hSigGausAngleYvsP2->Write();
 476
 477   for (Int_t iHist=0; iHist < nBinProf; iHist++){
 478      DeltaP1[iHist]->Write();
 479      DeltaP2[iHist]->Write();
 480      DeltaAngleX1[iHist]->Write();
 481      DeltaAngleY1[iHist]->Write();
 482      DeltaAngleX2[iHist]->Write();
 483      DeltaAngleY2[iHist]->Write();
 484   }
 485
 486   hDeltaPxvsPhi->Write();
 487   hDeltaPyvsPhi->Write();
 488   hDeltaPhivsPz->Write();
 489
 490   histoFile->Close();
 491
 492 }
 493
 494
 495 Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
 496 {
 497   // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass
 498
 499   Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
 500   Float_t e2 = TMath::Sqrt(MuonMass * MuonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
 501   return (TMath::Sqrt(2.0 * (MuonMass * MuonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
 502 }
 503
 504 Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
 505 {
 506   // return rapidity for particle
 507
 508   // Particle energy
 509   Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
 510   // Rapidity
 511   Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
 512   return Rapidity;
 513 }