1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 Revision 1.17 2000/06/09 14:58:37 jbarbosa
19 New digitisation per particle type
21 Revision 1.16 2000/04/19 12:55:43 morsch
22 Newly structured and updated version (JB, AM)
27 ////////////////////////////////////////////////
28 // Manager and hits classes for set:RICH //
29 ////////////////////////////////////////////////
37 #include <TObjArray.h>
40 #include <TParticle.h>
44 #include "AliRICHSegmentation.h"
45 #include "AliRICHHit.h"
46 #include "AliRICHCerenkov.h"
47 #include "AliRICHPadHit.h"
48 #include "AliRICHDigit.h"
49 #include "AliRICHTransientDigit.h"
50 #include "AliRICHRawCluster.h"
51 #include "AliRICHRecHit.h"
52 #include "AliRICHHitMapA1.h"
53 #include "AliRICHClusterFinder.h"
56 #include "AliPoints.h"
57 #include "AliCallf77.h"
60 // Static variables for the pad-hit iterator routines
61 static Int_t sMaxIterPad=0;
62 static Int_t sCurIterPad=0;
63 static TClonesArray *fClusters2;
64 static TClonesArray *fHits2;
69 //___________________________________________
72 // Default constructor for RICH manager class
84 //___________________________________________
85 AliRICH::AliRICH(const char *name, const char *title)
86 : AliDetector(name,title)
90 <img src="gif/alirich.gif">
94 fHits = new TClonesArray("AliRICHHit",1000 );
95 gAlice->AddHitList(fHits);
96 fPadHits = new TClonesArray("AliRICHPadHit",100000);
97 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
98 gAlice->AddHitList(fCerenkovs);
99 //gAlice->AddHitList(fHits);
104 fNdch = new Int_t[kNCH];
106 fDchambers = new TObjArray(kNCH);
108 fRecHits = new TObjArray(kNCH);
112 for (i=0; i<kNCH ;i++) {
113 (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
117 fNrawch = new Int_t[kNCH];
119 fRawClusters = new TObjArray(kNCH);
120 //printf("Created fRwClusters with adress:%p",fRawClusters);
122 for (i=0; i<kNCH ;i++) {
123 (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
127 fNrechits = new Int_t[kNCH];
129 for (i=0; i<kNCH ;i++) {
130 (*fRecHits)[i] = new TClonesArray("AliRICHRecHit",1000);
132 //printf("Created fRecHits with adress:%p",fRecHits);
135 SetMarkerColor(kRed);
138 AliRICH::AliRICH(const AliRICH& RICH)
144 //___________________________________________
148 // Destructor of RICH manager class
156 //___________________________________________
157 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
161 // Adds a hit to the Hits list
164 TClonesArray &lhits = *fHits;
165 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
167 //_____________________________________________________________________________
168 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
172 // Adds a RICH cerenkov hit to the Cerenkov Hits list
175 TClonesArray &lcerenkovs = *fCerenkovs;
176 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
177 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
179 //___________________________________________
180 void AliRICH::AddPadHit(Int_t *clhits)
184 // Add a RICH pad hit to the list
187 TClonesArray &lPadHits = *fPadHits;
188 new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
190 //_____________________________________________________________________________
191 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
195 // Add a RICH digit to the list
198 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
199 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
202 //_____________________________________________________________________________
203 void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
206 // Add a RICH digit to the list
209 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
210 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
213 //_____________________________________________________________________________
214 void AliRICH::AddRecHit(Int_t id, Float_t *rechit)
218 // Add a RICH reconstructed hit to the list
221 TClonesArray &lrec = *((TClonesArray*)(*fRecHits)[id]);
222 new(lrec[fNrechits[id]++]) AliRICHRecHit(id,rechit);
225 //___________________________________________
226 void AliRICH::BuildGeometry()
231 // Builds a TNode geometry for event display
235 const int kColorRICH = kGreen;
237 top=gAlice->GetGeometry()->GetNode("alice");
240 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
243 Float_t pos1[3]={0,471.8999,165.2599};
244 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
245 new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
246 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
249 node->SetLineColor(kColorRICH);
253 Float_t pos2[3]={171,470,0};
254 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
255 new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
256 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
259 node->SetLineColor(kColorRICH);
262 Float_t pos3[3]={0,500,0};
263 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
264 new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
265 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
268 node->SetLineColor(kColorRICH);
271 Float_t pos4[3]={-171,470,0};
272 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
273 new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
274 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
277 node->SetLineColor(kColorRICH);
280 Float_t pos5[3]={161.3999,443.3999,-165.3};
281 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
282 new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
283 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
285 node->SetLineColor(kColorRICH);
288 Float_t pos6[3]={0., 471.9, -165.3,};
289 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
290 new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
291 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
294 node->SetLineColor(kColorRICH);
297 Float_t pos7[3]={-161.399,443.3999,-165.3};
298 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
299 new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
300 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
301 node->SetLineColor(kColorRICH);
306 //___________________________________________
307 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
315 //___________________________________________
316 void AliRICH::MakeBranch(Option_t* option)
318 // Create Tree branches for the RICH.
320 const Int_t kBufferSize = 4000;
324 AliDetector::MakeBranch(option);
325 sprintf(branchname,"%sCerenkov",GetName());
326 if (fCerenkovs && gAlice->TreeH()) {
327 gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize);
328 printf("Making Branch %s for Cerenkov Hits\n",branchname);
331 sprintf(branchname,"%sPadHits",GetName());
332 if (fPadHits && gAlice->TreeH()) {
333 gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
334 printf("Making Branch %s for PadHits\n",branchname);
337 // one branch for digits per chamber
340 for (i=0; i<kNCH ;i++) {
341 sprintf(branchname,"%sDigits%d",GetName(),i+1);
343 if (fDchambers && gAlice->TreeD()) {
344 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
345 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
349 // one branch for raw clusters per chamber
350 for (i=0; i<kNCH ;i++) {
351 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
353 if (fRawClusters && gAlice->TreeR()) {
354 gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
355 printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
359 // one branch for rec hits per chamber
360 for (i=0; i<kNCH ;i++) {
361 sprintf(branchname,"%sRecHits%d",GetName(),i+1);
363 if (fRecHits && gAlice->TreeR()) {
364 gAlice->TreeR()->Branch(branchname,&((*fRecHits)[i]), kBufferSize);
365 printf("Making Branch %s for rec. hits in chamber %d\n",branchname,i+1);
370 //___________________________________________
371 void AliRICH::SetTreeAddress()
373 // Set branch address for the Hits and Digits Tree.
377 AliDetector::SetTreeAddress();
380 TTree *treeH = gAlice->TreeH();
381 TTree *treeD = gAlice->TreeD();
382 TTree *treeR = gAlice->TreeR();
386 branch = treeH->GetBranch("RICHPadHits");
387 if (branch) branch->SetAddress(&fPadHits);
390 branch = treeH->GetBranch("RICHCerenkov");
391 if (branch) branch->SetAddress(&fCerenkovs);
396 for (int i=0; i<kNCH; i++) {
397 sprintf(branchname,"%sDigits%d",GetName(),i+1);
399 branch = treeD->GetBranch(branchname);
400 if (branch) branch->SetAddress(&((*fDchambers)[i]));
405 for (i=0; i<kNCH; i++) {
406 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
408 branch = treeR->GetBranch(branchname);
409 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
413 for (i=0; i<kNCH; i++) {
414 sprintf(branchname,"%sRecHits%d",GetName(),i+1);
416 branch = treeR->GetBranch(branchname);
417 if (branch) branch->SetAddress(&((*fRecHits)[i]));
423 //___________________________________________
424 void AliRICH::ResetHits()
426 // Reset number of clusters and the cluster array for this detector
427 AliDetector::ResetHits();
430 if (fPadHits) fPadHits->Clear();
431 if (fCerenkovs) fCerenkovs->Clear();
435 //____________________________________________
436 void AliRICH::ResetDigits()
439 // Reset number of digits and the digits array for this detector
441 for ( int i=0;i<kNCH;i++ ) {
442 if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
443 if (fNdch) fNdch[i]=0;
447 //____________________________________________
448 void AliRICH::ResetRawClusters()
451 // Reset number of raw clusters and the raw clust array for this detector
453 for ( int i=0;i<kNCH;i++ ) {
454 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
455 if (fNrawch) fNrawch[i]=0;
459 //____________________________________________
460 void AliRICH::ResetRecHits()
463 // Reset number of raw clusters and the raw clust array for this detector
466 for ( int i=0;i<kNCH;i++ ) {
467 if ((*fRecHits)[i]) ((TClonesArray*)(*fRecHits)[i])->Clear();
468 if (fNrechits) fNrechits[i]=0;
472 //___________________________________________
473 void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
477 // Setter for the RICH geometry model
481 ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
484 //___________________________________________
485 void AliRICH::SetSegmentationModel(Int_t id, AliRICHSegmentation *segmentation)
489 // Setter for the RICH segmentation model
492 ((AliRICHChamber*) (*fChambers)[id])->SegmentationModel(segmentation);
495 //___________________________________________
496 void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
500 // Setter for the RICH response model
503 ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
506 void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
510 // Setter for the RICH reconstruction model (clusters)
513 ((AliRICHChamber*) (*fChambers)[id])->ReconstructionModel(reconst);
516 void AliRICH::SetNsec(Int_t id, Int_t nsec)
520 // Sets the number of padplanes
523 ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
527 //___________________________________________
529 void AliRICH::StepManager()
532 // Dummy step manager (should never be called)
536 void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
540 // Loop on chambers and on cathode planes
542 for (Int_t icat=1;icat<2;icat++) {
543 gAlice->ResetDigits();
544 gAlice->TreeD()->GetEvent(1); // spurious +1 ...
545 for (Int_t ich=0;ich<kNCH;ich++) {
546 AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
547 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
548 if (pRICHdigits == 0)
551 // Get ready the current chamber stuff
553 AliRICHResponse* response = iChamber->GetResponseModel();
554 AliRICHSegmentation* seg = iChamber->GetSegmentationModel();
555 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
557 rec->SetSegmentation(seg);
558 rec->SetResponse(response);
559 rec->SetDigits(pRICHdigits);
560 rec->SetChamber(ich);
561 if (nev==0) rec->CalibrateCOG();
562 rec->FindRawClusters();
565 fRch=RawClustAddress(ich);
569 gAlice->TreeR()->Fill();
571 for (int i=0;i<kNCH;i++) {
572 fRch=RawClustAddress(i);
573 int nraw=fRch->GetEntriesFast();
574 printf ("Chamber %d, raw clusters %d\n",i,nraw);
582 sprintf(hname,"TreeR%d",nev);
583 gAlice->TreeR()->Write(hname);
584 gAlice->TreeR()->Reset();
586 //gObjectTable->Print();
590 //______________________________________________________________________________
591 void AliRICH::Streamer(TBuffer &R__b)
593 // Stream an object of class AliRICH.
594 AliRICHChamber *iChamber;
595 AliRICHSegmentation *segmentation;
596 AliRICHResponse *response;
597 TClonesArray *digitsaddress;
598 TClonesArray *rawcladdress;
599 TClonesArray *rechitaddress;
601 if (R__b.IsReading()) {
602 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
603 AliDetector::Streamer(R__b);
605 R__b >> fPadHits; // diff
607 R__b >> fCerenkovs; // diff
609 R__b >> fRawClusters;
610 R__b >> fRecHits; //diff
611 R__b.ReadArray(fNdch);
612 R__b.ReadArray(fNrawch);
613 R__b.ReadArray(fNrechits);
616 // Stream chamber related information
617 for (Int_t i =0; i<kNCH; i++) {
618 iChamber=(AliRICHChamber*) (*fChambers)[i];
619 iChamber->Streamer(R__b);
620 segmentation=iChamber->GetSegmentationModel();
621 segmentation->Streamer(R__b);
622 response=iChamber->GetResponseModel();
623 response->Streamer(R__b);
624 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
625 rawcladdress->Streamer(R__b);
626 rechitaddress=(TClonesArray*) (*fRecHits)[i];
627 rechitaddress->Streamer(R__b);
628 digitsaddress=(TClonesArray*) (*fDchambers)[i];
629 digitsaddress->Streamer(R__b);
633 R__b.WriteVersion(AliRICH::IsA());
634 AliDetector::Streamer(R__b);
636 R__b << fPadHits; // diff
638 R__b << fCerenkovs; // diff
640 R__b << fRawClusters;
641 R__b << fRecHits; //diff
642 R__b.WriteArray(fNdch, kNCH);
643 R__b.WriteArray(fNrawch, kNCH);
644 R__b.WriteArray(fNrechits, kNCH);
647 // Stream chamber related information
648 for (Int_t i =0; i<kNCH; i++) {
649 iChamber=(AliRICHChamber*) (*fChambers)[i];
650 iChamber->Streamer(R__b);
651 segmentation=iChamber->GetSegmentationModel();
652 segmentation->Streamer(R__b);
653 response=iChamber->GetResponseModel();
654 response->Streamer(R__b);
655 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
656 rawcladdress->Streamer(R__b);
657 rechitaddress=(TClonesArray*) (*fRecHits)[i];
658 rechitaddress->Streamer(R__b);
659 digitsaddress=(TClonesArray*) (*fDchambers)[i];
660 digitsaddress->Streamer(R__b);
664 AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
667 // Initialise the pad iterator
668 // Return the address of the first padhit for hit
669 TClonesArray *theClusters = clusters;
670 Int_t nclust = theClusters->GetEntriesFast();
671 if (nclust && hit->fPHlast > 0) {
672 sMaxIterPad=Int_t(hit->fPHlast);
673 sCurIterPad=Int_t(hit->fPHfirst);
674 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
681 AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters)
684 // Iterates over pads
687 if (sCurIterPad <= sMaxIterPad) {
688 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
695 void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
697 // keep galice.root for signal and name differently the file for
698 // background when add! otherwise the track info for signal will be lost !
700 static Bool_t first=kTRUE;
702 char *addBackground = strstr(option,"Add");
704 FILE* points; //these will be the digits...
706 points=fopen("points.dat","w");
708 AliRICHChamber* iChamber;
709 AliRICHSegmentation* segmentation;
714 TObjArray *list=new TObjArray;
715 static TClonesArray *pAddress=0;
716 if(!pAddress) pAddress=new TClonesArray("TVector",1000);
719 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
720 AliRICHHitMap* pHitMap[10];
722 for (i=0; i<10; i++) {pHitMap[i]=0;}
723 if (addBackground ) {
726 cout<<"filename"<<fFileName<<endl;
727 pFile=new TFile(fFileName);
728 cout<<"I have opened "<<fFileName<<" file "<<endl;
729 fHits2 = new TClonesArray("AliRICHHit",1000 );
730 fClusters2 = new TClonesArray("AliRICHPadHit",10000);
734 // Get Hits Tree header from file
735 if(fHits2) fHits2->Clear();
736 if(fClusters2) fClusters2->Clear();
737 if(TrH1) delete TrH1;
741 sprintf(treeName,"TreeH%d",nev);
742 TrH1 = (TTree*)gDirectory->Get(treeName);
744 printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
746 // Set branch addresses
749 sprintf(branchname,"%s",GetName());
750 if (TrH1 && fHits2) {
751 branch = TrH1->GetBranch(branchname);
752 if (branch) branch->SetAddress(&fHits2);
754 if (TrH1 && fClusters2) {
755 branch = TrH1->GetBranch("RICHCluster");
756 if (branch) branch->SetAddress(&fClusters2);
760 // loop over cathodes
764 for (int icat=0; icat<1; icat++) {
766 for (i =0; i<kNCH; i++) {
767 iChamber=(AliRICHChamber*) (*fChambers)[i];
768 if (iChamber->Nsec()==1 && icat==1) {
771 segmentation=iChamber->GetSegmentationModel(icat+1);
773 pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
779 TTree *treeH = gAlice->TreeH();
780 Int_t ntracks =(Int_t) treeH->GetEntries();
781 for (Int_t track=0; track<ntracks; track++) {
783 treeH->GetEvent(track);
786 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
788 mHit=(AliRICHHit*)pRICH->NextHit())
793 Int_t nch = mHit->fChamber-1; // chamber number
794 if (nch >kNCH) continue;
795 iChamber = &(pRICH->Chamber(nch));
797 TParticle *current = (TParticle*)(*gAlice->Particles())[track];
799 Int_t particle = current->GetPdgCode();
801 //printf("Flag:%d\n",flag);
802 //printf("Track:%d\n",track);
803 //printf("Particle:%d\n",particle);
809 if(TMath::Abs(particle) == 211 || TMath::Abs(particle) == 111)
813 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
814 || TMath::Abs(particle)==311)
817 if (flag == 3 && TMath::Abs(particle)==2212)
820 if (flag == 4 && TMath::Abs(particle)==13)
823 if (flag == 5 && TMath::Abs(particle)==11)
826 if (flag == 6 && TMath::Abs(particle)==2112)
830 //printf ("Particle: %d, Flag: %d, Digitse: %d\n",particle,flag,digitse);
837 // Loop over pad hits
838 for (AliRICHPadHit* mPad=
839 (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
841 mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
843 Int_t cathode = mPad->fCathode; // cathode number
844 Int_t ipx = mPad->fPadX; // pad number on X
845 Int_t ipy = mPad->fPadY; // pad number on Y
846 Int_t iqpad = mPad->fQpad; // charge per pad
849 //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
852 segmentation=iChamber->GetSegmentationModel(cathode);
853 segmentation->GetPadCxy(ipx,ipy,thex,they);
854 new((*pAddress)[countadr++]) TVector(2);
855 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
856 trinfo(0)=(Float_t)track;
857 trinfo(1)=(Float_t)iqpad;
863 AliRICHTransientDigit* pdigit;
864 // build the list of fired pads and update the info
865 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
866 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
867 pHitMap[nch]->SetHit(ipx, ipy, counter);
869 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
871 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
872 trlist->Add(&trinfo);
874 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
876 (*pdigit).fSignal+=iqpad;
877 // update list of tracks
878 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
879 Int_t lastEntry=trlist->GetLast();
880 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
881 TVector &ptrk=*ptrkP;
882 Int_t lastTrack=Int_t(ptrk(0));
883 Int_t lastCharge=Int_t(ptrk(1));
884 if (lastTrack==track) {
886 trlist->RemoveAt(lastEntry);
888 trinfo(1)=lastCharge;
889 trlist->AddAt(&trinfo,lastEntry);
891 trlist->Add(&trinfo);
893 // check the track list
894 Int_t nptracks=trlist->GetEntriesFast();
896 printf("Attention - tracks: %d (>2)\n",nptracks);
897 //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
898 for (Int_t tr=0;tr<nptracks;tr++) {
899 TVector *pptrkP=(TVector*)trlist->At(tr);
900 TVector &pptrk=*pptrkP;
901 trk[tr]=Int_t(pptrk(0));
902 chtrk[tr]=Int_t(pptrk(1));
906 } //end loop over clusters
907 }// track type condition
911 // open the file with background
913 if (addBackground ) {
914 ntracks =(Int_t)TrH1->GetEntries();
915 //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
916 //printf("background - Start loop over tracks \n");
920 for (Int_t trak=0; trak<ntracks; trak++) {
921 if (fHits2) fHits2->Clear();
922 if (fClusters2) fClusters2->Clear();
923 TrH1->GetEvent(trak);
927 for(int j=0;j<fHits2->GetEntriesFast();++j)
929 mHit=(AliRICHHit*) (*fHits2)[j];
930 Int_t nch = mHit->fChamber-1; // chamber number
931 if (nch >6) continue;
932 iChamber = &(pRICH->Chamber(nch));
933 Int_t rmin = (Int_t)iChamber->RInner();
934 Int_t rmax = (Int_t)iChamber->ROuter();
936 // Loop over pad hits
937 for (AliRICHPadHit* mPad=
938 (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
940 mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2))
942 Int_t cathode = mPad->fCathode; // cathode number
943 Int_t ipx = mPad->fPadX; // pad number on X
944 Int_t ipy = mPad->fPadY; // pad number on Y
945 Int_t iqpad = mPad->fQpad; // charge per pad
946 if (trak==3 && nch==0 && icat==0) printf("bgr - trak,iqpad,ipx,ipy %d %d %d %d\n",trak,iqpad,ipx,ipy);
950 segmentation=iChamber->GetSegmentationModel(cathode);
951 segmentation->GetPadCxy(ipx,ipy,thex,they);
952 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
953 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
954 new((*pAddress)[countadr++]) TVector(2);
955 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
956 trinfo(0)=-1; // tag background
961 if (trak <4 && icat==0 && nch==0)
962 printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
963 pHitMap[nch]->TestHit(ipx, ipy),trak);
964 AliRICHTransientDigit* pdigit;
965 // build the list of fired pads and update the info
966 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
967 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
969 pHitMap[nch]->SetHit(ipx, ipy, counter);
971 printf("bgr new elem in list - counter %d\n",counter);
973 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
975 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
976 trlist->Add(&trinfo);
978 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
980 (*pdigit).fSignal+=iqpad;
981 // update list of tracks
982 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
983 Int_t lastEntry=trlist->GetLast();
984 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
985 TVector &ptrk=*ptrkP;
986 Int_t lastTrack=Int_t(ptrk(0));
990 trlist->Add(&trinfo);
992 // check the track list
993 Int_t nptracks=trlist->GetEntriesFast();
995 for (Int_t tr=0;tr<nptracks;tr++) {
996 TVector *pptrkP=(TVector*)trlist->At(tr);
997 TVector &pptrk=*pptrkP;
998 trk[tr]=Int_t(pptrk(0));
999 chtrk[tr]=Int_t(pptrk(1));
1001 } // end if nptracks
1003 } //end loop over clusters
1006 TTree *fAli=gAlice->TreeK();
1007 if (fAli) pFile =fAli->GetCurrentFile();
1013 //cout<<"Start filling digits \n "<<endl;
1014 Int_t nentries=list->GetEntriesFast();
1015 //printf(" \n \n nentries %d \n",nentries);
1017 // start filling the digits
1019 for (Int_t nent=0;nent<nentries;nent++) {
1020 AliRICHTransientDigit *address=(AliRICHTransientDigit*)list->At(nent);
1021 if (address==0) continue;
1022 Int_t ich=address->fChamber;
1023 Int_t q=address->fSignal;
1024 iChamber=(AliRICHChamber*) (*fChambers)[ich];
1025 AliRICHResponse * response=iChamber->GetResponseModel();
1026 Int_t adcmax= (Int_t) response->MaxAdc();
1027 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
1028 Float_t meanNoise = gRandom->Gaus(1.7, 0.25);
1029 Float_t noise = gRandom->Gaus(0, meanNoise);
1031 // magic number to be parametrised !!!
1032 if ( q <= 6.8) continue;
1033 if ( q >= adcmax) q=adcmax;
1034 digits[0]=address->fPadX;
1035 digits[1]=address->fPadY;
1038 TObjArray* trlist=(TObjArray*)address->TrackList();
1039 Int_t nptracks=trlist->GetEntriesFast();
1041 // this was changed to accomodate the real number of tracks
1042 if (nptracks > 10) {
1043 cout<<"Attention - tracks > 10 "<<nptracks<<endl;
1047 printf("Attention - tracks > 2 %d \n",nptracks);
1048 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
1049 //icat,ich,digits[0],digits[1],q);
1051 for (Int_t tr=0;tr<nptracks;tr++) {
1052 TVector *ppP=(TVector*)trlist->At(tr);
1054 tracks[tr]=Int_t(pp(0));
1055 charges[tr]=Int_t(pp(1));
1056 } //end loop over list of tracks for one pad
1057 if (nptracks < 10 ) {
1058 for (Int_t t=nptracks; t<10; t++) {
1065 fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
1068 pRICH->AddDigits(ich,tracks,charges,digits);
1070 gAlice->TreeD()->Fill();
1073 for(Int_t ii=0;ii<kNCH;++ii) {
1081 //TTree *TD=gAlice->TreeD();
1082 //Stat_t ndig=TD->GetEntries();
1083 //cout<<"number of digits "<<ndig<<endl;
1085 for (int k=0;k<kNCH;k++) {
1086 fDch= pRICH->DigitsAddress(k);
1087 int ndigit=fDch->GetEntriesFast();
1088 printf ("Chamber %d digits %d \n",k,ndigit);
1090 pRICH->ResetDigits();
1091 } //end loop over cathodes
1093 sprintf(hname,"TreeD%d",nev);
1094 gAlice->TreeD()->Write(hname);
1097 // gAlice->TreeD()->Reset();
1100 // gObjectTable->Print();
1103 AliRICH& AliRICH::operator=(const AliRICH& rhs)
1105 // Assignment operator
1110 Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
1113 // Calls the charge disintegration method of the current chamber and adds
1114 // the simulated cluster to the root treee
1117 Float_t newclust[6][500];
1121 // Integrated pulse height on chamber
1125 ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
1130 for (Int_t i=0; i<nnew; i++) {
1131 if (Int_t(newclust[3][i]) > 0) {
1134 clhits[1] = Int_t(newclust[5][i]);
1136 clhits[2] = Int_t(newclust[0][i]);
1138 clhits[3] = Int_t(newclust[1][i]);
1140 clhits[4] = Int_t(newclust[2][i]);
1142 clhits[5] = Int_t(newclust[3][i]);
1143 // Pad: chamber sector
1144 clhits[6] = Int_t(newclust[4][i]);