get tables from the aliroot directory if they are not in the current one
[u/mrichter/AliRoot.git] / 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 }