+++ /dev/null
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include "iostream.h"
-#include "/home/alice/guernane/date/tmp/ionlx/event.h"
-#include "/home/alice/guernane/date/tmp/ionlx/monitor.h"
-#endif
-
-//=======================================================================
-//
-// The following macro is used to read raw data from run#.raw files
-// generated by DATE acquisition package. It was written using
-// Guy Jacquet (IPN Lyon) monitoring code.
-//
-// 1-Feb-2000 Rachid GUERNANE, IPN Lyon, France
-//
-//
-
-void MUONRawDigit (Int_t evNumber1 = 0, Int_t evNumber2 = 0, char* dataSource="/home/alice/guernane/test/run3403_Linux.raw")
-{
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile* file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (file) file->Close();
- file = new TFile("galice.root", "UPDATE");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-
- AliMUONv2* MUON = (AliMUONv2*)gAlice->GetModule("MUON");
-
-//
-// Start Event Loop
-//
-
- static int ilen;
- int status;
-
-
- status = monitorSetDataSource(dataSource);
-
- if (status != 0) {
- fprintf( stderr,
- "Error in monitorSetDataSource: %s\n",
- monitorSetDataSource(status));
- }
-
- status = monitorDeclareMp("Retrieve Raw Data");
-
- if (status != 0) {
- fprintf( stderr,
- "Error in monitorDeclareMp: %s\n",
- monitorSetDataSource(status));
- }
-
- Int_t nev = 0;
- Int_t iev = 0;
-
- for (; iev <= evNumber2 && nev < evNumber2*20; nev++) {
-
-//
-// Raw event decoding
-//
-
- void *ptr;
- struct eventStruct *rawevent;
-
- status = monitorGetEventDynamic(&ptr);
-
- if (status != 0) {
- fprintf( stderr,
- "Error in monitorGetEventDynamic: %s\n",
- monitorDecodeError(status));
- exit(1);
- }
-
- rawevent = (struct eventStruct*)ptr;
-
- ilen = (rawevent->eventHeader.size - rawevent->eventHeader.headLen)/4;
-
- if (rawevent->eventHeader.type == PHYSICS_EVENT) {
- iev = rawevent->eventHeader.nbInRun;
- gAlice->GetEvent(iev-1);
- printf("\nEvent No.: %d\n", iev);
- printf("--------------------------------------------\n");
-
- int* lptr = (int*)&rawevent->rawData[0];
-
- MUON->GetRawDigits(iev-1, lptr, ilen);
-
- }
-
- free(ptr);
-
- }
-//
-// End Event loop
-//
-
- file->Close();
-
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-/* $Id$ */
-
-// Macro MUONTracker.C (TO BE COMPILED)
-// for testing the C++ reconstruction code
-// Output is using aliroot standard output MUON.Tracks.root
-// The output is a TClonesArray of AliMUONTracks.
-// Gines MARTINEZ, Subatech, sep 2003
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <TClonesArray.h>
-
-#include "AliRun.h"
-#include "AliESD.h"
-#include "AliMUON.h"
-#include "AliMUONData.h"
-#include "AliMUONEventReconstructor.h"
-#include "AliMUONTrack.h"
-#include "AliMUONTrackHit.h"
-#include "AliMUONTrackParam.h"
-#endif
-
-void MUONRecoTest(Text_t *FileName = "galice.root")
-{
- //
- cout << "MUONRecoTest" << endl;
- cout << "FileName ``" << FileName << "''" << endl;
-
- // Creating Run Loader and openning file containing Hits, Digits and RecPoints
- AliRunLoader * RunLoader = AliRunLoader::Open(FileName,"Event","UPDATE");
- if (RunLoader ==0x0) {
- printf(">>> Error : Error Opening %s file \n",FileName);
- return;
- }
- // Loading AliRun master
- RunLoader->LoadgAlice();
- gAlice = RunLoader->GetAliRun();
-
- // Loading MUON subsystem
- AliMUON* MUON = (AliMUON *) gAlice->GetDetector("MUON");
-
- //AliESD *event = new AliESD();
- MUON->Reconstruct();
- // MUON->FillESD(event);
-}
+++ /dev/null
-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)
- // The AliMUONTrackParam class from the reconstruction is directly checked
- // in this macro using the GEANT particle momentum downstream of the absorber.
- // Histograms are saved in file MUONTestAbso.root.
- // Use DrawTestAbso.C to display control histograms.
-
-
- // Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
- // Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n Create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-
- // Create some histograms
-
- 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 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 zEndAbso = 503; // z position of the end of the absorber
- Double_t rLimit = zEndAbso * TMath::Tan(3*TMath::Pi()/180); // 3 deg. limit (different absorber composition)
- Double_t printLevel = 0;
-
- if (printLevel >= 1)
- cout <<" **** z end Absorber : "<< zEndAbso <<" rLimit : "<< rLimit << endl;
-
- 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=evNumber1; nev<= evNumber2; nev++) { // loop over events
- if (printLevel >= 0 && nev%100 == 0) {
- cout << " **** Event # " << nev <<endl;
- }
-
- Int_t nparticles = gAlice->GetEvent(nev);
- if (printLevel >= 1) {
- cout << " Total number of nparticles " << nparticles <<endl;
- }
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- 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;
- particle = gAlice->Particle(iPart);
-
- if ( particle->GetPdgCode() == kMuonPlus ) {
- PxGen[0] = particle->Px();
- 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]);
- 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;
- }
- }
- if ( particle->GetPdgCode() == kMuonMinus ) {
- PxGen[1] = particle->Px();
- 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]);
- 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;
- }
-
- for (Int_t track=0; track<ntracks;track++) { // loop over tracks
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (MUON) {
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit()) // loop over GEANT muon hits
- {
- Int_t nch = mHit->Chamber(); // chamber number
-
-
- 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();
- Double_t py0=mHit->Py();
- Double_t pz0=mHit->Pz();
- Double_t nonBendingSlope=px0/pz0;
- Double_t bendingSlope=py0/pz0;
-
- // create an AliMUONTrackParam object in chamber #1
- AliMUONTrackParam *trackParam = new AliMUONTrackParam();
- 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);
- trackParam->SetNonBendingCoor(x);
- trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
- trackParam->SetBendingSlope(bendingSlope);
- trackParam->SetNonBendingSlope(nonBendingSlope);
- trackParam->SetZ(z);
-
- if (printLevel >= 1) {
- cout <<" **** before extrap " <<endl;
- 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;
- trackParam->Dump();
- }
-
- Int_t iMuon;
- if (ipart == kMuonPlus) iMuon = 0;
- else iMuon = 1;
-
- // calculate track parameters extrapolated to vertex.
- bendingSlope = trackParam->GetBendingSlope();
- nonBendingSlope = trackParam->GetNonBendingSlope();
- Double_t pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
- pZ[iMuon] = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
- 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]);
- radius[iMuon] = TMath::Sqrt(x*x+y*y);
-
- if (printLevel >= 1)
- cout << " px, py, pz, theta, radius = " <<pX[iMuon]<<" "<<pY[iMuon]<<" "<<pZ[iMuon]<<" "<<theta[iMuon]<<" "<<radius[iMuon]<<endl;
- delete trackParam;
-
- } // if MuonPlus or MuonMinus
- } // loop over MUON hits
- } // if MUON
- } // 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++){ // !!! 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 MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t MuonMass)
-{
- // return invariant mass for particle 1 & 2 whose masses are equal to MuonMass
-
- Float_t e1 = TMath::Sqrt(MuonMass * MuonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
- 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;
-}
+++ /dev/null
-#include "iostream.h"
-// example of macro which retrieves the Trigger Output from galice.root
-
-void MUONTestTrigger (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-// Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (file) file->Close();
- file = new TFile("galice.root","READ");
-
-// Get AliRun object from file or create it if not on file
- printf ("I'm after Map \n");
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
- printf ("I'm after gAlice \n");
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
- for (int nev=evNumber1; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev npart =" << nev << " , " << nparticles << "\n";
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- Int_t nbytes = 0;
- AliMUONGlobalTrigger *gloTrg;
- AliMUONLocalTrigger *locTrg;
-
-// Get pointers to Alice detectors and Triggers containers
- AliMUON *MUON = (AliMUON*)gAlice->GetModule("MUON");
- TClonesArray *globalTrigger = MUON->GlobalTrigger();
- TClonesArray *localTrigger = MUON->LocalTrigger();
-
- TTree *TR = gAlice->TreeR();
- TBranch *gbranch = TR->GetBranch("MUONGlobalTrigger");
- TBranch *lbranch = TR->GetBranch("MUONLocalTrigger");
- gbranch->SetAddress(&globalTrigger);
- lbranch->SetAddress(&localTrigger);
-
- Int_t nent = TR->GetEntries();
- cout << ">>> " << nent << " entries found in TreeR of event "
- << nev << "\n";
- Int_t nb = 0;
-
- for (Int_t n=1; n<nent; n++) {
- MUON->ResetTrigger();
- nbytes += TR->GetEvent(n);
-
- Int_t nglobals = globalTrigger->GetEntries(); // should be = 1
- Int_t nlocals = localTrigger->GetEntries(); // who knows ?
-
- for (Int_t i=0; i<nglobals; i++) { // inspect Global Trigger
- cout << " >>> Output for Global Trigger " << "\n";
-
- gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(i);
-
- cout << "fSinglePlusLpt = " << gloTrg->SinglePlusLpt() << "\n";
- cout << "fSinglePlusHpt = " << gloTrg->SinglePlusHpt() << "\n";
- cout << "fSinglePlusApt = " << gloTrg->SinglePlusApt() << "\n";
-
- cout << "fSingleMinusLpt = " << gloTrg->SingleMinusLpt() << "\n";
- cout << "fSingleMinusHpt = " << gloTrg->SingleMinusHpt() << "\n";
- cout << "fSingleMinusApt = " << gloTrg->SingleMinusApt() << "\n";
-
- cout << "fSingleUndefLpt = " << gloTrg->SingleUndefLpt() << "\n";
- cout << "fSingleUndefHpt = " << gloTrg->SingleUndefHpt() << "\n";
- cout << "fSingleUndefApt = " << gloTrg->SingleUndefApt() << "\n";
-
- cout << "fPairLikeLpt = " << gloTrg->PairLikeLpt() << "\n";
- cout << "fPairLikeHpt = " << gloTrg->PairLikeHpt() << "\n";
- cout << "fPairLikeApt = " << gloTrg->PairLikeApt() << "\n";
-
- cout << "fPairUnlikeLpt = " << gloTrg->PairUnlikeLpt() << "\n";
- cout << "fPairUnlikeHpt = " << gloTrg->PairUnlikeHpt() << "\n";
- cout << "fPairUnlikeApt = " << gloTrg->PairUnlikeApt() << "\n";
- } // end of loop on Global Trigger
-
- for (Int_t i=0; i<nlocals; i++) { // inspect Local Trigger
- cout << " >>> Output for Local Trigger # " << i << "\n";
-
- locTrg = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(i);
-
- cout << "fLoCircuit = " << locTrg->LoCircuit() << "\n";
- cout << "fLoStripX = " << locTrg->LoStripX() << "\n";
- cout << "fLoDev = " << locTrg->LoDev() << "\n";
- cout << "fLoStripY = " << locTrg->LoStripY() << "\n";
- cout << "fLoLpt = " << locTrg->LoLpt() << "\n";
- cout << "fLoHpt = " << locTrg->LoHpt() << "\n";
- cout << "fLoApt = " << locTrg->LoApt() << "\n";
-
- } // end of loop on Local Trigger
- } // end of loop on entries of TreeR
- } // loop on event
-
-}
-
-
-
-
+++ /dev/null
-void MUONcombi (Int_t evNumber=0)
-{
- const Float_t runWeight = 4.e8;
-
-// Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- } else {
- delete gAlice;
- gAlice = 0;
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root");
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-
-// Import the Kine and Hits Trees for the event evNumber in the file
- Int_t nparticles = gAlice->GetEvent(evNumber);
- if (nparticles <= 0) return;
-//
-//
- TH1F *dmass = new TH1F("dmass","Dimuon-Mass Distribution"
- ,100,0.,5.);
-
- TH1F *dmassc = new TH1F("dmassc","Dimuon-Mass Distribution"
- ,200,0.,10.);
- TH1F *dmassd = new TH1F("dmassd","Dimuon-Mass Distribution"
- ,200,0.,10.);
-
- TH1F *pt = new TH1F("pt","pT-single"
- ,50,0.,10.);
- TH1F *hCont[30];
-
-//
-// Generator Loop
-//
-
- Int_t i=0;
-
- AliGenCocktailEntry *Entry, *e1, *e2;
- AliGenCocktail* Cocktail = (AliGenCocktail*) gAlice->Generator();
-// Single Generator
- for (Entry=Cocktail->FirstGenerator();
- Entry;
- Entry=Cocktail->NextGenerator()
- ) {
- Entry->PrintInfo();
- }
-//
-// Initialize Combinator
-//
- AliDimuCombinator* Combinator = new AliDimuCombinator();
- Combinator->SetEtaCut(2.5, 4.);
- Combinator->SetPtMin(1.0);
-
- i=0;
-
-//
-// Single Muon Loop
-//
- Combinator->ResetRange();
-
-
-
- for(Muon=Combinator->FirstMuon();
- Muon;
- Muon=Combinator->NextMuon()) {
-//
- Int_t chfirst = Muon->GetFirstDaughter();
- Int_t chlast = Muon->GetLastDaughter();
- Int_t parent = Muon->GetFirstMother();
- Float_t ptm = Muon->Pt();
- Float_t eta = Muon->Eta();
-
-// printf("\n Particle %d Parent %d first child %d last child %d",
-// i,parent, chfirst, chlast);
-// printf("\n Particle pt, eta: %f , %f ", ptm, eta);
- i++;
- }
-//
-// Di-Muon Loop
-
- Float_t pt1,pt2;
- TParticle* Muon1;
- TParticle* Muon2;
-
- Combinator->ResetRange();
-
-//
-// Dimuon Loop controlled by Generator Loop
-//
- Float_t sig = 0;
- Int_t icont = 0;
- char name[30];
- for (Cocktail->FirstGeneratorPair(e1,e2);
- (e1&&e2);
- Cocktail->NextGeneratorPair(e1,e2)
- ){
-
- sprintf(name, "%s-%s", e1->GetName(), e2->GetName());
- hCont[icont] = new TH1F(name,"Dimuon-Mass Distribution",100,0.,5.);
- printf("\n ----------------------------------------------------");
- e1->PrintInfo();
- e2->PrintInfo();
- printf("\n ----------------------------------------------------");
- Combinator->SetFirstRange (e1->GetFirst(), e1->GetLast());
- Combinator->SetSecondRange(e2->GetFirst(), e2->GetLast());
- Combinator->SetRate(e1->Rate(), e2->Rate());
- for (Combinator->FirstMuonPairSelected(Muon1,Muon2);
- (Muon1 && Muon2);
- Combinator->NextMuonPairSelected(Muon1,Muon2))
- {
- pt1 = Muon1->Pt();
- pt2 = Muon2->Pt();
- Float_t mass = Combinator->Mass(Muon1, Muon2);
- Float_t wgt = runWeight*Combinator->Weight(Muon1, Muon2)*
- Combinator->DecayProbability(Muon1)*
- Combinator->DecayProbability(Muon2);
- pt->Fill(pt1, wgt);
- pt->Fill(pt2, wgt);
- Float_t smeared_mass = mass;
- Combinator->SmearGauss(0.02*mass, smeared_mass);
-
- if (TMath::Min(pt1,pt2) > -0.5*TMath::Max(pt1,pt2)+2.) {
- if (Combinator->Correlated(Muon1, Muon2)) {
- dmassc->Fill(smeared_mass, wgt);
- sig += wgt;
- } else {
-// account for the fact that we sum like-sign and unlike-sign
- wgt *= 0.5;
- }
- dmass->Fill(smeared_mass, wgt);
- hCont[icont]->Fill(smeared_mass, wgt);
- }
- } // Dimuon Loop
- icont++;
- }// Generator Loop
- printf("\n Signal %e \n \n", sig);
-
-//
-//Create a canvas, set the view range, show histograms
-//
- TCanvas *c1 = new TCanvas("c1","Dimuon Plots",400,10,600,700);
- TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(11);
- pad12->SetFillColor(11);
- pad13->SetFillColor(11);
- pad14->SetFillColor(11);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- pt->SetFillColor(42);
- pt->SetXTitle("pT (GeV)");
- pt->Draw();
-
- pad12->cd();
- pad12->SetLogy(0);
- dmass->SetFillColor(42);
- dmass->SetXTitle("m (GeV)");
- dmass->Draw();
-
- pad13->cd();
- pad13->SetLogy(0);
- dmassc->SetFillColor(42);
- dmassc->SetXTitle("m (GeV)");
- dmassc->Draw();
-
- pad14->cd();
- pad14->SetLogy(0);
- dmassd->SetFillColor(42);
- dmassd->SetXTitle("m (GeV)");
- dmassd->Draw();
-
-
-}
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-void MUONdigit (Int_t evNumber1=0, Int_t evNumber2=0, Int_t ibg=0, Int_t bgr=10)
-{
- //////////////////////////////////////
- // //
- // ROOT macro for ALICE Dimuon Arm: //
- // Digitization //
- // //
- //////////////////////////////////////
- //
- // Adds the tree TD for digits (charges deposited on pads)
- // to the ROOT file "galice.root"
- // containing the signal hit coordinates from simulation (tree TH).
- // Eventually (argument "ibg"), background hits are also taken into account,
- // from the ROOT file "bg.root".
- //
- // Arguments:
- // evNumber1 = first event number to digitize in file "galice.root"
- // evNumber2 = last event number to digitize in file "galice.root"
- // ibg = 0 if no background hits to be taken into account;
- // = 1 if background hits to be taken into account
- // in file "bg.root"
- // bgr : used only if "ibg" = 1
- // = number of events in the background file "bg.root";
- // the signal events are divided into "bgr" successive parts,
- // all events of each part are associated
- // with the same background event,
- // starting with event number 0,
- // incrementing it by 1 for the next part.
- // Strictly speaking, "bgr" can be smaller than
- // the number of events in the background file,
- // in which case one will only use
- // the first "bgr" events of the background file.
- // But it SHOULD NOT BE LARGER THAN
- // THE NUMBER OF EVENTS IN THE BACKGROUND FILE.
- //
- // Input file(s):
- // "galice.root" for signal
- // "bg.root" for background (used only if "ibg" = 1)
- //
- // Output file:
- // "galice.root"
- //
- //__________________________________________________________________________
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (file) file->Close();
- file = new TFile("galice.root","UPDATE");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) cout<<"AliRun object found on file"<<endl;
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-
- AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
- if (pMUON) {
-// creation
- AliMUONMerger* merger = new AliMUONMerger();
-// configuration
- if (ibg) {
- merger->SetMode(ibg);
- merger->SetBackgroundFileName("bg.root");
- }
- // pass
- pMUON->SetMerger(merger);
- }
-// Action !
-//
-// Loop over events
-//
- for (int nev=evNumber1; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev " << nev <<endl;
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
- Int_t nbgr_ev = Int_t(nev*bgr/(evNumber2+1));
-
- if (ibg) {
- merger->SetBackgroundEventNumber(nbgr_ev);
- }
-
- gAlice->SDigits2Digits("MUON");
-
-// Tree writing was moved to AliMUON::SDigits2Digits();
-
- } // event loop
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-#include "iostream.h"
-
-void MUONdoubles (Int_t evNumber1=0,Int_t evNumber2=0)
-{
- {
- hprof = new TProfile("hprof","Profile dmin vs r",24,0,120,0,50);
- }
- hdist1 = new TH1F("hdist1","distance",100,0,10);
- hdist2 = new TH1F("hdist2","distance",100,0,10);
- Float_t a1=3.7/TMath::Sqrt(4);
- Float_t a2=26/TMath::Sqrt(4);
- for (Int_t i=1; i<10000; i++) {
- Float_t x1,x2,y1,y2;
- x1 = (2*gRandom->Rndm(i)-1.)*a1;
- x2 = (2*gRandom->Rndm(i)-1.)*a1;
- y1 = (2*gRandom->Rndm(i)-1.)*a1;
- y2 = (2*gRandom->Rndm(i)-1.)*a1;
- Float_t d=TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
- hdist1->Fill(d,1);
- }
- Float_t xh[100], yh[100];
- for (Int_t k=0; k<1000; k++) {
- for (Int_t i=0; i<100; i++) {
- xh[i] = (2*gRandom->Rndm(i)-1.)*a2;
- yh[i] = (2*gRandom->Rndm(i)-1.)*a2;
- }
- Float_t x1,y1;
-
- for (Int_t i=0; i<10; i++) {
- Float_t dmin=1000;
- x1 = (2*gRandom->Rndm(i)-1.)*a2;
- y1 = (2*gRandom->Rndm(i)-1.)*a2;
- for (Int_t j=0; j<100; j++) {
- Float_t d=TMath::Sqrt((x1-xh[j])*(x1-xh[j])+(y1-yh[j])*(y1-yh[j]));
- if (d<dmin) dmin=d;
- }
- hdist2->Fill(dmin,1);
- }
- }
-
-
-
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root","UPDATE");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-//
-// Set reconstruction models
-//
-// Get pointers to Alice detectors and Digits containers
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
-//
-// Book some histograms
- TH1F *Hcentre[10], *Hall[10], *Hcfrac[10], *Htail[10], *Htfrac[10];
- TH1F *Hdcut[10], *Hdall[10], *Hdfrac[10];
-
- for (Int_t i=0; i<10; i++) {
- Hcentre[i] = new TH1F("Hcentre","Hit Distribution",26,0,260);
- Hcentre[i]->SetFillColor(0);
- Hcentre[i]->SetXTitle("R (cm)");
-
- Hall[i] = new TH1F("Hall","Hit Distribution",26,0,260);
- Hall[i]->SetFillColor(0);
- Hall[i]->SetXTitle("R (cm)");
-
- Hcfrac[i] = new TH1F("Hcfrag","Hit Distribution",26,0,260);
- Hcfrac[i]->SetFillColor(0);
- Hcfrac[i]->SetXTitle("R (cm)");
-
- Htail[i] = new TH1F("Htail","Hit Distribution",26,0,260);
- Htail[i]->SetFillColor(0);
- Htail[i]->SetXTitle("R (cm)");
-
- Htfrac[i] = new TH1F("Htfrag","Hit Distribution",26,0,260);
- Htfrac[i]->SetFillColor(0);
- Htfrac[i]->SetXTitle("R (cm)");
-
- Hdall[i] = new TH1F("Hdall","Hit Distribution",26,0,260);
- Hdall[i]->SetFillColor(0);
- Hdall[i]->SetXTitle("R (cm)");
-
- Hdcut[i] = new TH1F("Hdcut","Hit Distribution",26,0,260);
- Hdcut[i]->SetFillColor(0);
- Hdcut[i]->SetXTitle("R (cm)");
-
- Hdfrac[i] = new TH1F("Hdfrac","Hit Distribution",26,0,260);
- Hdfrac[i]->SetFillColor(0);
- Hdfrac[i]->SetXTitle("R (cm)");
-
- }
-
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev " << nev <<endl;
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- cout<<"ntracks "<<ntracks<<endl;
-
- Int_t nbytes = 0;
-
-
- TClonesArray *Particles = gAlice->Particles();
- TTree *TD = gAlice->TreeD();
- Int_t nent=TD->GetEntries();
- printf("Found %d entries in the tree (must be one per cathode per event!)\n",nent);
- Float_t x[2000];
- Float_t y[2000];
- Int_t nhit=0;
-
- if (MUON) {
-//
-// Loop on chambers and on cathode planes
-//
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-
-//
-// Loop over Hits
-//
- for (i=0; i<1000; i++) {
- Float_t r =100.*gRandom->Rndm(i);
- Float_t phi=2.*TMath::Pi()*gRandom->Rndm(i);
- Float_t xr=r*TMath::Sin(phi);
- Float_t yr=r*TMath::Cos(phi);
- Hdall[0]->Fill(r,1.);
- Float_t dmin=100000.;
-
- if (i==0) {
- for (Int_t track=0; track<ntracks; track++) {
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->fChamber; // chamber number
- if (nch!=1) continue;
-
- x[nhit] = mHit->fX; // x-pos of hit
- y[nhit] = mHit->fY; // y-pos
- nhit++;
- } // hit loop
- } // track loop
- } else {
- for (Int_t nh=0; nh<nhit; nh++) {
- Float_t d=TMath::Sqrt((x[nh]-xr)*(x[nh]-xr)+(y[nh]-yr)*(y[nh]-yr));
-// printf ("\n r,d %f %f",r,d);
- Float_t dx=TMath::Abs(x[nh]-xr);
- Float_t dy=TMath::Abs(y[nh]-yr);
- if (d < dmin) dmin=d;
-
- if (dx<0.5 && dy < 0.75) Hdcut[0]->Fill(r,1.);
- } // hit loop
- } // i==0
- if (r>20) hprof->Fill(r,dmin,1);
- } // random loop
- Int_t icat=0;
- gAlice->ResetDigits();
- gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
-
- for (Int_t ich=0;ich<1;ich++) {
- AliMUONChamber iChamber= MUON->Chamber(ich);
- TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
- if (MUONdigits == 0) continue;
- //
- // Get ready the current chamber stuff
- //
- AliMUONResponse* response = iChamber.GetResponseModel();
- AliMUONSegmentation* seg = iChamber.GetSegmentationModel(icat+1);
- HitMap = new AliMUONHitMapA1(seg, MUONdigits);
- HitMap->FillHits();
- Int_t nxmax=seg->Npx();
- Int_t nymax=seg->Npy();
-//
-// generate random positions on the chamber
-//
- for (Int_t i=0; i<4000; i++) {
- Int_t ix = 2*Int_t((gRandom->Rndm(i)-0.5)*nxmax);
- Int_t iy = 2*Int_t((gRandom->Rndm(i)-0.5)*nymax);
-
-// test for centre overlap
-//
- Float_t xp,yp;
- seg->GetPadCxy(ix,iy,xp,yp);
- Float_t r=TMath::Sqrt(xp*xp+yp*yp);
- if (TMath::Abs(xp)>3 && TMath::Abs(yp)>3)
- Hall[ich]->Fill(r,1.);
- if (HitMap->TestHit(ix,iy) != 0) {
- Hcentre[ich]->Fill(r,1.);
- }
-// test for neighbour overlap
-//
- Int_t nn;
- Int_t X[20], Y[20];
- seg->Neighbours(ix,iy,&nn,X,Y);
- Bool_t hit=false;
-
- for (Int_t j=0; j<nn; j++) {
- if (HitMap->TestHit(X[j],Y[j]) != 0) hit=true;
- }
- if (hit && HitMap->TestHit(ix,iy) == 0) Htail[ich]->Fill(r,1.);
- } //random loop
- delete HitMap;
- } // chamber loop
- } // event loop
- } // if MUON
-
- for (Int_t ich=0; ich<10; ich++) {
- Hcfrac[ich]->Divide(Hcentre[ich],Hall[ich]);
- Htfrac[ich]->Divide(Htail[ich],Hall[ich]);
- Hdfrac[ich]->Divide(Hdcut[ich],Hdall[ich]);
- }
-
- TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
- pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(0);
- pad12->SetFillColor(0);
- pad13->SetFillColor(0);
- pad14->SetFillColor(0);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- Hcentre[0]->Draw();
-
- pad12->cd();
- Hall[0]->Draw();
-
- pad13->cd();
- Hcfrac[0]->Draw();
-
- pad14->cd();
- Htfrac[0]->Draw();
-
- TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
- pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
- pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
- pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
- pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
- pad21->SetFillColor(0);
- pad22->SetFillColor(0);
- pad23->SetFillColor(0);
- pad24->SetFillColor(0);
- pad21->Draw();
- pad22->Draw();
- pad23->Draw();
- pad24->Draw();
-
- pad21->cd();
- hdist1->Draw();
-
- pad22->cd();
- hdist2->Draw();
-
- pad23->cd();
- Hdfrac[0]->Draw();
-
- pad24->cd();
- hprof->Draw();
- file->Close();
-}
-
+++ /dev/null
-void MUONhitrectest (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-// Create some histograms
-
- TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
- TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
- TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
- TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
- TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
- TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
- TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
- TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
- TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
- TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
-
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev " << nev <<endl;
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- cout<<"ntracks "<<ntracks<<endl;
-
- Int_t nbytes = 0;
-
-
-
-// Get pointers to Alice detectors and Digits containers
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- TClonesArray *Particles = gAlice->Particles();
- TTree *TR = gAlice->TreeR();
- Int_t nent=TR->GetEntries();
-
- if (MUON) {
- Float_t xc[2000], yc[2000], itrc[2000];
- TClonesArray *MUONrawclust = MUON->RawClustAddress(0);
- nbytes += TR->GetEvent(0);
- Int_t nrawcl = MUONrawclust->GetEntries();
- AliMUONRawCluster *mRaw;
- for (Int_t iraw=0; iraw < nrawcl; iraw++) {
- mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
- xc[iraw]=mRaw->fX[1];
- yc[iraw]=mRaw->fY[0];
- itrc[iraw]=mRaw->fTracks[1];
- } // cluster
-
- gAlice->ResetHits();
- for (Int_t itrack=0; itrack<ntracks ;itrack++)
- {
- Int_t nbytes += TH->GetEvent(itrack);
- Int_t nhit=MUON->Hits()->GetEntriesFast();
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->Chamber(); // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- Int_t ip = mHit->Particle();
- if (ip != kMuonPlus && ip != kMuonMinus) continue;
- if (nch != 1) continue;
- Float_t dmin=1.e8;
- Int_t imin=-1;
-//
-// Loop over clusters
- Int_t npos=0;
-
- for (Int_t i=0; i<nrawcl; i++) {
- Float_t dx = xc[i]-x;
- Float_t dy = yc[i]-y;
- Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
- if (dr<dmin) {
- imin=i;
- dmin=dr;
- }
- if (TMath::Abs(dy) < 0.05 && TMath::Abs(dx) < 0.2 ) npos++;
- }
- mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
- Int_t track = mRaw->fTracks[1];
- Float_t xrec = mRaw->fX[1];
- Float_t yrec = mRaw->fY[0];
- Float_t x_res = xrec-x;
- Float_t y_res = yrec-y;
- if (y_res > 0.05 ) printf("\n track %d %f\n", itrack, y_res);
-
- hresx->Fill(x_res,1.);
- hresy->Fill(y_res,1.);
- hpos->Fill(Float_t(npos), 1.);
- hresym->Fill(y_res*1.e4,1.);
- if (npos == 0) {
- printf("No cluster %d %f %f %d\n",
- itrack, x, y, nhit);
- h21->Fill(x,y,1.);
- }
- } // hits
- } // tracks
- } // end if MUON
- } // event loop
- TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
- c1->Divide(2,2);
- c1->cd(1);
- hresx->SetFillColor(42);
- hresx->SetXTitle("xrec-x");
- hresx->Draw();
-
- c1->cd(2);
- hresy->SetFillColor(42);
- hresy->SetXTitle("yrec-y");
- hresy->Draw();
-
- c1->cd(3);
- hpos->SetFillColor(42);
- hpos->SetXTitle("Possibilities");
- hpos->Draw();
-
- c1->cd(4);
- hresym->SetFillColor(42);
- hresym->SetXTitle("yrec-y");
- hresym->Fit("gaus");
- hresym->Draw();
-
- TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
- c2->Divide(2,2);
-
- c2->cd(1);
- h21->SetFillColor(42);
- h21->SetXTitle("x");
- h21->SetYTitle("y");
- h21->Draw();
-
- c2->cd(2);
- hmult->SetFillColor(42);
- hmult->SetXTitle("multiplicity");
- hmult->Draw();
-
- c2->cd(3);
- hchi1->SetFillColor(42);
- hchi1->SetXTitle("chi2");
- hchi1->Draw();
-
- c2->cd(4);
- hchi2->SetFillColor(42);
- hchi2->SetXTitle("chi2");
- hchi2->Draw();
-
-}
-
-
-
+++ /dev/null
-void MUONhits (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
- Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
- 56.1, 56.1, 59.6, 59.6};
- Float_t rmax[14] = {81.6, 81.6, 109.3, 109.3, 154.4, 154.4, 197.8,
- 197.8, 229.5, 229.5, 254., 254, 270., 270.};
-
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
- file->ls();
-//
- TDatabasePDG* DataBase = new TDatabasePDG();
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-// Create some histograms
-
- TH1F *Hits[14], *HitDensity[14];
- TH1F *Hits_g[14], *HitDensity_g[14];
- TH1F *Hits_n[14], *HitDensity_n[14];
- TH1F *Mom[14], *Mom_g[14], *Mom_n[14];
-
- for (Int_t i=0; i<14; i++) {
- Hits[i] = new TH1F("hits1","Hit Distribution",30,0,300);
- Hits[i]->SetFillColor(0);
- Hits[i]->SetXTitle("R (cm)");
-
- HitDensity[i] = new TH1F("dhits1","Hit Density",30,0,300);
- HitDensity[i]->SetFillColor(0);
- HitDensity[i]->SetXTitle("R (cm)");
-
- Hits_g[i] = new TH1F("hits1","Hit Distribution",30,0,300);
- Hits_g[i]->SetFillColor(0);
- Hits_g[i]->SetXTitle("R (cm)");
-
- HitDensity_g[i] = new TH1F("dhits1","Hit Density",30,0,300);
- HitDensity_g[i]->SetFillColor(0);
- HitDensity_g[i]->SetXTitle("R (cm)");
-
- Mom[i] = new TH1F("mom","Energy",70,-6,1);
- Mom[i]->SetFillColor(0);
- Mom[i]->SetXTitle("E (GeV)");
-
- Mom_g[i] = new TH1F("mom","Energy",70,-6,1);
- Mom_g[i]->SetFillColor(0);
- Mom_g[i]->SetXTitle("E (GeV)");
-
- Mom_n[i] = new TH1F("mom","Energy",70,-6,1);
- Mom_n[i]->SetFillColor(0);
- Mom_n[i]->SetXTitle("E (GeV)");
-
- Hits_n[i] = new TH1F("hits1","Hit Distribution",30,0,300);
- Hits_n[i]->SetFillColor(0);
- Hits_n[i]->SetXTitle("R (cm)");
-
- HitDensity_n[i] = new TH1F("dhits1","Hit Density",30,0,300);
- HitDensity_n[i]->SetFillColor(0);
- HitDensity_n[i]->SetXTitle("R (cm)");
- }
-
- TH1F *theta = new TH1F("theta","Theta distribution",180,0,180);
-
- TH1F *emult = new TH1F("emult","Event Multiplicity",100,0,1000);
- AliMUONChamber* iChamber;
- AliSegmentation* seg;
-
-//
-// Loop over events
-//
- Float_t dpx[2], dpy[2];
- Int_t mdig =0;
- Int_t nhit=0;
-
- for (Int_t nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- printf("\n track %d %d \n ", nev, MUON);
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-//
-// Loop over events
-//
- Int_t EvMult=0;
-
- for (Int_t track=0; track<ntracks;track++) {
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (MUON) {
-//
-// Loop over hits
-//
- Int_t NperCh=0;
-
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->Chamber(); // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- Float_t Eloss = mHit->Eloss();
- Float_t Theta = mHit->Theta();
- Float_t Particle = mHit->Particle();
- Float_t P = mHit->Momentum();
- TParticlePDG* Part = DataBase->GetParticle(Particle);
- Double_t mass = Part->Mass();
-
- if (nch >13) continue;
- if (nch ==1) EvMult++;
- if (mHit->Age() > 5.e-6) continue;
-
- Float_t r=TMath::Sqrt(x*x+y*y);
- if (nch ==0) continue;
-
- if (r < rmin[nch-1]) continue;
- if (r > rmax[nch-1]) continue;
-
- Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
-
- Float_t Ekin=TMath::Sqrt(P*P+mass*mass)-mass;
- if (Particle == 2112) {
- Hits_n[nch-1]->Fill(r,(float) 1);
- HitDensity_n[nch-1]->Fill(r,wgt);
- Mom_n[nch-1]->Fill(TMath::Log10(Ekin),1);
- } else if (Particle == 22) {
- Hits_g[nch-1]->Fill(r,(float) 1);
- HitDensity_g[nch-1]->Fill(r,wgt);
- Mom_g[nch-1]->Fill(TMath::Log10(Ekin),1);
- } else {
- Hits[nch-1]->Fill(r,(float) 1);
- HitDensity[nch-1]->Fill(r,wgt);
- Mom[nch-1]->Fill(TMath::Log10(Ekin),1);
- }
- } // hit loop
- } // if MUON
- } // track loop
- }
-
-//Create a canvas, set the view range, show histograms
- Int_t k;
- TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
- c1->Divide(2,4);
- for (k=0; k<7; k++) {
- c1->cd(k+1);
- HitDensity[2*k]->Draw();
- }
-
- TCanvas *c2 = new TCanvas("c2","Hit Densities (gamma)",400,10,600,700);
- c2->Divide(2,4);
- for (k=0; k<7; k++) {
- c2->cd(k+1);
- HitDensity_g[2*k]->Draw();
- }
-
- TCanvas *c3 = new TCanvas("c3","Hit Densities (neutron)",400,10,600,700);
- c3->Divide(2,4);
- for (k=0; k<7; k++) {
- c3->cd(k+1);
- HitDensity_n[2*k]->Draw();
- }
-
- TCanvas *c4 = new TCanvas("c4","Energy (charged)",400,10,600,700);
- c4->Divide(2,4);
- for (k=0; k<7; k++) {
- c4->cd(k+1);
- Mom[2*k]->Draw();
- }
-
- TCanvas *c5 = new TCanvas("c5","Energy (gamma)",400,10,600,700);
- c5->Divide(2,4);
- for (k=0; k<7; k++) {
- c5->cd(k+1);
- Mom_g[2*k]->Draw();
- }
-
- TCanvas *c6 = new TCanvas("c6","Energy (gamma)",400,10,600,700);
- c6->Divide(2,4);
- for (k=0; k<7; k++) {
- c6->cd(k+1);
- Mom_n[2*k]->Draw();
- }
-}
-
-
-
-
-
-
-
+++ /dev/null
-// macro to make invariant mass plots
-// for combinations of 2 muons with opposite charges,
-// from root file "MUONtrackReco.root" containing the result of track reconstruction,
-// generated by the macro "MUONrecoNtuple.C".
-// Histograms are stored on the "MUONmassPlot.root" file.
-// A model for macros using the Ntuple in the file "MUONtrackReco.root"
-// may be found in "MUONtrackRecoModel.C":
-// it has been obtained by reading with Root a file "MUONtrackReco.root",
-// and executing the command:
-// MUONtrackReco->MakeCode("MUONtrackRecoModel.C")
-
-// Arguments:
-// FirstEvent (default 0)
-// LastEvent (default 0)
-// ResType (default 553)
-// 553 for Upsilon, anything else for J/Psi
-// Chi2Cut (default 100)
-// to keep only tracks with chi2 per d.o.f. < Chi2Cut
-// PtCutMin (default 1)
-// to keep only tracks with transverse momentum > PtCutMin
-// PtCutMax (default 10000)
-// to keep only tracks with transverse momentum < PtCutMax
-// massMin (default 9.17 for Upsilon)
-// & massMax (default 9.77 for Upsilon)
-// to calculate the reconstruction efficiency for resonances with invariant mass
-// massMin < mass < massMax.
-
-// compile MUONrecoNtuple.C and load shared library in this macro
-// Add parameters and histograms for analysis
-
-void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77)
-{
- cout << "MUONmassPlot " << endl;
- cout << "FirstEvent " << FirstEvent << endl;
- cout << "LastEvent " << LastEvent << endl;
- cout << "ResType " << ResType << endl;
-// cout << "Nsig " << Nsig << endl;
- cout << "Chi2Cut " << Chi2Cut << endl;
- cout << "PtCutMin " << PtCutMin << endl;
- cout << "PtCutMax " << PtCutMax << endl;
- cout << "massMin " << massMin << endl;
- cout << "massMax " << massMax << endl;
-
-
- //////////////////////////////////////////////////////////
- // This file has been automatically generated
- // (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
- // from TTree MUONtrackReco/MUONtrackReco
- // found on file: MUONtrackReco.root
- //////////////////////////////////////////////////////////
-
-
- //Reset ROOT and connect tree file
- gROOT->Reset();
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER");
-// gROOT->LoadMacro("loadlibs.C");
-// loadlibs();
-// compile MUONrecoNtuple and load shared library
- gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf");
- }
-
- TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
- if (!f) {
- f = new TFile("MUONtrackReco.root");
- }
- TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
- MUONtrackReco->SetMakeClass(1);
-
-//Declaration of leaves types
- Int_t fEvent;
- UInt_t fUniqueID;
- UInt_t fBits;
- Int_t Tracks_;
- Int_t Tracks_fCharge[500];
- Float_t Tracks_fPxRec[500];
- Float_t Tracks_fPyRec[500];
- Float_t Tracks_fPzRec[500];
- Float_t Tracks_fZRec[500];
- Float_t Tracks_fZRec1[500];
- Int_t Tracks_fNHits[500];
- Float_t Tracks_fChi2[500];
- Float_t Tracks_fPxGen[500];
- Float_t Tracks_fPyGen[500];
- Float_t Tracks_fPzGen[500];
- UInt_t Tracks_fUniqueID[500];
- UInt_t Tracks_fBits[500];
-
- //Set branch addresses
- //MUONtrackReco->SetBranchAddress("Header",&Header);
- MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
- MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
- MUONtrackReco->SetBranchAddress("fBits",&fBits);
-// MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
- MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
- MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
- MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
- MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
- MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
- MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
- MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
- MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
- MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
- MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
- MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
- MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
- MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
- MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
-
-// This is the loop skeleton
-// To read only selected branches, Insert statements like:
-// MUONtrackReco->SetBranchStatus("*",0); // disable all branches
-// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
-
- Int_t nentries = MUONtrackReco->GetEntries();
-
- Int_t nbytes = 0;
- // for (Int_t i=0; i<nentries;i++) {
- // nbytes += MUONtrackReco->GetEntry(i);
- // }
-
- /////////////////////////////////////////////////////////////////
- // Here comes the specialized part for MUONmassPlot
- /////////////////////////////////////////////////////////////////
-
- // File for histograms and histogram booking
- TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
- TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
- TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
- TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
- TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
- if (ResType == 553) {
- TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
- } else {
- TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
- }
-
- TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
- TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5.,-2);
- TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
- TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
- Int_t EventInMass = 0;
- Float_t muonMass = 0.105658389;
- Float_t UpsilonMass = 9.46037;
- Float_t JPsiMass = 3.097;
- // Loop over events
- for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
- // get current event
- Int_t NumberOfTrack = 0;
- cout <<" event: " << event <<endl;
- cout << " number of tracks: " << Tracks_ <<endl;
- nbytes += MUONtrackReco->GetEntry(event);
- // loop over all reconstructed tracks (also first track of combination)
- for (Int_t t1 = 0; t1 < Tracks_; t1++) {
- // transverse momentum
- Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
-
- // total momentum
- Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] );
- // Rapidity
- Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass);
-
- // chi2 per d.o.f.
- Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
- printf(" px %f py %f pz %f NHits %d Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1);
- // condition for good track (Chi2Cut and PtCut)
- if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
- // fill histos hPtMuon and hChi2PerDof
- NumberOfTrack++;
- hPtMuon->Fill(pt1);
- hPMuon->Fill(p1);
- hChi2PerDof->Fill(ch1);
- hRapMuon->Fill(rapMuon1);
- // loop over second track of combination
- for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
- // transverse momentum
- Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
- // chi2 per d.o.f.
- Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
- // condition for good track (Chi2Cut and PtCut)
- if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) {
- // condition for opposite charges
- if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
- // invariant mass
- Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass);
- // fill histos hInvMassAll and hInvMassRes
- hInvMassAll->Fill(invMass);
- hInvMassRes->Fill(invMass);
- Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2];
- Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2];
- Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2];
- if (invMass > massMin && invMass < massMax) {
- EventInMass++;
- Float_t rapRes = 0;
- if (ResType == 553) {
- rapRes = rapParticle(pxRes,pyRes,pzRes,UpsilonMass);
- } else {
- rapRes = rapParticle(pxRes,pyRes,pzRes,JPsiMass);
- }
- hRapResonance->Fill(rapRes);
- hPtResonance->Fill(TMath::Sqrt(pxRes*pxRes+pyRes*pyRes));
- }
-
- } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
- } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
- } // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
- } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
- } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
- hNumberOfTrack->Fill(NumberOfTrack);
- } // for (Int_t event = FirstEvent;
-
- histoFile->Write();
- histoFile->Close();
-
- cout << "MUONmassPlot " << endl;
- cout << "FirstEvent " << FirstEvent << endl;
- cout << "LastEvent " << LastEvent << endl;
- cout << "ResType " << ResType << endl;
-// cout << "Nsig " << Nsig << endl;
- cout << "Chi2Cut " << Chi2Cut << endl;
- cout << "PtCutMin " << PtCutMin << endl;
- cout << "PtCutMax " << PtCutMax << endl;
- cout << "massMin " << massMin << endl;
- cout << "massMax " << massMax << endl;
- cout << "EventInMass " << EventInMass << endl;
-}
-
-
-Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
-{
- // return invariant mass for particle 1 & 2 whose masses are equal to muonMass
-
- Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
- 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;
-}
+++ /dev/null
-#include "iostream.h"
-
-void MUONmultest (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-// Create some histograms
-
- TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
- TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
- TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
- TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
- TH1F *hresx = new TH1F("hresx","Residuals",100,-4,4);
- TH1F *hresy = new TH1F("hresy","Residuals",100,-.1,.1);
- TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nev " << nev <<endl;
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- cout<<"ntracks "<<ntracks<<endl;
-
- Int_t nbytes = 0;
-
- AliMUONRawCluster *mRaw;
-
-// Get pointers to Alice detectors and Digits containers
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- TClonesArray *Particles = gAlice->Particles();
- TTree *TR = gAlice->TreeR();
- Int_t nent=TR->GetEntries();
- printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
- if (MUON) {
- for (Int_t ich=0;ich<1;ich++) {
- TClonesArray *MUONrawclust = MUON->RawClustAddress(ich);
- TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
- for (Int_t icat=1; icat<2; icat++) {
- MUON->ResetRawClusters();
- nbytes += TR->GetEvent(icat);
- Int_t nrawcl = MUONrawclust->GetEntries();
- printf(
- "Found %d raw-clusters for cathode %d in chamber %d \n"
- ,nrawcl,icat,ich+1);
-
-
- gAlice->ResetDigits();
-
- Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
- gAlice->TreeD()->GetEvent(nent-2+icat-1);
- Int_t ndigits = MUONdigits->GetEntriesFast();
- printf(
- "Found %d digits for cathode %d in chamber %d \n"
- ,ndigits,icat,ich+1);
-
-
- Int_t TotalMult =0;
-
- for (Int_t iraw=0; iraw < nrawcl; iraw++) {
- mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
- Int_t mult=mRaw->fMultiplicity;
- h1->Fill(mult,float(1));
- TotalMult+=mult;
- Int_t itrack=mRaw->fTracks[0];
- Float_t xrec=mRaw->fX;
- Float_t yrec=mRaw->fY;
- Float_t R=TMath::Sqrt(xrec*xrec+yrec*yrec);
- if (R > 55.2) continue;
- if (itrack ==1) continue;
- Float_t res[2];
- Int_t nres=0;
- nbytes=0;
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(itrack);
-
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->fChamber; // chamber number
- Float_t x = mHit->fX; // x-pos of hit
- Float_t y = mHit->fY; // y-pos
- if (nch==(ich+1)){
- hresx->Fill(xrec-x,float(1));
- hresy->Fill(yrec-y,float(1));
- if ((yrec-y)*1e4 <500 )
- hresym->Fill((yrec-y)*1e4,float(1));
- if (TMath::Abs(yrec-y)>.02) {
- h22->Fill(mRaw->fX,mRaw->fY,float(1));
- hmult->Fill(mult,float(1));
- }
- } // chamber
- } //hit
- } //iraw
- printf("Total Cluster Multiplicity %d \n" , TotalMult);
- } // icat
- } // ich
- } // end if MUON
- } // event loop
- TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
- pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(11);
- pad12->SetFillColor(11);
- pad13->SetFillColor(11);
- pad14->SetFillColor(11);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- hresx->SetFillColor(42);
- hresx->SetXTitle("xrec-x");
- hresx->Draw();
-
- pad12->cd();
- hresy->SetFillColor(42);
- hresy->SetXTitle("yrec-y");
- hresy->Draw();
-
- pad13->cd();
- h1->SetFillColor(42);
- h1->SetXTitle("multiplicity");
- h1->Draw();
-
- pad14->cd();
- hresym->SetFillColor(42);
- hresym->SetXTitle("yrec-y");
- hresym->Fit("gaus");
- hresym->Draw();
- TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
- pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
- pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
- pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
- pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
- pad21->SetFillColor(11);
- pad22->SetFillColor(11);
- pad23->SetFillColor(11);
- pad24->SetFillColor(11);
- pad21->Draw();
- pad22->Draw();
- pad23->Draw();
- pad24->Draw();
-
- pad21->cd();
- h22->SetFillColor(42);
- h22->SetXTitle("x");
- h22->SetYTitle("y");
- h22->Draw();
-
- pad22->cd();
- hmult->SetFillColor(42);
- hmult->SetXTitle("multiplicity");
- hmult->Draw();
-
-}
-
-
-
+++ /dev/null
-void MUONocc (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
- file->ls();
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-// Create some histograms
-
- TH1F *Hits[10], *HitDensity[10], *Occ0[10],*Occ1[10], *Mult[10];
- for (Int_t i=0; i<10; i++) {
- Hits[i] = new TH1F("hits1","Hit Distribution",26,0,260);
- Hits[i]->SetFillColor(0);
- Hits[i]->SetXTitle("R (cm)");
-
- HitDensity[i] = new TH1F("dhits1","Hit Density Distribution",26,0,260);
- HitDensity[i]->SetFillColor(0);
- HitDensity[i]->SetXTitle("R (cm)");
-
- Occ0[i] = new TH1F("occ0","Occupancy Density",26,0,260);
- Occ0[i] -> SetFillColor(0);
- Occ0[i] -> SetXTitle("R (cm)");
- Occ1[i] = new TH1F("occ1","Occupancy",26,0,260);
- Occ1[i] -> SetFillColor(0);
- Occ1[i] -> SetXTitle("R (cm)");
- Occ1[i] -> SetYTitle("Occupancy");
- Mult[i] = new TH1F("mult","Mult distribution",26,0,260);
- Mult[i] -> SetFillColor(0);
- Mult[i] -> SetXTitle("R (cm)");
- }
-
-
- TH1F *theta = new TH1F("theta","Theta distribution",180,0,180);
-
- TH1F *emult = new TH1F("emult","Event Multiplicity",100,0,1000);
- AliMUONChamber* iChamber;
- AliMUONSegmentation* seg;
-
-//
-// Loop over events
-//
- Float_t dpx[2], dpy[2];
- Int_t mdig =0;
- Int_t nhit=0;
-
- for (Int_t nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- printf("\n track %d %d \n ", nev, MUON);
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-//
-// Loop over events
-//
- Int_t EvMult=0;
-
- for (Int_t track=0; track<ntracks;track++) {
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (MUON) {
-//
-// Loop over hits
-//
- Int_t NperCh=0;
-
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->Chamber(); // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- Float_t Eloss = mHit->Eloss();
- Float_t Theta = mHit->Theta();
- theta->Fill(Theta,(float) 1);
- if (nch >10) continue;
- if (nch ==1) EvMult++;
-//
-//
- iChamber = & MUON->Chamber(nch-1);
- response= iChamber.GetResponseModel();
- seg = iChamber.GetSegmentationModel(1);
- NperCh++;
- Int_t i,j;
- seg->GetPadIxy(x,y,i,j);
- Int_t isec = seg->Sector(i,j);
- if (isec ==1 && nch==ic) nhit++;
-
- Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
- Float_t r=TMath::Sqrt(x*x+y*y);
- Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
- wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
- Hits[nch-1]->Fill(r,(float) 1);
- HitDensity[nch-1]->Fill(r,wgt);
- Occ0[nch-1]->Fill(r,wgt*a);
- } // hit loop
- } // if MUON
- } // track loop
- emult->Fill(Float_t(EvMult), (Float_t) 1);
-
- Int_t iseg=1;
-
- for (Int_t ich=0; ich<10; ich++) {
- iChamber = & MUON->Chamber(ich);
- seg=iChamber.GetSegmentationModel(iseg);
- gAlice->ResetDigits();
- gAlice->TreeD()->GetEvent(iseg);
- TClonesArray *MUONDigits = MUON->DigitsAddress(ich);
- Int_t Ndigits=MUONDigits->GetEntriesFast();
- AliMUONDigit* dig;
- printf("\n Reading %d digits\n", Ndigits);
- if (MUONDigits) {
- for (Int_t ndig=0; ndig<Ndigits; ndig++)
- {
- dig = (AliMUONDigit*)MUONDigits->UncheckedAt(ndig);
- Int_t i=dig->PadX();
- Int_t j=dig->PadY();
- Float_t x,y;
- seg->GetPadCxy(i,j,x,y);
- Int_t isec = seg->Sector(i,j);
- Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
- Float_t r=TMath::Sqrt(x*x+y*y);
-// if (r<25)
-// printf("\n Sector,a %d %f", isec,a);
- if (isec==1) mdig++;
-
- Float_t wgt;
- wgt=1/(2.*10.*TMath::Pi()*r)/(evNumber2+1)*a;
-// Take into account inefficiency due to frames
- wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
- Occ1[ich]->Fill(r,wgt);
- } // digit loop
- Mult[ich]->Divide(Occ1[ich],Occ0[ich]);
- } // chamber loop
- } // if MUONDigits
- } // event loop
-//
- printf("\n hits, digits %d %d\n ", nhit, mdig);
-
-
-
-//Create a canvas, set the view range, show histograms
- TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
- pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(0);
- pad12->SetFillColor(0);
- pad13->SetFillColor(0);
- pad14->SetFillColor(0);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- HitDensity[0]->Draw();
- pad12->cd();
- HitDensity[1]->Draw();
- pad13->cd();
- HitDensity[2]->Draw();
- pad14->cd();
- HitDensity[3]->Draw();
-
- TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
- pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
- pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
- pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
- pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
- pad21->SetFillColor(0);
- pad22->SetFillColor(0);
- pad23->SetFillColor(0);
- pad24->SetFillColor(0);
- pad21->Draw();
- pad22->Draw();
- pad23->Draw();
- pad24->Draw();
-
- pad21->cd();
- HitDensity[4]->Draw();
- pad22->cd();
- HitDensity[5]->Draw();
- pad23->cd();
- HitDensity[6]->Draw();
- pad24->cd();
- HitDensity[7]->Draw();
-
- TCanvas *c3 = new TCanvas("c3","Hit Densities",400,10,600,700);
- pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
- pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
- pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
- pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
- pad31->SetFillColor(0);
- pad32->SetFillColor(0);
- pad33->SetFillColor(0);
- pad34->SetFillColor(0);
- pad31->Draw();
- pad32->Draw();
- pad33->Draw();
- pad34->Draw();
-
- pad31->cd();
- HitDensity[8]->Draw();
- pad32->cd();
- HitDensity[9]->Draw();
-
- TCanvas *c4 = new TCanvas("c4","Occupancies",400,10,600,700);
- pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
- pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
- pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
- pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
- pad41->SetFillColor(0);
- pad42->SetFillColor(0);
- pad43->SetFillColor(0);
- pad44->SetFillColor(0);
- pad41->Draw();
- pad42->Draw();
- pad43->Draw();
- pad44->Draw();
-
-
- pad41->cd();
- Occ1[0]->Draw();
- pad42->cd();
- Occ1[2]->Draw();
- pad43->cd();
- Occ1[4]->Draw();
- pad44->cd();
- Occ1[6]->Scale(1.25);
- Occ1[6]->Draw();
-
- TCanvas *c5 = new TCanvas("c5","Occupancies",400,10,600,700);
- pad51 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
- pad52 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
- pad53 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
- pad54 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
- pad51->SetFillColor(0);
- pad52->SetFillColor(0);
- pad53->SetFillColor(0);
- pad54->SetFillColor(0);
- pad51->Draw();
- pad52->Draw();
- pad53->Draw();
- pad54->Draw();
-
-
- pad51->cd();
- Occ1[8]->Scale(1.25);
-
- Occ1[8]->Draw();
-}
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-void MUONpadtest (Int_t evNumber1=0,Int_t evNumber2=1)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
- const Int_t NpadX = 252; // number of pads on X
- const Int_t NpadY = 374; // number of pads on Y
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-
-// Create some histograms
-
- Int_t xmin=-NpadX/2;
- Int_t xmax= NpadX/2;
- Int_t ymin=-NpadY/2;
- Int_t ymax= NpadY/2;
-
- TH2F *hc = new TH2F("hc","Chamber1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
- TH2F *h = new TH2F("h","Chamber1 hit distribution",100,-100,100,100,-100,100);
- TH1F *charge = new TH1F("charge","Charge distribution",100,0.,1000.);
-
-// Start loop over events
-
- Int_t Nh=0;
- Int_t Nh1=0;
- AliMUON *pMUON = (AliDetector*) gAlice->GetDetector("MUON");
- for (int nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- cout<<"nev "<<nev<<endl;
- cout<<"nparticles "<<nparticles<<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
-// Get pointers to MUON detector and Hits containers
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-
-// Start loop on tracks in the hits containers
-
- // Int_t Nh=0;
- Int_t Nc=0;
- for (Int_t track=0; track<ntracks;track++) {
-// printf("ntracks %d\n",ntracks);
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (pMUON) {
- TClonesArray *Clusters = pMUON->PadHits(); // Cluster branch
- TClonesArray *Hits = pMUON->Hits(); // Hits branch
-// printf("%d %d \n",Clusters,Hits);
- }
- //see hits distribution
- Int_t nhits = Hits->GetEntriesFast();
- if (nhits) Nh+=nhits;
-// printf("nhits %d\n",nhits);
- for (Int_t hit=0;hit<nhits;hit++) {
-// printf("hit# %d\n",hit);
- mHit = (AliMUONHit*) Hits->UncheckedAt(hit);
- Int_t nch = mHit->Chamber(); // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- // Fill the histograms
- if( nch==1) {
- Float_t rhit=TMath::Sqrt(x*x+y*y);
- if( rhit<= 55 ) Nh1+=nhits;
- h->Fill(x,y,(float) 1);
- }
-
-
- // see signal distribution
- for (AliMUONPadHit* mPad =
- (AliMUONPadHit*)pMUON->FirstPad(mHit,Clusters);
- mPad;
- mPad = (AliMUONPadHit*)pMUON->NextPad(Clusters))
- {
- Int_t nhit = mPad->HitNumber(); // hit number
- Int_t qtot = mPad->Q(); // charge
- Int_t ipx = mPad->PadX(); // pad number on X
- Int_t ipy = mPad->PadY(); // pad number on Y
- Int_t iqpad = mPad->QPad(); // charge per pad
-// printf("%d %d %d\n",ipx,ipy,iqpad);
- if (qtot > 0 && nch == 1){
- charge->Fill(float(iqpad),(float) 1);
- }
- if(nch == 1) {
- hc->Fill(ipx,ipy,(float) iqpad);
- } // if rpad
- } // padhits
- } // hitloops
- }
- // cout<<"Nh "<<Nh<<endl;
- }
-//Create a canvas, set the view range, show histograms
- TCanvas *c1 = new TCanvas("c1","Alice MUON pad hits",400,10,600,700);
- hc->SetXTitle("ix (npads)");
- hc->Draw();
- TCanvas *c2 = new TCanvas("c2","Alice MUON hits",400,10,600,700);
- h->SetFillColor(42);
- h->SetXTitle("x (cm)");
- h->Draw();
- TCanvas *c3 = new TCanvas("c3","Charge distribution",400,10,600,700);
- charge->SetFillColor(42);
- charge->SetXTitle("ADC units");
- charge->Draw();
-
-}
+++ /dev/null
-#include "iostream.h"
-
-void MUONraw (Int_t evNumber1=0,Int_t evNumber2=0)
-{
- //////////////////////////////////////
- // //
- // ROOT macro for ALICE Dimuon Arm: //
- // Clusterization of digits //
- // //
- //////////////////////////////////////
- //
- // Adds the tree TR for raw clusters
- // to the ROOT file "galice.root"
- // containing the digits (tree TD).
- //
- // Arguments:
- // evNumber1 = first event number to act on in file "galice.root"
- // evNumber2 = last event number to act on in file "galice.root"
- //
- // Input/output file:
- // "galice.root"
- //
- //__________________________________________________________________________
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!file) file = new TFile("galice.root","UPDATE");
-
-// Get AliRun object from file or create it if not on file
-
- if (!gAlice) {
- gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
- }
-
-// Set reconstruction models
-//
-// Get pointers to Alice detectors and Digits containers
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- RecModel = new AliMUONClusterFinderAZ(0,1); //AZ
- //RecModel = new AliMUONClusterFinderAZ(1,0); //AZ
- for (Int_t i=0; i<10; i++) {
- MUON->SetReconstructionModel(i,(AliMUONClusterFinderVS*)RecModel);
- }
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- gAlice->RunReco("MUON", evNumber1, evNumber2);
-}
-
+++ /dev/null
-// Macro MUONreco.C for testing the C++ reconstruction code
-
-// Arguments:
-// FirstEvent (default 0)
-// LastEvent (default 0)
-// RecGeantHits (1 to reconstruct GEANT hits) (default 0)
-// FileName (for signal) (default "galice.root")
-// BkgGeantFileName (for background),
-// needed only if RecGeantHits = 1 and background to be added
-#include <iostream.h>
-void MUONreco (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
-{
- //
- cout << "MUONreco" << endl;
- cout << "FirstEvent " << FirstEvent << endl;
- cout << "LastEvent " << LastEvent << endl;
- cout << "RecGeantHits " << RecGeantHits << endl;
- cout << "FileName ``" << FileName << "''" << endl;
- cout << "BkgGeantFileName ``" << BkgGeantFileName << "''" << endl;
- // Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
- // Connect the Root Galice file containing Geometry, Kine, Hits
- // and eventually RawClusters
- TFile *file = (TFile*) gROOT->GetListOfFiles()->FindObject(FileName);
- if (!file) {
- printf("\n Creating file %s\n", FileName);
- file = new TFile(FileName);
- }
- else printf("\n File %s found in file list\n", FileName);
-
- // Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*) file->Get("gAlice");
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n Create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
- // Initializations
- // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
- AliMUONEventReconstructor *Reco = new AliMUONEventReconstructor();
- //Reco->SetTrackMethod(2); //AZ - use Kalman
- Reco->SetRecGeantHits(RecGeantHits);
-
- // The right place for changing AliMUONEventReconstructor parameters
- // with respect to the default ones
-// Reco->SetMaxSigma2Distance(100.0);
- //Reco->SetPrintLevel(10);
- //Reco->SetPrintLevel(2);
- //Reco-> SetChamberThicknessInX0(0.05);
- //Reco->SetEfficiency(1.); //0.95);
- cout << "AliMUONEventReconstructor: actual parameters" << endl;
- Reco->Dump();
-
- // Loop over events
- for (Int_t event = FirstEvent; event <= LastEvent; event++) {
- cout << "Event: " << event << endl;
- // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
- Int_t nparticles = gAlice->GetEvent(event);
- cout << "nparticles: " << nparticles << endl;
- // prepare background file and/or event if necessary
- if (RecGeantHits == 1) {
- if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName);
- if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
- }
- Reco->EventReconstruct();
- // Write this event in a tree
- Reco->FillEvent();
- // Dump current event
- Reco->EventDump();
- } // Event loop
-}
+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-/* $Id$ */
-
-
-
-// Macro MUONrecoNtuple.C (TO BE COMPILED)
-// for testing the C++ reconstruction code
-// and producing an Ntuple with reconstructed tracks
-// in output file "MUONtrackReco.root".
-// An example for using the Ntuple is in the macro MUONmassPlot.C
-
-// Arguments:
-// FirstEvent (default 0)
-// LastEvent (default 0)
-// RecGeantHits (1 to reconstruct GEANT hits) (default 0)
-// FileName (for signal) (default "galice.root")
-// BkgGeantFileName (for background),
-// needed only if RecGeantHits = 1 and background to be added
-
-// IMPORTANT NOTICE FOR USERS:
-// under "root" or "root.exe", execute the following commands:
-// 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER -I$ROOTSYS/include")" to get the right path at compilation time
-// 2. ".x loadlibs.C" to load the shared libraries
-// 3. ".x MUONrecoNtuple.C+()" with the right arguments, without forgetting the "+" which implies the compilation of the macro before its execution
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <iostream.h>
-
-#include <TClassTable.h>
-#include <TClonesArray.h>
-#include <TFile.h>
-#include <TParticle.h>
-#include <TROOT.h>
-#include <TTree.h>
-
-#include "AliHeader.h"
-#include "AliRun.h"
-
-#include "AliMUON.h"
-#include "AliMUONEventReconstructor.h"
-#include "AliMUONTrack.h"
-#include "AliMUONTrackHit.h"
-#include "AliMUONTrackParam.h"
-#endif
-
-// Classes for Ntuple ////////////////////////////////////////////////////
-
-class AliMUONTrackRecNtuple : public TObject {
- public:
- // for direct access to members
- Int_t fCharge; // charge of reconstructed track (+/- 1)
- Float_t fPxRec; // Px of reconstructed track at vertex (GeV/c)
- Float_t fPyRec; // Py of reconstructed track at vertex (GeV/c)
- Float_t fPzRec; // Pz of reconstructed track at vertex (GeV/c)
- Float_t fZRec; // Z of reconstructed track at vertex (cm)
- Float_t fZRec1; // Z of reconstructed track at first hit (cm)
- Int_t fNHits; // number of hits
- Float_t fChi2; // chi2 of fit
- Float_t fPxGen; // Px of best compatible generated track at vertex (GeV/c)
- Float_t fPyGen; // Py of best compatible generated track at vertex (GeV/c)
- Float_t fPzGen; // Pz of best compatible generated track at vertex (GeV/c)
- AliMUONTrackRecNtuple(){;} // Constructor
- virtual ~AliMUONTrackRecNtuple(){;} // Destructor
- protected:
- private:
- ClassDef(AliMUONTrackRecNtuple, 1) // AliMUONTrackRecNtuple
- };
-
-class AliMUONHeaderRecNtuple : public TObject {
- public:
- // for direct access to members
- Int_t fEvent; // event number
- AliMUONHeaderRecNtuple(){;} // Constructor
- virtual ~AliMUONHeaderRecNtuple(){;} // Destructor
- protected:
- private:
- ClassDef(AliMUONHeaderRecNtuple, 1) // AliMUONHeaderRecNtuple
- };
-
-ClassImp(AliMUONTrackRecNtuple) // Class implementation in ROOT context
-ClassImp(AliMUONHeaderRecNtuple) // Class implementation in ROOT context
-
- //__________________________________________________________________________
-void AliMUONEventRecNtupleFill(AliMUONEventReconstructor *Reco, Int_t FillWrite = 0)
-{
- // Fill Ntuple for reconstructed event pointed to by "Reco"
- // if "FillWrite" different from -1.
- // Ntuple is created automatically before filling first event.
- // If "FillWrite" = -1, write and close the file
-
- static Bool_t firstTime = kTRUE;
- static TTree *ntuple;
- static AliMUONHeaderRecNtuple *header = new AliMUONHeaderRecNtuple();
- static TClonesArray *recTracks = new TClonesArray("AliMUONTrackRecNtuple",5);
-
- Int_t trackIndex;
- AliMUONTrack *track;
- AliMUONTrackParam *trackParam;
- AliMUONTrackRecNtuple *recTrackNt;
- Double_t bendingSlope, nonBendingSlope, pYZ;
-
- if (FillWrite == -1) {
- printf(">>> Writing Ntuple of reconstructed tracks\n");
- // better to create the file before the Ntuple ????
- TFile *file = new TFile("MUONtrackReco.root","recreate");
- ntuple->Write();
- file->Write();
- file->Close();
- }
-
- if (firstTime) {
- firstTime = kFALSE;
- // first call: create tree for Ntuple...
- printf(">>> Creating Ntuple of reconstructed tracks\n");
- ntuple = new TTree("MUONtrackReco", "MUONtrackReco");
- ntuple->Branch("Header","AliMUONHeaderRecNtuple", &header);
- ntuple->Branch("Tracks", &recTracks);
- }
-
- // header
-
- header->fEvent = gAlice->GetHeader()->GetEvent();
-
- TClonesArray *recoTracksPtr = Reco->GetRecTracksPtr();
- recoTracksPtr->Compress(); // for simple loop without "Next" since no hole
- recTracks->Clear(); // to reset the TClonesArray of tracks to be put in the ntuple
- // Loop over reconstructed tracks
- for (trackIndex = 0; trackIndex < Reco->GetNRecTracks(); trackIndex++) {
- track = (AliMUONTrack*) ((*recoTracksPtr)[trackIndex]);
- recTrackNt = (AliMUONTrackRecNtuple*)
- new ((*recTracks)[trackIndex]) AliMUONTrackRecNtuple();
- // track parameters at Vertex
- trackParam = track->GetTrackParamAtVertex();
- recTrackNt->fCharge =
- Int_t(TMath::Sign(1., trackParam->GetInverseBendingMomentum()));
- bendingSlope = trackParam->GetBendingSlope();
- nonBendingSlope = trackParam->GetNonBendingSlope();
- pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
- recTrackNt->fPzRec = - pYZ / TMath::Sqrt(1.0 + bendingSlope * bendingSlope); // spectro. (z<0)
- recTrackNt->fPxRec = recTrackNt->fPzRec * nonBendingSlope;
- recTrackNt->fPyRec = recTrackNt->fPzRec * bendingSlope;
- recTrackNt->fZRec = trackParam->GetZ();
- // track parameters at first hit
- trackParam = ((AliMUONTrackHit*)
- (track->GetTrackHitsPtr()->First()))->GetTrackParam();
- recTrackNt->fZRec1 = trackParam->GetZ();
- // chi2
- recTrackNt->fChi2 = track->GetFitFMin();
- // number of hits
- recTrackNt->fNHits = track->GetNTrackHits();
- printf("test> Px %f Py %f Pz %f \n", recTrackNt->fPxRec, recTrackNt->fPyRec, recTrackNt->fPzRec);
- // track parameters at vertex of best compatible generated track:
- // in fact muon with the right charge
- TTree* mtreeK=gAlice->TreeK();
- TBranch *brparticle = mtreeK->GetBranch("Particles");
- Int_t nPart = brparticle->GetEntries();
- TParticle *particle = new TParticle();
- mtreeK->SetBranchAddress("Particles",&particle);
- for (Int_t iPart = 0; iPart < nPart; iPart++) {
- brparticle->GetEntry(iPart);
- //cout << "Code Particle: " << particle->GetPdgCode() << "\n";
- if ((particle->GetPdgCode() * recTrackNt->fCharge) == -13) {
- recTrackNt->fPxGen = particle->Px();
- recTrackNt->fPyGen = particle->Py();
- recTrackNt->fPzGen = particle->Pz();
- printf("Gen: Px %f Py %f Pz %f \n", recTrackNt->fPxGen, recTrackNt->fPyGen, recTrackNt->fPzGen);
- }
- }
- } // for (trackIndex = 0;...
-
- printf(">>> Filling Ntuple of reconstructed tracks\n");
- ntuple->Fill();
-
- return;
-}
-
-void MUONrecoNtuple (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
-{
- //
- cout << "MUON_recoNtuple" << endl;
- cout << "FirstEvent " << FirstEvent << endl;
- cout << "LastEvent " << LastEvent << endl;
- cout << "RecGeantHits " << RecGeantHits << endl;
- cout << "FileName ``" << FileName << "''" << endl;
- cout << "BkgGeantFileName ``" << BkgGeantFileName << "''" << endl;
-// // Dynamically link some shared libs
-// if (gClassTable->GetID("AliRun") < 0) {
-// gROOT->LoadMacro("loadlibs.C");
-// loadlibs();
-// }
-
-
- // Creating Run Loader and openning file containing Hits, Digits and RecPoints
- AliRunLoader * RunLoader = AliRunLoader::Open(FileName,"Event","UPDATE");
- if (RunLoader ==0x0) {
- printf(">>> Error : Error Opening %s file \n",FileName);
- return;
- }
- // Loading AliRun master
- RunLoader->LoadgAlice();
- gAlice = RunLoader->GetAliRun();
- RunLoader->LoadKinematics("READ");
-
- // Loading MUON subsystem
- AliMUON * MUON = (AliMUON *) gAlice->GetDetector("MUON");
- AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
- MUONLoader->LoadHits("READ");
- MUONLoader->LoadRecPoints("READ");
- AliMUONData * muondata = MUON->GetMUONData();
-
- Int_t ievent, nevents;
- nevents = RunLoader->GetNumberOfEvents();
-
- // Initializations
- //AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
-
- //MUON->SetTreeAddress();
- AliMUONEventReconstructor *Reco = new AliMUONEventReconstructor();
-
- Reco->SetRecGeantHits(RecGeantHits);
-
- // The right place for changing AliMUONEventReconstructor parameters
- // with respect to the default ones
-// Reco->SetMaxSigma2Distance(100.0);
-// Reco->SetPrintLevel(20);
- Reco->SetPrintLevel(0);
-// Reco->SetBendingResolution(0.0);
-// Reco->SetNonBendingResolution(0.0);
- cout << "AliMUONEventReconstructor: actual parameters" << endl;
- Reco->Dump();
-// gObjectTable->Print();
- // Loop over events
- if (LastEvent>nevents) LastEvent = nevents;
- for (Int_t event = FirstEvent; event < LastEvent; event++) {
- cout << "Event: " << event << endl;
- RunLoader->GetEvent(event);
- muondata->SetTreeAddress("H,RC");
- // Int_t nparticles = gAlice->GetEvent(event);
- // cout << "nparticles: " << nparticles << endl;
- // prepare background file and/or event if necessary
- if (RecGeantHits == 1) {
- if (event == FirstEvent) Reco->SetBkgGeantFile(BkgGeantFileName);
- if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
- }
- Reco->EventReconstruct();
- // Dump current event
- Reco->EventDump();
- // Fill Ntuple
- AliMUONEventRecNtupleFill(Reco, 0);
-// gObjectTable->Print();
- MUON->ResetRawClusters();
- } // Event loop
- // Write Ntuple
- AliMUONEventRecNtupleFill(Reco, -1);
-}
+++ /dev/null
-// Macro for muon reconstruction
-//
-// This macro is a simplified way to call the macro MUONrecoNtuple.C which produces
-// an Ntuple with reconstructed tracks in output file "MUONtrackReco.root".
-////
-// Arguments:
-// FirstEvent (default 0)
-// LastEvent (default 0)
-// RecGeantHits (1 to reconstruct GEANT hits) (default 0)
-// FileName (for signal) (default "galice.root")
-// BkgGeantFileName (for background),
-// needed only if RecGeantHits = 1 and background to be added
-
-void MUONrecoNtupleRun (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0, Text_t *FileName = "galice.root", Text_t *BkgGeantFileName = "")
-{
-
- TStopwatch timer;
- timer.Start();
-
- gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER -I$ROOTSYS/include") ;
- gROOT->ProcessLine(".x loadlibs.C");
- // force macros compilation (++), load the resulting shared library and execute the macro MUONrecoNtuple.C
- gROOT->ProcessLine(".x MUONrecoNtuple.C++(FirstEvent,LastEvent,RecGeantHits,FileName,BkgGeantFileName) ");
-
- timer.Stop();
- cout << "**** MUONrecoNtupleRun.C , timer" <<endl;
- timer.Print();
-
-}
-
-
-
-
-
+++ /dev/null
-#include <iostream.h>
-
-void MUONrecodisplay(Int_t evNumber=0)
-{
-////////////////////////////////////////////////////////////////////
-// The display pops up a context menu when the right-button is //
-// clicked on the main pad. The following functions are available //
-// * SetDrawHits() - switches on or off Geant hits ; //
-// * CutMomentum() - displays only tracks within Pmin - Pmax //
-// * ListTracks() - prints ID and momentum info. for all //
-// tracks within momentum range Pmin,Pmax; //
-// * Highlight() - shows only one selected reco. track //
-// and its best matching Geant track; //
-// *UnHighlight() - self explaining; //
-// *SetDrawHits() - switch on or off Geant hits //
-////////////////////////////////////////////////////////////////////
- if (evNumber<0) return;
-
-
-//Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-
-// Connect the Root Galice file containing ...
-
- TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!galice_file) galice_file = new TFile("galice.root");
- if (!galice_file) {
- cout << "File galice.root not found\n";
- return;
- }
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)galice_file->Get("gAlice");
- if (gAlice) {
- cout << "AliRun object found on file\n";
- } else {
- cout << "AliRun object not found on file !!!\n";
- return;
- }
- }
-
- AliMUONRecoDisplay *display = new AliMUONRecoDisplay(evNumber);
-
- display->DisableDetector("ITS");
- display->DisableDetector("TPC");
- display->DisableDetector("TOF");
- display->DisableDetector("RICH");
- display->DisableDetector("ZDC");
- display->DisableDetector("CASTOR");
- display->DisableDetector("TRD");
- display->DisableDetector("FMD");
- display->DisableDetector("PHOS");
- display->DisableDetector("PMD");
-
- display->ShowNextEvent(0);
-
-}
+++ /dev/null
- TH1F *mass1 = new TH1F("mass1","Invariant Mass",120,0,12);
- TH1F *mass2 = new TH1F("mass2","Invariant Mass",100,9.0,10.5);
- TH1F *mass3 = new TH1F("mass3","Invariant Mass",100,2.5,3.5);
-void MUONstraggling (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
- static Float_t xmuon, ymuon;
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n Create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-
-// Create some histograms
-
-
- AliMUONChamber* iChamber;
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- for (Int_t nev=0; nev<= evNumber2; nev++) {
- Int_t nparticles = gAlice->GetEvent(nev);
- //cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
-
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-//
-// Loop over tracks
-//
-
- Float_t pups[4]={0,0,0,0};
- for (Int_t track=0; track<ntracks;track++) {
- Int_t Nelt=0;
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
-
-
- if (MUON) {
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->Chamber(); // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- Float_t z = mHit->Z(); // y-pos
- Float_t p = mHit->Momentum();
- Float_t px = mHit->Px();
- Float_t py = mHit->Py();
- Float_t pz = mHit->Pz();
-
- if (nch != 1) continue;
-
- Int_t ipart = mHit->Particle();
- TParticle *Part;
- Int_t ftrack = mHit->Track();
- Part = gAlice->Particle(ftrack);
- Int_t ipart = Part->GetPdgCode();
- TParticle *Mother;
- Float_t px0 = Part->Px();
- Float_t py0 = Part->Py();
- Float_t pz0 = Part->Pz();
- Float_t p0 = Part->P();
-
- Float_t e0=Part->Energy();
-
- if (ipart == kMuonPlus || ipart == kMuonMinus) {
-//
-// Branson Correction
-//
- Float_t zch = z;
- Float_t r = TMath::Sqrt(x*x+y*y);
- Float_t zb;
-
- if (r<26.3611) {
- zb = 466.;
- } else {
- zb = 441.;
- }
-
- Float_t xb = x-(zch-zb)*px/pz;
- Float_t yb = y-(zch-zb)*py/pz;
- Float_t pzB = p * zb/TMath::Sqrt(zb*zb+yb*yb+xb*xb);
- Float_t pxB = pzB * xb/zb;
- Float_t pyB = pzB * yb/zb;
- Float_t theta
- = TMath::ATan2(TMath::Sqrt(pxB*pxB+pyB*pyB), pzB) *
- 180./TMath::Pi();
-
-
-//
-// Energy Correction
-//
-//
-
- Float_t corr=(p+CorrectP(p, theta))/p;
- pups[0]+=p*corr;
- pups[1]+=pxB*corr;
- pups[2]+=pyB*corr;
- pups[3]+=pzB*corr;
- }
- } // hits
- } // if MUON
- } // tracks
- Float_t mass =
- TMath::Sqrt(pups[0]*pups[0]
- -pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]);
- mass1->Fill(mass, 1.);
- mass2->Fill(mass, 1.);
- mass3->Fill(mass, 1.);
- } // event
-//Create a canvas, set the view range, show histograms
- TCanvas *c1 = new TCanvas("c1","Mass Spectra",400,10,600,700);
-
- mass2->SetFillColor(42);
- mass2->SetXTitle("Mass (GeV)");
- mass2->Draw();
-}
-
-
-
-Float_t CorrectP(Float_t p, Float_t theta)
-{
- Float_t c;
-
- if (theta<3.) {
-//W
- if (p<15) {
- c = 2.916+0.0156*p-0.00000866*p*p;
- } else {
- c = 2.998+0.0137*p;
- }
- } else {
-//Cu+Fe+Cc+C
- if (p<15) {
- c = 2.464+0.00877*p-0.0000106*p*p;
- } else {
- c = 2.496+0.00817*p;
- }
- }
-//
-//
- c*=0.95;
-//
- return c;
-}
-
-
-
-
-
-
+++ /dev/null
-void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
- gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
- gSystem->Load("$ROOTSYS/lib/libEG.so");
- gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
- gSystem->Load("libGeant3Dummy.so"); //a dummy version of Geant3
- gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
- gSystem->Load("libgalice.so"); // the standard Alice classes
- }
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
- file->ls();
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-// Create some histograms
-
-
- TH1F *zv = new TH1F("zv","z vertex"
- ,100,468.,503.);
- TH1F *xv = new TH1F("xv","x vertex"
- ,100,0.,50.);
- TH1F *yv = new TH1F("yv","y vertex"
- ,100,0.,50.);
-
- TH1F *ip = new TH1F("ip","ipart"
- ,100,0.,50.);
-
- AliMUONchamber* iChamber;
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- for (Int_t nev=0; nev<= evNumber2; nev++) {
- cout << "nev " << nev <<endl;
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- printf("\n nev %d \n ", nev);
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
-//
-// Loop over tracks
-//
- for (Int_t track=0; track<ntracks;track++) {
-
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (MUON) {
- for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONhit*)MUON->NextHit())
- {
- Int_t nch = mHit->fChamber; // chamber number
- Float_t x = mHit->fX; // x-pos of hit
- Float_t y = mHit->fY; // y-pos
- Float_t z = mHit->fZ; // y-pos
-
- if (nch > 1) continue;
-
- Int_t ipart = mHit->fParticle;
- TClonesArray *fPartArray = gAlice->Particles();
- TParticle *Part;
- Int_t ftrack = mHit->fTrack;
- Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
- ip->Fill((float)ipart,(float) 1);
-
- if (ipart != 3) continue;
- Float_t xv = Part->Vx(); // vertex
- Float_t yv = Part->Vy(); // vertex
- Float_t zv = Part->Vz(); // z vertex
- xv->Fill(xv,(float) 1);
- yv->Fill(yv,(float) 1);
- zv->Fill(zv,(float) 1);
- }
-
- }
- }
- }
-
-
-//Create a canvas, set the view range, show histograms
- TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
- pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(11);
- pad12->SetFillColor(11);
- pad13->SetFillColor(11);
- pad14->SetFillColor(11);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- ip->SetFillColor(42);
- ip->SetXTitle("ipart");
- ip->Draw();
-
- pad12->cd();
- xv->SetFillColor(42);
- xv->SetXTitle("xvert");
- xv->Draw();
-
- pad13->cd();
- yv->SetFillColor(42);
- yv->SetXTitle("yvert");
- yv->Draw();
-
- pad14->cd();
- zv->SetFillColor(42);
- zv->SetXTitle("zvert");
- zv->Draw();
-
-}
+++ /dev/null
-void MUONtestsim(Int_t ncham, Int_t evNumber1=0,Int_t evNumber2=0)
-{
-/////////////////////////////////////////////////////////////////////////
-// This macro is a small example of a ROOT macro
-// illustrating how to read the output of GALICE
-// and do some analysis.
-//
-/////////////////////////////////////////////////////////////////////////
-
-// Dynamically link some shared libs
-
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else {
- printf("\n galice.root found in file list");
- }
- file->ls();
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-// Create some histograms
-
-
- TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
- ,100,0.,5000.);
- TH1F *pcharge1 = new TH1F("pcharge1","Pad Charge distribution"
- ,100,0.,200.);
- TH1F *xresid1 = new TH1F("xresid1","x-Residuals"
- ,100,-10.,10);
- TH1F *yresid1 = new TH1F("yresid1","y-Residuals"
- ,100,-0.2,0.2);
- TH1F *npads1 = new TH1F("npads1" ,"Pads per Hit"
- ,20,-0.5,19.5);
- TH2F *xresys1 = new TH2F("xresys1","x-Residuals systematics"
- ,50,-2.0,2.0,100,-1.0,1.0);
- TH2F *yresys1 = new TH2F("yresys1","y-Residuals systematics"
- ,50,-0.4,0.4,100,-0.1,0.1);
-
- TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
- ,100,0.,5000.);
- TH1F *pcharge2 = new TH1F("pcharge2","Pad Charge distribution"
- ,100,0.,200.);
- TH1F *xresid2 = new TH1F("xresid2","x-Residuals"
- ,100,-1,1);
- TH1F *yresid2 = new TH1F("yresid2","y-Residuals"
- ,100,-10, 10.);
- TH1F *npads2 = new TH1F("npads2" ,"Pads per Hit"
- ,20,-0.5,19.5);
- TH2F *xresys2 = new TH2F("xresys2","x-Residuals systematics"
- ,50,-2.0,2.0,100,-1.0,1.0);
- TH2F *yresys2 = new TH2F("yresys2","y-Residuals systematics"
- ,50,-0.4,0.4,100,-0.1,0.1);
-
-
- AliMUONChamber* iChamber;
-//
-// Loop over events
-//
- Int_t Nh=0;
- Int_t Nh1=0;
- for (Int_t nev=0; nev<= evNumber2; nev++) {
- cout << "nev " << nev <<endl;
- Int_t nparticles = gAlice->GetEvent(nev);
- cout << "nparticles " << nparticles <<endl;
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetDetector("MUON");
- printf("\n track %d %d \n ", nev, MUON);
-
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- Int_t Nc=0;
- Int_t npad[2];
- Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
-//
-// Loop over events
-//
- for (Int_t track=0; track<ntracks;track++) {
-
- gAlice->ResetHits();
- Int_t nbytes += TH->GetEvent(track);
- if (MUON) {
- for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)MUON->NextHit())
- {
- Int_t nch = mHit->fChamber; // chamber number
- Float_t x = mHit->X(); // x-pos of hit
- Float_t y = mHit->Y(); // y-pos
- Float_t z = mHit->Z();
-
-//
-//
- iChamber = & MUON->Chamber(nch-1);
- response=iChamber->ResponseModel();
-//
-//
- if (nch == ncham) {
- for (Int_t i=0; i<2; i++) {
- xmean[i]=0;
- ymean[i]=0;
- Q[i]=0;
- npad[i]=0;
- }
-
- for (AliMUONPadHit* mPad=
- (AliMUONPadHit*)MUON->FirstPad(mHit,MUON->PadHits());
- mPad;
- mPad=(AliMUONPadHit*)MUON->NextPad(MUON->PadHits()))
- {
-
- Int_t nseg = mPad->fCathode; // segmentation
- Int_t nhit = mPad->fHitNumber; // hit number
- Int_t qtot = mPad->fQ; // charge
- Int_t ipx = mPad->fPadX; // pad number on X
- Int_t ipy = mPad->fPadY; // pad number on Y
- Int_t iqpad = mPad->fQpad; // charge per pad
- Int_t secpad = mPad->fRSec; // r-pos of pad
-//
-//
- segmentation=iChamber->SegmentationModel(nseg);
- Int_t ind=nseg-1;
- Float_t xg,yg,zg;
- segmentation->GetPadC(ipx, ipy, xg, yg, zg);
- Int_t sec = segmentation->Sector(ipx, ipy);
- Float_t dpx=segmentation->Dpx(sec);
- Float_t dpy=segmentation->Dpy(sec);
-
- if (nseg == 1) {
- pcharge1->Fill(iqpad,(float) 1);
- ccharge1->Fill(qtot ,(float) 1);
- } else {
- pcharge2->Fill(iqpad,(float) 1);
- ccharge2->Fill(qtot ,(float) 1);
- }
-// Calculate centre of gravity
-//
- if (iqpad == 0 && ind==0) {
- printf("\n attention iqpad 0");
- }
-
- if (iqpad > 0) {
- Float_t xc,yc,zc;
- npad[ind]++;
- segmentation->GetPadC(ipx,ipy,xc,yc,zc);
- xmean[ind]+=Float_t(iqpad*xc);
- ymean[ind]+=Float_t(iqpad*yc);
- Q[ind]+=iqpad;
- }
-
- } //pad hit loop
- for (Int_t i=0; i<2; i++) {
- segmentation = iChamber->SegmentationModel(i+1);
- if (Q[i] >0) {
- xmean[i] = xmean[i]/Q[i];
- xres[i] = xmean[i]-x;
- ymean[i] = ymean[i]/Q[i];
- yres[i] = ymean[i]-y;
-// Systematic Error
-//
- Int_t icx, icy;
- Float_t zonpad;
-
- segmentation->
- GetPadI(x,y,z,icx,icy);
- segmentation->
- GetPadC(icx,icy,xonpad[i],yonpad[i],zonpad);
- xonpad[i]-=x;
- yonpad[i]-=y;
- } // charge not 0
- } // plane loop
- xresid1->Fill(xres[0],(float) 1);
- yresid1->Fill(yres[0],(float) 1);
- npads1->Fill(npad[0],(float) 1);
- if (npad[0] >=2 && Q[0] > 6 ) {
- xresys1->Fill(xonpad[0],xres[0],(float) 1);
- yresys1->Fill(yonpad[0],yres[0],(float) 1);
- }
-
- xresid2->Fill(xres[1],(float) 1);
- yresid2->Fill(yres[1],(float) 1);
- npads2->Fill(npad[1],(float) 1);
- if (npad[1] >=2 && Q[1] > 6) {
- xresys2->Fill(xonpad[1],xres[1],(float) 1);
- yresys2->Fill(yonpad[1],yres[1],(float) 1);
- }
- } // chamber 1
- } // hit loop
- } // if MUON
- } // track loop
- } // event loop
-
-//Create a canvas, set the view range, show histograms
- TCanvas *c1 = new TCanvas("c1","Charge and Residuals (1)",400,10,600,700);
- pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
- pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
- pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
- pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
- pad11->SetFillColor(11);
- pad12->SetFillColor(11);
- pad13->SetFillColor(11);
- pad14->SetFillColor(11);
- pad11->Draw();
- pad12->Draw();
- pad13->Draw();
- pad14->Draw();
-
- pad11->cd();
- ccharge1->SetFillColor(42);
- ccharge1->SetXTitle("ADC units");
- ccharge1->Draw();
-
- pad12->cd();
- pcharge1->SetFillColor(42);
- pcharge1->SetXTitle("ADC units");
- pcharge1->Draw();
-
- pad13->cd();
- xresid1->SetFillColor(42);
- xresid1->Draw();
-
- pad14->cd();
- yresid1->SetFillColor(42);
- yresid1->Draw();
-
- TCanvas *c2 = new TCanvas("c2","Cluster Size (1)",400,10,600,700);
- pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
- pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
- pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
- pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
- pad21->SetFillColor(11);
- pad22->SetFillColor(11);
- pad23->SetFillColor(11);
- pad24->SetFillColor(11);
- pad21->Draw();
- pad22->Draw();
- pad23->Draw();
- pad24->Draw();
-
- pad21->cd();
- npads1->SetFillColor(42);
- npads1->SetXTitle("Cluster Size");
- npads1->Draw();
-
- pad23->cd();
- xresys1->SetXTitle("x on pad");
- xresys1->SetYTitle("x-xcog");
- xresys1->Draw();
-
- pad24->cd();
- yresys1->SetXTitle("y on pad");
- yresys1->SetYTitle("y-ycog");
- yresys1->Draw();
-
- TCanvas *c3 = new TCanvas("c3","Charge and Residuals (2)",400,10,600,700);
- pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
- pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
- pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
- pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
- pad31->SetFillColor(11);
- pad32->SetFillColor(11);
- pad33->SetFillColor(11);
- pad34->SetFillColor(11);
- pad31->Draw();
- pad32->Draw();
- pad33->Draw();
- pad34->Draw();
-
- pad31->cd();
- ccharge2->SetFillColor(42);
- ccharge2->SetXTitle("ADC units");
- ccharge2->Draw();
-
- pad32->cd();
- pcharge2->SetFillColor(42);
- pcharge2->SetXTitle("ADC units");
- pcharge2->Draw();
-
- pad33->cd();
- xresid2->SetFillColor(42);
- xresid2->Draw();
-
- pad34->cd();
- yresid2->SetFillColor(42);
- yresid2->Draw();
-
- TCanvas *c4 = new TCanvas("c4","Cluster Size (2)",400,10,600,700);
- pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
- pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
- pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
- pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
- pad41->SetFillColor(11);
- pad42->SetFillColor(11);
- pad43->SetFillColor(11);
- pad44->SetFillColor(11);
- pad41->Draw();
- pad42->Draw();
- pad43->Draw();
- pad44->Draw();
-
- pad41->cd();
- npads2->SetFillColor(42);
- npads2->SetXTitle("Cluster Size");
- npads2->Draw();
-
- pad43->cd();
- xresys2->SetXTitle("x on pad");
- xresys2->SetYTitle("x-xcog");
- xresys2->Draw();
-
- pad44->cd();
- yresys2->SetXTitle("y on pad");
- yresys2->SetYTitle("y-ycog");
- yresys2->Draw();
-
-}
+++ /dev/null
-void MUONtrackRecoModel() {
-//////////////////////////////////////////////////////////
-// This file has been automatically generated
-// (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
-// from TTree MUONtrackReco/MUONtrackReco
-// found on file: MUONtrackReco.root
-//////////////////////////////////////////////////////////
-
-
-//Reset ROOT and connect tree file
- gROOT->Reset();
- TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
- if (!f) {
- f = new TFile("MUONtrackReco.root");
- }
- TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
-
-//Declaration of leaves types
- Int_t fEvent;
- UInt_t fUniqueID;
- UInt_t fBits;
- Int_t Tracks_;
- Int_t Tracks_fCharge[5];
- Float_t Tracks_fPxRec[5];
- Float_t Tracks_fPyRec[5];
- Float_t Tracks_fPzRec[5];
- Float_t Tracks_fZRec[5];
- Float_t Tracks_fZRec1[5];
- Int_t Tracks_fNHits[5];
- Float_t Tracks_fChi2[5];
- Float_t Tracks_fPxGen[5];
- Float_t Tracks_fPyGen[5];
- Float_t Tracks_fPzGen[5];
- UInt_t Tracks_fUniqueID[5];
- UInt_t Tracks_fBits[5];
-
-//Set branch addresses
- //MUONtrackReco->SetBranchAddress("Header",&Header);
- MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
- MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
- MUONtrackReco->SetBranchAddress("fBits",&fBits);
- MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
- MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
- MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
- MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
- MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
- MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
- MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
- MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
- MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
- MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
- MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
- MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
- MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
- MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
-
-// This is the loop skeleton
-// To read only selected branches, Insert statements like:
-// MUONtrackReco->SetBranchStatus("*",0); // disable all branches
-// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
-
- Int_t nentries = MUONtrackReco->GetEntries();
-
- Int_t nbytes = 0;
-// for (Int_t i=0; i<nentries;i++) {
-// nbytes += MUONtrackReco->GetEntry(i);
-// }
-}
+++ /dev/null
-void MUONtracking (Int_t evNumber1=0, Int_t evNumber2=99, Int_t idres=116, Int_t ireadgeant=1, Int_t ibgr=1, Int_t nev_bgd=10)
-{
- //////////////////////////////////////
- // //
- // ROOT macro for ALICE Dimuon Arm: //
- // Track reconstruction //
- // //
- //////////////////////////////////////
- //
- // Reconstructs tracks from events in the ROOT file "galice.root".
- // Track reconstruction is performed (argument "ireadgeant")
- // either directly from GEANT hits (tree TH),
- // or from raw clusters (tree TR) constructed from digits (tree TD).
- // Eventually (argument "ibgr"), background GEANT hits
- // are also taken into account from the ROOT file "galice_bgr.root".
- //
- // Arguments:
- // evNumber1 = first event number to act on in file "galice.root"
- // evNumber2 = last event number to act on in file "galice.root"
- // idres : used for statistics
- // = 116 for Upsilon
- // = 114 for J/Psi
- // ireadgeant = 1 to reconstruct tracks directly from GEANT hits
- // = 0 to reconstruct tracks from raw clusters
- // ibg : used only if "ireadgeant" = 1
- // = 0 if no background GEANT hits to be taken into account;
- // = 1 if background GEANT hits
- // to be taken into account in file "galice_bgr.root"
- // used only if "ireadgeant" = 1
- // nev_bgd : used only if "ireadgeant" = 1 and "ibg" = 1
- // = number of events in the background file "galice_bgr.root";
- // successive signal events are mixed
- // with different background events in file "galice_bgr.root",
- // starting from event number 0,
- // incrementing it by 1 till event number ("nev_bgd" - 1),
- // continuing with event number 0 and so on.
- // Strictly speaking, "nev_bgd" can be smaller than
- // the number of events in the background file,
- // in which case one will only use
- // the first "nev_bgd" events of the background file.
- // But it SHOULD NOT BE LARGER THAN
- // THE ACTUAL NUMBER OF EVENTS IN THE BACKGROUND FILE.
- //
- // Input file(s):
- // "galice.root" for GEANT hits or raw clusters
- // "galice_bgr.root" for background GEANT hits
- // (used only if "ireadgeant" = 1 and "ibg" = 1)
- //
- // Output file:
- // "reconst.root" for ROOT ntuples
- //
- //__________________________________________________________________________
-
-// Dynamically link some shared libs
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-// Connect the Root Galice file containing Geometry, Kine and Hits
-
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
-
- if (!file) {
- printf("\n Creating galice.root \n");
- file = new TFile("galice.root");
- } else { printf("\n galice.root found in file list");
- }
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)(file->Get("gAlice"));
- if (gAlice) printf("AliRun object found on file\n");
- if (!gAlice) {
- printf("\n create new gAlice object");
- gAlice = new AliRun("gAlice","Alice test program");
- }
- }
-
-// seff = efficiency per chamber (ireadgeant=1)
- Double_t seff = 1;
- // Double_t seff = 1.;
-// sb0 = magn. field in dipole, sbl3 = magn. field in L3
-// necessary for trackfinding only.
- Double_t sb0 = 0.7;
- Double_t sbl3 = 0.2;
-
-// ifit = 0 trackfinding only
-// ifit = 1 trackfinding + fit
- Int_t ifit = 1;
-// idebug = 0,1,2 print level for reco_muon.F
- Int_t idebug = 1;
-
- AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
- AliMUONTrackReconstructor *Reconstruction = new AliMUONTrackReconstructor();
- Int_t nparticles = gAlice->GetEvent(evNumber1);
- if (nparticles <= 0) return;
-
- Reconstruction->Init(seff,sb0,sbl3);
-
-// Loop over events
-//
- Int_t inev_bgd=0;
-
- for (Int_t nev= evNumber1; nev<= evNumber2; nev++)
- {
- printf("nev=%d\n",nev);
- if (nev != evNumber1) Int_t nparticles = gAlice->GetEvent(nev);
- if (nev < evNumber1) continue;
- if (nparticles <= 0) return;
- Reconstruction->FinishEvent();
- if (ireadgeant==1 && ibgr==1) {
- if (inev_bgd==nev_bgd) inev_bgd=0;
- Reconstruction->Reconst(ifit,idebug,inev_bgd,nev,idres,ireadgeant,"Add","galice_bgr.root");
- Reconstruction->FinishEvent();
- inev_bgd++;
- }
- else {
- // printf("ireadgeant=%d\n",ireadgeant);
- Reconstruction->Reconst(ifit,idebug,inev_bgd,nev,idres,ireadgeant,"rien1","rien2");
- Reconstruction->FinishEvent();
-
- }
-
- } // event loop
-
- Reconstruction->Close();
-}
-
-
-
+++ /dev/null
-#include <iostream.h>
-
-
-//***************************************************************************
-
-void ReadTree()
-{
-// This is a basic example on how the tree of reconstructed events
-// should be accessed
-//Dynamically link some shared libs
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C");
- loadlibs();
- }
-
-
-// Connect the Root Galice file containing ...
- TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("tree_reco.root");
- if (!file) file = new TFile("tree_reco.root");
- if (!file) {
- cout << "File tree_reco.root not found\n";
- return;
- }
-
- TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
- if (!galice_file) galice_file = new TFile("galice.root");
- if (!galice_file) {
- cout << "File galice.root not found\n";
- return;
- }
-
-// Get AliRun object from file or create it if not on file
- if (!gAlice) {
- gAlice = (AliRun*)galice_file->Get("gAlice");
- if (gAlice) {
- cout << "AliRun object found on file\n";
- } else {
- cout << "AliRun object not found on file !!!\n";
- return;
- }
- }
-
- TTree *tree = (TTree*)file->Get("TreeRecoEvent");
- TBranch *branch = tree->GetBranch("Event");
- AliMUONRecoEvent *event = 0;
- branch->SetAddress(&event);
- Int_t nentries = (Int_t)tree->GetEntries(); // no. of reco events on file
-
- for (Int_t evNumber=0; evNumber<nentries; evNumber++) {
- tree->GetEntry(evNumber);
- cout << "Event : " << event->GetNoEvent() << " with" << event->GetNoTracks() << " tracks\n";
-// print reconstr. event information
-// event->EventInfo();
- }
-}
+++ /dev/null
-void masse() {
-//////////////////////////////////////////////////////////
-// This file has been automatically generated
-// (Wed Jun 16 15:35:08 1999 by ROOT version 2.21/07)
-// from TTree ntuple/Reconst ntuple
-// found on file: reconst.root
-//////////////////////////////////////////////////////////
-
-//Reset ROOT and connect tree file
-gROOT->Reset();
-
-TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("reconst.root");
-if (!f) {
- f = new TFile("reconst.root");
-}
-
-TTree *ntuple = (TTree*)gDirectory->Get("ntuple");
-
-TPostScript ps("mass.ps",112);
-//ps->Open("mass.ps");
-//Declaration of leaves types
-Int_t ntrack=500;
-
-Int_t ievr;
-Int_t ntrackr;
-Int_t istatr[ntrack];
-Int_t isignr[ntrack];
-Float_t pxr[ntrack];
-Float_t pyr[ntrack];
-Float_t pzr[ntrack];
-Float_t zvr[ntrack];
-Float_t chi2r[ntrack];
-
-//Set branch addresses
-ntuple->SetBranchAddress("ievr",&ievr);
-ntuple->SetBranchAddress("ntrackr",&ntrackr);
-ntuple->SetBranchAddress("istatr",&istatr);
-ntuple->SetBranchAddress("isignr",&isignr);
-ntuple->SetBranchAddress("pxr",&pxr);
-ntuple->SetBranchAddress("pyr",&pyr);
-ntuple->SetBranchAddress("pzr",&pzr);
-ntuple->SetBranchAddress("zvr",&zvr);
-ntuple->SetBranchAddress("chi2r",&chi2r);
-
-// This is the loop skeleton
-// To read only selected branches, Insert statements like:
-// ntuple->SetBranchStatus("*",0); // disable all branches
-// ntuple->SetBranchStatus("branchname",1); // activate branchname
-
-TH1F *h1 = new TH1F("h1","mass",240,0.,12.);
-TH1F *h2 = new TH1F("h2","mass good muon",240,0.,12.);
-// Upsilon mass cut
-TH1F *h5 = new TH1F("h5"," Mass Upsilon",70,8.,11.);
-// J/psi mass cut
-//TH1F *h5 = new TH1F("h5"," Mass J/psi",120,1.,5.);
-TH1F *h3 = new TH1F("h3","Pt",100,0.,20.);
-TH1F *h4 = new TH1F("h4","Pt bon muon",100,0.,20.);
-
-gStyle->SetOptFit();
-//gStyle->SetOptStat(0000);
-gStyle->SetOptStat(10001110);
-
-Float_t mass = 1.;
-Float_t amassmu = 0.10566;
-Float_t Chi2Cut = 100.;
-Float_t PtCut = 0.;
-//Float_t MassMin = 8.96;
-//Float_t MassMax = 9.96;
-// Upsilon
-Float_t MassMin = 8.5;
-Float_t MassMax = 10.5;
-// J/psi
-//Float_t MassMin = 2.5;
-//Float_t MassMax = 3.8;
-
-//Float_t Chi2Cut = 999.;
-
-Int_t nentries = ntuple->GetEntries();
-
-Int_t nbytes = 0;
-Int_t nbtrack[50]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
-
-Int_t compt = 0;
-Int_t Nres = 0;
-
-// Good muons
-for (Int_t iev=0; iev<nentries; iev++) {
- nbytes += ntuple->GetEvent(iev);
-// ntuple->SetBranchStatus("main",1);
- ntuple->SetBranchStatus("*",1);
-
- //printf("\n");
- // printf("ntrackr=%d\n",ntrackr);
-
- if(ntrackr < 50) nbtrack[ntrackr]=nbtrack[ntrackr]+1 ;
-
-
- for (Int_t i=0; i<(ntrackr-1); i++) {
-
- Float_t px1 = pxr[i] ;
- Float_t py1 = pyr[i] ;
- Float_t pz1 = pzr[i] ;
- Float_t Pt1 = TMath::Sqrt(px1*px1+py1*py1) ;
-
-// if (chi2r[i] < Chi2Cut && istatr[i] == 1 && Pt1 > PtCut) {
- h4->Fill(Pt1);
-
- for (Int_t j=i+1; j<ntrackr; j++) {
- Float_t px2 = pxr[j] ;
- Float_t py2 = pyr[j] ;
- Float_t pz2 = pzr[j] ;
- Float_t Pt2 = TMath::Sqrt(px2*px2+py2*py2) ;
- if (chi2r[j]<Chi2Cut && istatr[j] == 1 && Pt2 > PtCut ) {
- if (isignr[i]!=isignr[j]) {
- Float_t p12 = px1*px1+py1*py1+pz1*pz1;
- Float_t p22 = px2*px2+py2*py2+pz2*pz2;
- Float_t e1 = TMath::Sqrt(amassmu*amassmu+p12);
- Float_t e2 = TMath::Sqrt(amassmu*amassmu+p22);
- Float_t amass = TMath::Sqrt((e1+e2)*(e1+e2)-(px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
-// printf("amass=%f\n",amass);
- h2->Fill(amass);
- if (amass<11. && amass>8.) compt++;
- }
- }
- }
- }
-
-
-
-
- // All tracks
- for (Int_t i=0; i<(ntrackr-1); i++) {
- Float_t px1 = pxr[i] ;
- Float_t py1 = pyr[i] ;
- Float_t pz1 = pzr[i] ;
- Float_t Pt1 = TMath::Sqrt(px1*px1+py1*py1) ;
- if (chi2r[i] < Chi2Cut && Pt1 > PtCut) {
- h3->Fill(Pt1);
- for (Int_t j=i+1; j<ntrackr; j++) {
- Float_t px2 = pxr[j] ;
- Float_t py2 = pyr[j] ;
- Float_t pz2 = pzr[j] ;
- Float_t Pt2 = TMath::Sqrt(px2*px2+py2*py2) ;
- if (chi2r[j]<Chi2Cut && Pt2 > PtCut) {
- if (isignr[i]!=isignr[j]) {
- Float_t p12 = px1*px1+py1*py1+pz1*pz1;
- Float_t p22 = px2*px2+py2*py2+pz2*pz2;
- Float_t e1 = TMath::Sqrt(amassmu*amassmu+p12);
- Float_t e2 = TMath::Sqrt(amassmu*amassmu+p22);
- Float_t amass = TMath::Sqrt((e1+e2)*(e1+e2)-(px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
- // printf("amass=%f\n",amass);
- h1->Fill(amass);
- if (amass > MassMin && amass < MassMax) {
- // if (ntrackr==2) h5->Fill(amass);
- h5->Fill(amass);
-// if (amass < 9.715 && amass > 9.205)
- if (amass < 9.766 && amass > 9.154)
- {
- printf("amass=%f\n",amass);
- Nres++;
- }
- }
- }
- }
- }
- }
- }
-}
-
-// printf("compt= %d\n",compt);
-
-for (Int_t i = 0 ; i<50 ; i++)
-{
-
- printf(" nb of events= %d with %d tracks \n",nbtrack[i], i) ;
-
-}
-g1= new TF1("g1","gaus",9.30,9.75) ;
-//g1= new TF1("g1","gaus",8.7,10.2) ;
-//g1= new TF1("g1","gaus",2.95,3.25) ;
-h5->Fit("g1","R0") ;
-h5->SetXTitle("Mass (GeV/c^2!)");
-h5->Draw();
-printf("Nombre de resonances : %d\n",Nres);
-
-ps->Close();
-}
-
-
-
-
-
-
-
-
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff