- // keep galice.root for signal and name differently the file for
- // background when add! otherwise the track info for signal will be lost !
-
- static Bool_t first=kTRUE;
- static TFile *file;
- char *addBackground = strstr(option,"Add");
-
-
- AliMUONChamber* iChamber;
- AliSegmentation* segmentation;
-
-
- Int_t trk[50];
- Int_t chtrk[50];
- TObjArray *list=new TObjArray;
- static TClonesArray *pAddress=0;
- if(!pAddress) pAddress=new TClonesArray("TVector",1000);
- Int_t digits[5];
-
- AliMUON *pMUON = (AliMUON *) gAlice->GetModule("MUON");
- AliHitMap** hitMap= new AliHitMap* [AliMUONConstants::NCh()];
- for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {hitMap[i]=0;}
- if (addBackground ) {
- if(first) {
- fFileName=filename;
- cout<<"filename"<<fFileName<<endl;
- file=new TFile(fFileName);
- cout<<"I have opened "<<fFileName<<" file "<<endl;
- fHits2 = new TClonesArray("AliMUONHit",1000 );
- fPadHits2 = new TClonesArray("AliMUONPadHit",10000);
- }
- first=kFALSE;
- file->cd();
- //file->ls();
- // Get Hits Tree header from file
- if(fHits2) fHits2->Clear();
- if(fPadHits2) fPadHits2->Clear();
- if(fTrH1) delete fTrH1;
- fTrH1=0;
-
- char treeName[20];
- sprintf(treeName,"TreeH%d",bgrEvent);
- fTrH1 = (TTree*)gDirectory->Get(treeName);
- //printf("fTrH1 %p of treename %s for event %d \n",fTrH1,treeName,bgrEvent);
-
- if (!fTrH1) {
- printf("ERROR: cannot find Hits Tree for event:%d\n",bgrEvent);
- }
- // Set branch addresses
- TBranch *branch;
- char branchname[20];
- sprintf(branchname,"%s",GetName());
- if (fTrH1 && fHits2) {
- branch = fTrH1->GetBranch(branchname);
- if (branch) branch->SetAddress(&fHits2);
- }
- if (fTrH1 && fPadHits2) {
- branch = fTrH1->GetBranch("MUONCluster");
- if (branch) branch->SetAddress(&fPadHits2);
- }
-// test
- //Int_t ntracks1 =(Int_t)fTrH1->GetEntries();
- //printf("background - ntracks1 - %d\n",ntracks1);
- }
- //
- // loop over cathodes
- //
- AliHitMap* hm;
- Int_t countadr=0;
- for (int icat=0; icat<2; icat++) {
- Int_t counter=0;
- Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
- for (Int_t i =0; i<AliMUONConstants::NCh(); i++) {
- iChamber=(AliMUONChamber*) (*fChambers)[i];
- if (iChamber->Nsec()==1 && icat==1) {
- continue;
- } else {
- segmentation=iChamber->SegmentationModel(icat+1);
- }
- hitMap[i] = new AliMUONHitMapA1(segmentation, list);
- nmuon[i]=0;
- }
- //printf("Start loop over tracks \n");
-//
-// Loop over tracks
-//
-
- TTree *treeH = gAlice->TreeH();
- Int_t ntracks =(Int_t) treeH->GetEntries();
- Int_t jj;
- Int_t nmaxmuon = 5;
-
- Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
- for (jj=0; jj<AliMUONConstants::NCh(); jj++)
- xhit[jj] = new Float_t[nmaxmuon];
- Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
- for (jj=0; jj<AliMUONConstants::NCh(); jj++)
- yhit[jj] = new Float_t[nmaxmuon];
-
- for (Int_t track=0; track<ntracks; track++) {
- gAlice->ResetHits();
- treeH->GetEvent(track);
-//
-// Loop over hits
- for(AliMUONHit* mHit=(AliMUONHit*)pMUON->FirstHit(-1);
- mHit;
- mHit=(AliMUONHit*)pMUON->NextHit())
- {
- Int_t nch = mHit->fChamber-1; // chamber number
- if (nch > AliMUONConstants::NCh()-1) continue;
-// if (nch > 9) continue;
- iChamber = &(pMUON->Chamber(nch));
- // new 17.07.99
- if (addBackground) {
-
- if (mHit->fParticle == kMuonPlus
- || mHit->fParticle == kMuonMinus) {
- // mark muon hits
- if (nmuon[nch] < nmaxmuon) {
- xhit[nch][nmuon[nch]]=mHit->X();
- yhit[nch][nmuon[nch]]=mHit->Y();
- nmuon[nch]++;
- }
- // ignore muon if too many compared to nmaxmuon
- else printf("AliMUON::Digitize: nmuon %d ==> ignored\n",nmuon[nch]);
- }
- }
-
-
-
-
-//
-// Loop over pad hits
- for (AliMUONPadHit* mPad=
- (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
- mPad;
- mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits))
- {
- Int_t cathode = mPad->fCathode; // cathode number
- Int_t ipx = mPad->fPadX; // pad number on X
- Int_t ipy = mPad->fPadY; // pad number on Y
- Int_t iqpad = Int_t(mPad->fQpad);// charge per pad
-// printf("\n Pad: %d %d %d %d", ipx, ipy, cathode,nch);
-//
-//
-
- if (cathode != (icat+1)) continue;
- // fill the info array
-// Float_t thex, they, thez;
- segmentation=iChamber->SegmentationModel(cathode);
-// segmentation->GetPadC(ipx,ipy,thex,they,thez);
-// Float_t rpad=TMath::Sqrt(thex*thex+they*they);
-// if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
-
- new((*pAddress)[countadr++]) TVector(2);
- TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
- trinfo(0)=(Float_t)track;
- trinfo(1)=(Float_t)iqpad;
-
- digits[0]=ipx;
- digits[1]=ipy;
- digits[2]=iqpad;
- digits[3]=iqpad;
- if (mHit->fParticle == kMuonPlus ||
- mHit->fParticle == kMuonMinus) {
- digits[4]=mPad->fHitNumber;
- } else digits[4]=-1;
-
- AliMUONTransientDigit* pdigit;
- // build the list of fired pads and update the info
- if (!hitMap[nch]->TestHit(ipx, ipy)) {
-
- list->AddAtAndExpand(
- new AliMUONTransientDigit(nch,digits),counter);
-
- hitMap[nch]->SetHit(ipx, ipy, counter);
- counter++;
- pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
- // list of tracks
- TObjArray *trlist=(TObjArray*)pdigit->TrackList();
- trlist->Add(&trinfo);
- } else {
- pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
- // update charge
- (*pdigit).fSignal+=iqpad;
- (*pdigit).fPhysics+=iqpad;
- // update list of tracks
- TObjArray* trlist=(TObjArray*)pdigit->TrackList();
- Int_t lastEntry=trlist->GetLast();
- TVector *pTrack=(TVector*)trlist->At(lastEntry);
- TVector &ptrk=*pTrack;
- Int_t lastTrack = Int_t(ptrk(0));
- Int_t lastCharge = Int_t(ptrk(1));
- if (lastTrack==track) {
- lastCharge+=iqpad;
- trlist->RemoveAt(lastEntry);
- trinfo(0)=lastTrack;
- trinfo(1)=lastCharge;
- trlist->AddAt(&trinfo,lastEntry);
- } else {
- trlist->Add(&trinfo);
- }
- // check the track list
- Int_t nptracks=trlist->GetEntriesFast();
- if (nptracks > 2) {
- for (Int_t tr=0;tr<nptracks;tr++) {
- TVector *ppTrack=(TVector*)trlist->At(tr);
- TVector &pptrk=*ppTrack;
- trk[tr] = Int_t(pptrk(0));
- chtrk[tr] = Int_t(pptrk(1));
- }
- } // end if nptracks
- } // end if pdigit
- } //end loop over clusters
- } // hit loop
- } // track loop
-
- // open the file with background
-
- if (addBackground) {
- ntracks =(Int_t)fTrH1->GetEntries();
-//
-// Loop over tracks
-//
- for (Int_t track=0; track<ntracks; track++) {
-
- if (fHits2) fHits2->Clear();
- if (fPadHits2) fPadHits2->Clear();
-
- fTrH1->GetEvent(track);
-//
-// Loop over hits
- AliMUONHit* mHit;
- for(int i=0;i<fHits2->GetEntriesFast();++i)
- {
- mHit=(AliMUONHit*) (*fHits2)[i];
- Int_t nch = mHit->fChamber-1; // chamber number
- if (nch >9) continue;
- iChamber = &(pMUON->Chamber(nch));
-// Int_t rmin = (Int_t)iChamber->RInner();
-// Int_t rmax = (Int_t)iChamber->ROuter();
- Float_t xbgr=mHit->X();
- Float_t ybgr=mHit->Y();
- Bool_t cond=kFALSE;
-
- for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
- Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
- +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
- if (dist<100) cond=kTRUE;
- }
- if (!cond) continue;
-
-//
-// Loop over pad hits
- for (AliMUONPadHit* mPad=
- (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits2);
- mPad;
- mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits2))
- {
- // mPad = (AliMUONPadHit*) (*fPadHits2)[j];
- Int_t cathode = mPad->fCathode; // cathode number
- Int_t ipx = mPad->fPadX; // pad number on X
- Int_t ipy = mPad->fPadY; // pad number on Y
- Int_t iqpad = Int_t(mPad->fQpad);// charge per pad
-
- if (cathode != (icat+1)) continue;
-// printf("\n Pad: %d %d %d", ipx, ipy, cathode);
-
-// Float_t thex, they, thez;
-// segmentation=iChamber->SegmentationModel(cathode);
-// segmentation->GetPadC(ipx,ipy,thex,they,thez);
-// Float_t rpad=TMath::Sqrt(thex*thex+they*they);
-// if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
- new((*pAddress)[countadr++]) TVector(2);
- TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
- trinfo(0)=-1; // tag background
- trinfo(1)=-1;
-
- digits[0]=ipx;
- digits[1]=ipy;
- digits[2]=iqpad;
- digits[3]=0;
- digits[4]=-1;
-
- AliMUONTransientDigit* pdigit;
- // build the list of fired pads and update the info
- if (!hitMap[nch]->TestHit(ipx, ipy)) {
- list->AddAtAndExpand(new AliMUONTransientDigit(nch,digits),counter);
-
- hitMap[nch]->SetHit(ipx, ipy, counter);
- counter++;
-
- pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
- // list of tracks
- TObjArray *trlist=(TObjArray*)pdigit->
- TrackList();
- trlist->Add(&trinfo);
- } else {
- pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
- // update charge
- (*pdigit).fSignal+=iqpad;
-
- // update list of tracks
- TObjArray* trlist=(TObjArray*)pdigit->
- TrackList();
- Int_t lastEntry=trlist->GetLast();
- TVector *pTrack=(TVector*)trlist->
- At(lastEntry);
- TVector &ptrk=*pTrack;
- Int_t lastTrack=Int_t(ptrk(0));
- if (lastTrack==-1) {
- continue;
- } else {
- trlist->Add(&trinfo);
- }
- // check the track list
- Int_t nptracks=trlist->GetEntriesFast();
- if (nptracks > 0) {
- for (Int_t tr=0;tr<nptracks;tr++) {
- TVector *ppTrack=(TVector*)trlist->At(tr);
- TVector &pptrk=*ppTrack;
- trk[tr]=Int_t(pptrk(0));
- chtrk[tr]=Int_t(pptrk(1));
- }
- } // end if nptracks
- } // end if pdigit
- } //end loop over clusters
- } // hit loop
- } // track loop
- //Int_t nentr2=list->GetEntriesFast();
- //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
- TTree *fAli=gAlice->TreeK();
- TFile *file=NULL;
-
- if (fAli) file =fAli->GetCurrentFile();
- file->cd();
- } // if addBackground
- delete [] xhit;
- delete [] yhit;
-
- Int_t tracks[10];
- Int_t charges[10];
- Int_t nentries=list->GetEntriesFast();
-
- for (Int_t nent=0;nent<nentries;nent++) {
- AliMUONTransientDigit *address=(AliMUONTransientDigit*)list->At(nent);
- if (address==0) continue;
- Int_t ich=address->fChamber;
- Int_t q=address->fSignal;
- iChamber=(AliMUONChamber*) (*fChambers)[ich];
-//
-// Digit Response (noise, threshold, saturation, ...)
-// if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise;
- AliMUONResponse * response=iChamber->ResponseModel();
- q=response->DigitResponse(q);
-
- if (!q) continue;
-
- digits[0]=address->fPadX;
- digits[1]=address->fPadY;
- digits[2]=q;
- digits[3]=address->fPhysics;
- digits[4]=address->fHit;
-
- TObjArray* trlist=(TObjArray*)address->TrackList();
- Int_t nptracks=trlist->GetEntriesFast();
- //printf("nptracks, trlist %d %p\n",nptracks,trlist);
-
- // this was changed to accomodate the real number of tracks
- if (nptracks > 10) {
- cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
- nptracks=10;
- }
- if (nptracks > 2) {
- printf("Attention - nptracks > 2 %d \n",nptracks);
- printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
- }
- for (Int_t tr=0;tr<nptracks;tr++) {
- TVector *ppP=(TVector*)trlist->At(tr);
- if(!ppP ) printf("ppP - %p\n",ppP);
- TVector &pp =*ppP;
- tracks[tr]=Int_t(pp(0));
- charges[tr]=Int_t(pp(1));
- //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
- } //end loop over list of tracks for one pad
- // Sort list of tracks according to charge
- if (nptracks > 1) {
- SortTracks(tracks,charges,nptracks);
- }
- if (nptracks < 10 ) {
- for (Int_t i=nptracks; i<10; i++) {
- tracks[i]=0;
- charges[i]=0;
- }
- }
-
- // fill digits
- pMUON->AddDigits(ich,tracks,charges,digits);
- // delete trlist;
- }
- //cout<<"I'm out of the loops for digitisation"<<endl;
- // gAlice->GetEvent(nev);
- gAlice->TreeD()->Fill();
- pMUON->ResetDigits();
- list->Delete();
-
-
- for(Int_t ii=0;ii<AliMUONConstants::NCh();++ii) {
- if (hitMap[ii]) {
- hm=hitMap[ii];
- delete hm;
- hitMap[ii]=0;
- }
- }
- delete [] nmuon;
- } //end loop over cathodes
- delete [] hitMap;
- char hname[30];
- sprintf(hname,"TreeD%d",nev);
- gAlice->TreeD()->Write(hname);
- // reset tree
- gAlice->TreeD()->Reset();
- delete list;
-
- pAddress->Delete();
- // gObjectTable->Print();