-void MUONdigit (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1)
+#include "iostream.h"
+
+void MUONdigit (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal =25)
{
/////////////////////////////////////////////////////////////////////////
-// //
-// This macros converts the pad-hit information into digits //
-// //
+// 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
+ 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");
file->ls();
-// file->Map();
- printf ("\n File loaded !");
+ printf ("I'm after Map \n");
// Get AliRun object from file or create it if not on file
if (!gAlice) {
gAlice = (AliRun*)file->Get("gAlice");
- if (gAlice) printf("\n AliRun object found on file\n");
+ if (gAlice) printf("AliRun object found on file\n");
if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
}
- printf ("\n gAlice created !\n");
-
- AliMUON *MUON = gAlice->GetDetector("MUON");
-
- AliMUONchamber* iChamber;
- AliMUONsegmentation* segmentation;
-
- Int_t Npx[10];
- Int_t Npy[10];
-
- Int_t nxmax=1026;
- Int_t nymax=1026;
- AliMUONlist *elem[10][1026][1026];
- TObjArray *obj=new TObjArray;
-
- Int_t trk[500];
- Int_t chtrk[500];
-
-
- Int_t digits[3];
- TVector *trinfo;
- trinfo=new TVector(2);
-
+ printf ("I'm after gAlice \n");
+
+ AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
//
-// Loop over events
+// Event Loop
//
-
- Int_t Nh=0;
- Int_t Nh1=0;
+ Int_t nbgr_ev=0;
+
for (int nev=0; nev<= evNumber2; nev++) {
Int_t nparticles = gAlice->GetEvent(nev);
cout << "nev " <<nev<<endl;
if (nev < evNumber1) continue;
if (nparticles <= 0) return;
- TTree *TH = gAlice->TreeH();
- Int_t ntracks = TH->GetEntries();
- Int_t Nc=0;
-
- Int_t counter=0;
- //
- // loop over cathodes
- for (int icat=0; icat<nCathode; icat++) {
- //
- // initialize the array of pointers
- printf("\n Initialize array of pointers \n");
- for (int i=0; i<10; i++) {
- for (int j=0; j<nymax; j++) {
- for (int k=0; k<nxmax; k++) {
- elem[i][j][k]=0;
- }
- }
- }
- printf("\n Finished !\n");
-
-//
-// Loop over tracks
-//
- printf("\n ntracks: %d",ntracks);
-
- 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
-
-
- iChamber = &(MUON->Chamber(nch-1));
- response=iChamber->GetResponseModel();
-
- Int_t nsec=iChamber->Nsec();
- Int_t rmin = (Int_t)iChamber->frMin;
- Int_t rmax = (Int_t)iChamber->frMax;
- if (nch > 10) continue;
-
-//
-//
- for (AliMUONcluster* mPad=
- (AliMUONcluster*)MUON->FirstPad(mHit);
- mPad;
- mPad=(AliMUONcluster*)MUON->NextPad())
- {
- Int_t cathode = mPad->fCathode; // cathode number
- Int_t nhit = mPad->fHitNumber; // hit number
- Int_t qtot = mPad->fQ; // charge
- Int_t ipx = mPad->fPadX; // pad X
- Int_t ipy = mPad->fPadY; // pad Y
- Int_t iqpad = mPad->fQpad; // charge per pad
- Int_t rsec = mPad->fRSec; // sector#
-//
-//
- if (cathode != (icat+1)) continue;
- segmentation=iChamber->
- GetSegmentationModel(cathode);
- Int_t Npx[nch-1] = segmentation->Npx();
- Int_t Npy[nch-1] = segmentation->Npy();
-
-
- Int_t npx=Npx[nch-1];
- Int_t npy=Npy[nch-1];
- Float_t thex, they;
- segmentation->GetPadCxy(ipx,ipy,thex,they);
-
- Float_t rpad=TMath::Sqrt(thex*thex+they*they);
-
-
- // check boundaries
- if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
-
- // fill the info array
+ Int_t nbgr_ev=Int_t(nev/nsignal);
+// printf("nbgr_ev %d\n",nbgr_ev);
+ //if (MUON) MUON->Digitise(nev,nbgr_ev,"Add"," ","galice_bgr.root");
+ if (MUON) MUON->Digitise(nev,nbgr_ev,"rien"," ","galice_bgr.root");
+// char hname[30];
+// sprintf(hname,"TreeD%d",nev);
+// file->ls();
+ } // event loop
+ file->Close();
+}
- trinfo(0)=(Float_t)track;
- trinfo(1)=(Float_t)iqpad;
- digits[0]=ipx;
- digits[1]=ipy;
- digits[2]=iqpad;
- // build the list of fired pads and update the info
-
- AliMUONlist *pdigit=elem[nch-1][ipy+npy][ipx+npx];
- if (pdigit==0) {
- obj->AddAtAndExpand(
- new AliMUONlist(rpad,digits),counter);
- counter++;
- Int_t last=obj->GetLast();
- pdigit=(AliMUONlist*)obj->At(last);
- elem[nch-1][ipy+npy][ipx+npx]=pdigit;
- // list of tracks
- TObjArray *trlist=
- (TObjArray*)pdigit->TrackList();
- trlist->Add(trinfo);
- } else {
- // update charge
- (*pdigit).fSignal+=iqpad;
- // update list of tracks
- TObjArray* trlist=(TObjArray*)pdigit
- ->TrackList();
- Int_t nptracks=trlist->GetEntriesFast();
-
- for (Int_t tr=0;tr<nptracks;tr++) {
- TVector *ptrk=(TVector*)trlist->At(tr);
-// TVector &vtrk = *ptrk;
- Int_t entries=ptrk->GetNrows();
- trk[tr]=Int_t(ptrk(0));
- chtrk[tr]=Int_t(ptrk(1));
- if (trk[tr]==track) {
- chtrk[tr]+=iqpad;
- trlist->RemoveAt(tr);
- trinfo(0)=trk[tr];
- trinfo(1)=chtrk[tr];
- trlist->AddAt(trinfo,tr);
- } else {
- trlist->Add(trinfo);
- }
- } // end loop over list of tracks for one pad
- } // end if pdigit
- } //end loop over clust
- } // hit loop
- } // if MUON
- } // track loop
- Int_t tracks[10];
- Int_t charges[10];
- cout<<"Start filling digits now !"<<endl;
- const Int_t zero_supm = 6;
- const Int_t adc_satm = 1024;
- Int_t nd=0;
-
- // start filling the digits
- for (Int_t id=0; id<10; id++) {
- Int_t nd=0;
- for (Int_t iy=0;iy<nymax;iy++) {
- for (Int_t ix=0;ix<nxmax;ix++) {
- AliMUONlist *address=elem[id][iy][ix];
- if (address==0) continue;
- // do zero-suppression and signal truncation
- Int_t q=address->fSignal;
- //
- if ( q <= zero_supm ) continue;
- if ( q > adc_satm) q=adc_satm;
- digits[0]=address->fPadX;
- digits[1]=address->fPadY;
- digits[2]=address->fSignal;
- TObjArray* trlist=(TObjArray*)address->TrackList();
- Int_t nptracks=trlist->GetEntriesFast();
- // this was changed to accomodate the real number of tracks
- if (nptracks > 10) {
- cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
- nptracks=10;
- }
- for (Int_t tr=0;tr<nptracks;tr++) {
- TVector *pp=(TVector*)trlist->At(tr);
- Int_t entries=pp->GetNrows();
- TVector &v2 = *pp;
- tracks[tr]=(Int_t)v2(0);
- charges[tr]=(Int_t)v2(1);
- } //end loop over list of tracks for one pad
- // fill digits
- MUON->AddDigits(id,tracks,charges,digits);
- nd++;
- } // end loop over ix
- } // end loop over iy
- cout<<"I'm out of the loops over pads"<<endl;
- cout<<" id nd "<<id<<" "<<nd<<endl;
- } //loop over chambers
- cout<<"I'm out of the loops for digitisation"<<endl;
- gAlice->TreeD()->Fill();
- TTree *TD=gAlice->TreeD();
- int ndig=TD->GetEntries();
- cout<<"number of digits "<<ndig<<endl;
- TClonesArray *fDch;
- for (int i=0;i<10;i++) {
- fDch= MUON->DigitsAddress(i);
- int ndig=fDch->GetEntriesFast();
- printf (" i, ndig %d %d \n",i,ndig);
- }
- MUON->ResetDigits();
- obj->Clear();
- } //end loop over cathodes
- char hname[30];
- sprintf(hname,"TreeD%d",nev);
- gAlice->TreeD()->Write(hname);
-
- file->ls();
-// file->Map();
- } // event loop
- obj->Delete();
- delete obj;
- file->Close();
- cout<<"END digitisation "<<endl;
-
-}