Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/*
17$Log$
18*/
19
ddae0931 20////////////////////////////////////////////////
21// Manager and hits classes for set:RICH //
22////////////////////////////////////////////////
fe4da5cc 23
ddae0931 24#include <TBRIK.h>
25#include <TTUBE.h>
fe4da5cc 26#include <TNode.h>
27#include <TRandom.h>
ddae0931 28#include <TObject.h>
29#include <TVector.h>
30#include <TObjArray.h>
fe4da5cc 31
fe4da5cc 32#include "AliRICH.h"
ddae0931 33#include "AliRICHHitMap.h"
fe4da5cc 34#include "AliRun.h"
ddae0931 35#include "AliMC.h"
36#include "iostream.h"
37#include "AliCallf77.h"
38
39// Static variables for the pad-hit iterator routines
40static Int_t sMaxIterPad=0;
41static Int_t sCurIterPad=0;
42static TClonesArray *fClusters2;
43static TClonesArray *fHits2;
44
fe4da5cc 45ClassImp(AliRICH)
ddae0931 46
47//___________________________________________
fe4da5cc 48AliRICH::AliRICH()
49{
ddae0931 50 fIshunt = 0;
51 fHits = 0;
52 fClusters = 0;
53 fNclusters = 0;
54 fDchambers = 0;
55 fRecClusters= 0;
56 fCerenkovs = 0;
57 fNdch = 0;
fe4da5cc 58}
ddae0931 59
60//___________________________________________
fe4da5cc 61AliRICH::AliRICH(const char *name, const char *title)
ddae0931 62 : AliDetector(name,title)
fe4da5cc 63{
ddae0931 64//Begin_Html
65/*
66 <img src="gif/alirich.gif">
67*/
68//End_Html
69
70 fHits = new TClonesArray("AliRICHhit",1000 );
71 fClusters = new TClonesArray("AliRICHcluster",10000);
72 fCerenkovs = new TClonesArray("AliRICHCerenkov",20000);
73 fNclusters = 0;
74 fIshunt = 0;
75
76 fNdch = new Int_t[7];
77
78 fDchambers = new TObjArray(7);
79
80 Int_t i;
81
82 for (i=0; i<7 ;i++) {
83 (*fDchambers)[i] = new TClonesArray("AliRICHdigit",10000);
84 fNdch[i]=0;
85 }
86
87 fRecClusters=new TObjArray(7);
88 for (i=0; i<7;i++)
89 (*fRecClusters)[i] = new TObjArray(1000);
90
91//
92// Transport angular cut
93 fAccCut=0;
94 fAccMin=2;
95 fAccMax=9;
96
97 SetMarkerColor(kRed);
fe4da5cc 98}
99
ddae0931 100//___________________________________________
fe4da5cc 101AliRICH::~AliRICH()
102{
ddae0931 103 fIshunt = 0;
104 delete fHits;
105 delete fClusters;
106 delete fCerenkovs;
107 for (Int_t i=0;i<7;i++)
108 delete (*fRecClusters)[i];
109 delete fRecClusters;
110
fe4da5cc 111}
112
ddae0931 113//___________________________________________
fe4da5cc 114void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
115{
ddae0931 116 TClonesArray &lhits = *fHits;
117 new(lhits[fNhits++]) AliRICHhit(fIshunt,track,vol,hits);
fe4da5cc 118}
fe4da5cc 119//_____________________________________________________________________________
ddae0931 120void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
121{
122 TClonesArray &lcerenkovs = *fCerenkovs;
123 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
fe4da5cc 124}
ddae0931 125//___________________________________________
126void AliRICH::AddCluster(Int_t *clhits)
127{
128 TClonesArray &lclusters = *fClusters;
129 new(lclusters[fNclusters++]) AliRICHcluster(clhits);
fe4da5cc 130}
fe4da5cc 131//_____________________________________________________________________________
ddae0931 132void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
133{
134 //
135 // Add a RICH digit to the list
136 //
137
138 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
139 new(ldigits[fNdch[id]++]) AliRICHdigit(tracks,charges,digits);
fe4da5cc 140}
141
ddae0931 142
fe4da5cc 143//_____________________________________________________________________________
ddae0931 144void AliRICH::AddRecCluster(Int_t iCh, Int_t iCat, AliRICHRecCluster* Cluster)
fe4da5cc 145{
ddae0931 146 //
147 // Add a RICH reconstructed cluster to the list
148 //
149 TObjArray* ClusterList = RecClusters(iCh,iCat);
150 ClusterList->Add(Cluster);
fe4da5cc 151}
152
fe4da5cc 153
154
ddae0931 155//___________________________________________
156void AliRICH::BuildGeometry()
157
fe4da5cc 158{
ddae0931 159 //
160 // Builds a TNode geometry for event display
161 //
162 TNode *Node, *Top;
163
164 const int kColorRICH = kGreen;
165 //
166 Top=gAlice->GetGeometry()->GetNode("alice");
167
168
169 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
170
171 Top->cd();
172 Float_t pos1[3]={0,471.8999,165.2599};
173 Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90));
174 Node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
175
176
177 Node->SetLineColor(kColorRICH);
178 fNodes->Add(Node);
179 Top->cd();
180
181 Float_t pos2[3]={171,470,0};
182 Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90,-20,90,70,0,0));
183 Node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
184
185
186 Node->SetLineColor(kColorRICH);
187 fNodes->Add(Node);
188 Top->cd();
189 Float_t pos3[3]={0,500,0};
190 Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90,0,90,90,0,0));
191 Node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
192
193
194 Node->SetLineColor(kColorRICH);
195 fNodes->Add(Node);
196 Top->cd();
197 Float_t pos4[3]={-171,470,0};
198 Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2], new TRotMatrix("rot996","rot996",90,20,90,110,0,0));
199 Node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
200
201
202 Node->SetLineColor(kColorRICH);
203 fNodes->Add(Node);
204 Top->cd();
205 Float_t pos5[3]={161.3999,443.3999,-165.3};
206 Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70));
207 Node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
208
209 Node->SetLineColor(kColorRICH);
210 fNodes->Add(Node);
211 Top->cd();
212 Float_t pos6[3]={0., 471.9, -165.3,};
213 Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90));
214 Node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
215
216
217 Node->SetLineColor(kColorRICH);
218 fNodes->Add(Node);
219 Top->cd();
220 Float_t pos7[3]={-161.399,443.3999,-165.3};
221 Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110));
222 Node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
223 Node->SetLineColor(kColorRICH);
224 fNodes->Add(Node);
225
fe4da5cc 226}
227
ddae0931 228//___________________________________________
229Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
230{
231 return 9999;
232}
fe4da5cc 233
ddae0931 234//___________________________________________
235void AliRICH::MakeBranch(Option_t* option)
fe4da5cc 236{
ddae0931 237 // Create Tree branches for the RICH.
fe4da5cc 238
ddae0931 239 const Int_t buffersize = 4000;
240 char branchname[20];
fe4da5cc 241
ddae0931 242
243 AliDetector::MakeBranch(option);
244 sprintf(branchname,"%sCerenkov",GetName());
245 if (fCerenkovs && gAlice->TreeH()) {
246 gAlice->TreeH()->Branch(branchname,&fCerenkovs, buffersize);
247 printf("Making Branch %s for Cerenkov Hits\n",branchname);
fe4da5cc 248 }
ddae0931 249
250 sprintf(branchname,"%sCluster",GetName());
251 if (fClusters && gAlice->TreeH()) {
252 gAlice->TreeH()->Branch(branchname,&fClusters, buffersize);
253 printf("Making Branch %s for clusters\n",branchname);
fe4da5cc 254 }
255
ddae0931 256// one branch for digits per chamber
257 Int_t i;
258
259 for (i=0; i<7 ;i++) {
260 sprintf(branchname,"%sDigits%d",GetName(),i+1);
261
262 if (fDchambers && gAlice->TreeD()) {
263 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), buffersize);
264 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
265 }
fe4da5cc 266 }
ddae0931 267// one branch for rec clusters
268 for (i=0; i<7; i++) {
269 sprintf(branchname,"%sRecClus%d",GetName(),i+1);
270 if (fRecClusters && gAlice->TreeD()) {
271 gAlice->TreeR()
272 ->Branch(branchname,"TObjArray",
273 &((*fRecClusters)[i]), buffersize,0);
274 printf("Making Branch %s for clusters in chamber %d\n",
275 branchname,i+1);
fe4da5cc 276 }
ddae0931 277 }
fe4da5cc 278}
279
ddae0931 280//___________________________________________
281void AliRICH::SetTreeAddress()
fe4da5cc 282{
ddae0931 283 // Set branch address for the Hits and Digits Tree.
284 char branchname[20];
285 AliDetector::SetTreeAddress();
286
287 TBranch *branch;
288 TTree *treeH = gAlice->TreeH();
289 TTree *treeD = gAlice->TreeD();
290
291 if (treeH) {
292 if (fClusters) {
293 branch = treeH->GetBranch("RICHCluster");
294 if (branch) branch->SetAddress(&fClusters);
fe4da5cc 295 }
ddae0931 296 if (fCerenkovs) {
297 branch = treeH->GetBranch("RICHCerenkov");
298 if (branch) branch->SetAddress(&fCerenkovs);
fe4da5cc 299 }
ddae0931 300
301 }
302
303 if (treeD) {
304 for (int i=0; i<7; i++) {
305 sprintf(branchname,"%sDigits%d",GetName(),i+1);
306 if (fDchambers) {
307 branch = treeD->GetBranch(branchname);
308 if (branch) branch->SetAddress(&((*fDchambers)[i]));
309 }
fe4da5cc 310 }
fe4da5cc 311 }
ddae0931 312}
313//___________________________________________
314void AliRICH::ResetHits()
315{
316 // Reset number of clusters and the cluster array for this detector
317 AliDetector::ResetHits();
318 fNclusters = 0;
319 if (fClusters) fClusters->Clear();
320 if (fCerenkovs) fCerenkovs->Clear();
fe4da5cc 321}
322
ddae0931 323//____________________________________________
324void AliRICH::ResetDigits()
325{
326 //
327 // Reset number of digits and the digits array for this detector
328 //
329 for ( int i=0;i<7;i++ ) {
330 if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
331 if (fNdch) fNdch[i]=0;
332 }
333}
334//____________________________________________
335void AliRICH::ResetRecClusters()
fe4da5cc 336{
ddae0931 337 //
338 // Reset the rec clusters
339 //
340 for ( int i=0;i<7;i++ ) {
341 if ((*fRecClusters)[i]) (*fRecClusters)[i]->Clear();
fe4da5cc 342 }
ddae0931 343}
344//___________________________________________
fe4da5cc 345
ddae0931 346void AliRICH::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2)
fe4da5cc 347{
ddae0931 348 Int_t i=2*(id-1);
349 ((AliRICHchamber*) (*fChambers)[i]) ->SetPADSIZ(isec,p1,p2);
350 ((AliRICHchamber*) (*fChambers)[i+1])->SetPADSIZ(isec,p1,p2);
351}
fe4da5cc 352
ddae0931 353//___________________________________________
354void AliRICH::SetMUCHSP(Int_t id, Float_t p1)
fe4da5cc 355{
ddae0931 356 Int_t i=2*(id-1);
357 ((AliRICHchamber*) (*fChambers)[i])->SetMUCHSP(p1);
358 ((AliRICHchamber*) (*fChambers)[i+1])->SetMUCHSP(p1);
359}
fe4da5cc 360
ddae0931 361//___________________________________________
362void AliRICH::SetMUSIGM(Int_t id, Float_t p1, Float_t p2)
363{
364 Int_t i=2*(id-1);
365 ((AliRICHchamber*) (*fChambers)[i])->SetMUSIGM(p1,p2);
366 ((AliRICHchamber*) (*fChambers)[i+1])->SetMUSIGM(p1,p2);
367}
fe4da5cc 368
ddae0931 369//___________________________________________
370void AliRICH::SetRSIGM(Int_t id, Float_t p1)
371{
372 Int_t i=2*(id-1);
373 ((AliRICHchamber*) (*fChambers)[i])->SetRSIGM(p1);
374 ((AliRICHchamber*) (*fChambers)[i+1])->SetRSIGM(p1);
375}
fe4da5cc 376
ddae0931 377//___________________________________________
378void AliRICH::SetMAXADC(Int_t id, Float_t p1)
fe4da5cc 379{
ddae0931 380 Int_t i=2*(id-1);
381 ((AliRICHchamber*) (*fChambers)[i])->SetMAXADC(p1);
382 ((AliRICHchamber*) (*fChambers)[i+1])->SetMAXADC(p1);
fe4da5cc 383}
384
ddae0931 385//___________________________________________
386void AliRICH::SetSMAXAR(Float_t p1)
fe4da5cc 387{
ddae0931 388 fMaxStepGas=p1;
fe4da5cc 389}
390
ddae0931 391//___________________________________________
392void AliRICH::SetSMAXAL(Float_t p1)
393{
394 fMaxStepAlu=p1;
395}
fe4da5cc 396
ddae0931 397//___________________________________________
398void AliRICH::SetDMAXAR(Float_t p1)
fe4da5cc 399{
ddae0931 400 fMaxDestepGas=p1;
fe4da5cc 401}
402
ddae0931 403//___________________________________________
404void AliRICH::SetDMAXAL(Float_t p1)
405{
406 fMaxDestepAlu=p1;
407}
408//___________________________________________
409void AliRICH::SetRICHACC(Bool_t acc, Float_t angmin, Float_t angmax)
410{
411 fAccCut=acc;
412 fAccMin=angmin;
413 fAccMax=angmax;
414}
415//___________________________________________
416void AliRICH::SetSegmentationModel(Int_t id, Int_t isec, AliRICHsegmentation *segmentation)
417{
418 ((AliRICHchamber*) (*fChambers)[id])->SegmentationModel(isec, segmentation);
fe4da5cc 419
ddae0931 420}
421//___________________________________________
422void AliRICH::SetResponseModel(Int_t id, Response_t res, AliRICHresponse *response)
fe4da5cc 423{
ddae0931 424 ((AliRICHchamber*) (*fChambers)[id])->ResponseModel(res, response);
fe4da5cc 425}
426
ddae0931 427void AliRICH::SetNsec(Int_t id, Int_t nsec)
428{
429 ((AliRICHchamber*) (*fChambers)[id])->SetNsec(nsec);
430}
431
432
433//___________________________________________
434
435void AliRICH::StepManager()
fe4da5cc 436{
ddae0931 437 printf("Dummy version of RICH step -- it should never happen!!\n");
438 const Float_t kRaddeg = 180/TMath::Pi();
439 AliMC* pMC = AliMC::GetMC();
440 Int_t nsec, ipart;
441 Float_t x[4], p[4];
442 Float_t pt, th0, th1;
443 char proc[5];
444 if(fAccCut) {
445 if((nsec=pMC->NSecondaries())>0) {
446 pMC->ProdProcess(proc);
447 if((pMC->TrackPid()==113 || pMC->TrackPid()==114) && !strcmp(proc,"DCAY")) {
448
449 // Check angular acceptance
450 //* --- and have muons from resonance decays in the wanted window ---
451 if(nsec != 2) {
452 printf(" AliRICH::StepManager: Strange resonance Decay into %d particles\n",nsec);
453 pMC->StopEvent();
454 } else {
455 pMC->GetSecondary(0,ipart,x,p);
456 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
457 th0 = TMath::ATan2(pt,p[2])*kRaddeg;
458 pMC->GetSecondary(1,ipart,x,p);
459 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
460 th1 = TMath::ATan2(pt,p[2])*kRaddeg;
461 if(!(fAccMin < th0 && th0 < fAccMax) ||
462 !(fAccMin < th1 && th1 < fAccMax))
463 pMC->StopEvent();
464 }
465 }
466 }
fe4da5cc 467 }
ddae0931 468}
469void AliRICH::ReconstructClusters()
470{
471//
472// Initialize the necessary correspondance table
473//
f91473f6 474 static const Int_t kMaxNPadx = 600;
475 static const Int_t kMaxNPady = 600;
476 Int_t elem[kMaxNPadx*2][kMaxNPady*2];
ddae0931 477//
478// Loop on chambers and on cathode planes
479//
480 for (Int_t ich=0;ich<7;ich++)
481 for (Int_t icat=0;icat<1;icat++) {
482 //
483 // Get ready the current chamber stuff
484 //
485
486 printf ("Olarilole");
487 AliRICHchamber* iChamber= &(this->Chamber(ich));
488 AliRICHsegmentation* segmentation =
489 iChamber->GetSegmentationModel(icat+1);
490 if (!segmentation)
491 continue;
492 TClonesArray *RICHdigits = this->DigitsAddress(ich);
493 if (RICHdigits == 0)
494 continue;
495 cout << "Npx " << segmentation->Npx()
496 << " Npy " << segmentation->Npy() << endl;
497 //
498 // Ready the digits
499 //
500 gAlice->ResetDigits();
501 gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
502 Int_t ndigits = RICHdigits->GetEntriesFast();
503 if (ndigits == 0)
504 continue;
505 printf("Found %d digits for cathode %d in chamber %d \n",
506 ndigits,icat,ich+1);
507 AliRICHdigit *mdig;
508 AliRICHRecCluster *Cluster;
509 //
510 // Build the correspondance table
511 //
f91473f6 512 memset(elem,0,sizeof(Int_t)*kMaxNPadx*kMaxNPady*4);
ddae0931 513 Int_t digit;
514 for (digit=0; digit<ndigits; digit++)
515 {
516 mdig = (AliRICHdigit*)RICHdigits->UncheckedAt(digit);
f91473f6 517 elem[kMaxNPadx+mdig->fPadX][kMaxNPady+mdig->fPadY] = digit+1;
ddae0931 518 // because default is 0
519 }
520 //
521 // Declare some useful variables
522 //
523 Int_t Xlist[10];
524 Int_t Ylist[10];
525 Int_t Nlist;
526 Int_t nclust=0;
527 //
528 // loop over digits
529 //
530 for (digit=0;digit<ndigits;digit++) {
531 mdig = (AliRICHdigit*)RICHdigits->UncheckedAt(digit);
532 //
533 // if digit still available, start clustering
534 //
f91473f6 535 if (elem[kMaxNPadx+mdig->fPadX][kMaxNPady+mdig->fPadY]) {
ddae0931 536 Cluster = new AliRICHRecCluster(digit, ich, icat);
f91473f6 537 elem[kMaxNPadx+mdig->fPadX][kMaxNPady+mdig->fPadY]=0;
ddae0931 538 //
539 // loop over the current list of digits
540 // and look for neighbours
541 //
542 for(Int_t clusDigit=Cluster->FirstDigitIndex();
543 clusDigit!=Cluster->InvalidDigitIndex();
544 clusDigit=Cluster->NextDigitIndex()) {
545 AliRICHdigit* pDigit=(AliRICHdigit*)RICHdigits
546 ->UncheckedAt(clusDigit);
547 segmentation->Neighbours(pDigit->fPadX,pDigit->fPadY,
548 &Nlist, Xlist, Ylist);
549 for (Int_t Ilist=0;Ilist<Nlist;Ilist++) {
f91473f6 550 if (elem[kMaxNPadx+Xlist[Ilist]][kMaxNPady
ddae0931 551 +Ylist[Ilist]]) {
552 //
553 // Add the digit at the end and reset the table
554 //
f91473f6 555 Cluster->AddDigit(elem[kMaxNPadx+Xlist[Ilist]][kMaxNPady+Ylist[Ilist]]-1);
556 elem[kMaxNPadx+Xlist[Ilist]][kMaxNPady
ddae0931 557 +Ylist[Ilist]] =0;
558 } // if elem
559 } // for Ilist
560 } // for pDigit
561 //
562 // Store the cluster (good time to do Cluster polishing)
563 //
564 segmentation->FitXY(Cluster,RICHdigits);
565 nclust ++;
566 AddRecCluster(ich,icat,Cluster);
567 }
568 }
569 printf("===> %d Clusters\n",nclust);
570 } // for icat
571}
572
573
574//______________________________________________________________________________
575void AliRICH::Streamer(TBuffer &R__b)
576{
577 // Stream an object of class AliRICH.
578 AliRICHchamber *iChamber;
579 AliRICHsegmentation *segmentation;
580 AliRICHresponse *response;
581 TClonesArray *digitsaddress;
582
583 if (R__b.IsReading()) {
584 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
585 AliDetector::Streamer(R__b);
586 R__b >> fNclusters;
587 R__b >> fClusters; // diff
588 R__b >> fDchambers;
589 R__b.ReadArray(fNdch);
590 //
591 R__b >> fAccCut;
592 R__b >> fAccMin;
593 R__b >> fAccMax;
594 //
595 R__b >> fChambers;
596// Stream chamber related information
597 for (Int_t i =0; i<7; i++) {
598 iChamber=(AliRICHchamber*) (*fChambers)[i];
599 iChamber->Streamer(R__b);
600 if (iChamber->Nsec()==1) {
601 segmentation=iChamber->GetSegmentationModel(1);
602 segmentation->Streamer(R__b);
603 } else {
604 segmentation=iChamber->GetSegmentationModel(1);
605 segmentation->Streamer(R__b);
606 segmentation=iChamber->GetSegmentationModel(2);
607 segmentation->Streamer(R__b);
608 }
609 response=iChamber->GetResponseModel(mip);
610 response->Streamer(R__b);
611 response=iChamber->GetResponseModel(cerenkov);
612 response->Streamer(R__b);
613
614 digitsaddress=(TClonesArray*) (*fDchambers)[i];
615 digitsaddress->Streamer(R__b);
616 }
617
618 } else {
619 R__b.WriteVersion(AliRICH::IsA());
620 AliDetector::Streamer(R__b);
621 R__b << fNclusters;
622 R__b << fClusters; // diff
623 R__b << fDchambers;
624 R__b.WriteArray(fNdch, 7);
625 //
626 R__b << fAccCut;
627 R__b << fAccMin;
628 R__b << fAccMax;
629 //
630 R__b << fChambers;
631// Stream chamber related information
632 for (Int_t i =0; i<7; i++) {
633 iChamber=(AliRICHchamber*) (*fChambers)[i];
634 iChamber->Streamer(R__b);
635 if (iChamber->Nsec()==1) {
636 segmentation=iChamber->GetSegmentationModel(1);
637 segmentation->Streamer(R__b);
638 } else {
639 segmentation=iChamber->GetSegmentationModel(1);
640 segmentation->Streamer(R__b);
641 segmentation=iChamber->GetSegmentationModel(2);
642 segmentation->Streamer(R__b);
643 }
644 response=iChamber->GetResponseModel(mip);
645 response->Streamer(R__b);
646 response=iChamber->GetResponseModel(cerenkov);
647 response->Streamer(R__b);
648
649 digitsaddress=(TClonesArray*) (*fDchambers)[i];
650 digitsaddress->Streamer(R__b);
651 }
fe4da5cc 652 }
ddae0931 653}
654AliRICHcluster* AliRICH::FirstPad(AliRICHhit* hit,TClonesArray *clusters )
655{
656//
657 // Initialise the pad iterator
658 // Return the address of the first padhit for hit
659 TClonesArray *theClusters = clusters;
660 Int_t nclust = theClusters->GetEntriesFast();
661 if (nclust && hit->fPHlast > 0) {
662 sMaxIterPad=hit->fPHlast;
663 sCurIterPad=hit->fPHfirst;
664 return (AliRICHcluster*) clusters->UncheckedAt(sCurIterPad-1);
665 } else {
666 return 0;
fe4da5cc 667 }
ddae0931 668
fe4da5cc 669}
670
ddae0931 671AliRICHcluster* AliRICH::NextPad(TClonesArray *clusters)
fe4da5cc 672{
ddae0931 673 sCurIterPad++;
674 if (sCurIterPad <= sMaxIterPad) {
675 return (AliRICHcluster*) clusters->UncheckedAt(sCurIterPad-1);
676 } else {
677 return 0;
678 }
fe4da5cc 679}
680
ddae0931 681ClassImp(AliRICHcluster)
682
683//___________________________________________
684AliRICHcluster::AliRICHcluster(Int_t *clhits)
685{
686 fHitNumber=clhits[0];
687 fCathode=clhits[1];
688 fQ=clhits[2];
689 fPadX=clhits[3];
690 fPadY=clhits[4];
691 fQpad=clhits[5];
692 fRSec=clhits[6];
693}
694ClassImp(AliRICHdigit)
fe4da5cc 695//_____________________________________________________________________________
ddae0931 696AliRICHdigit::AliRICHdigit(Int_t *digits)
fe4da5cc 697{
ddae0931 698 //
699 // Creates a RICH digit object to be updated
700 //
701 fPadX = digits[0];
702 fPadY = digits[1];
703 fSignal = digits[2];
704
fe4da5cc 705}
fe4da5cc 706//_____________________________________________________________________________
ddae0931 707AliRICHdigit::AliRICHdigit(Int_t *tracks, Int_t *charges, Int_t *digits)
fe4da5cc 708{
ddae0931 709 //
710 // Creates a RICH digit object
711 //
712 fPadX = digits[0];
713 fPadY = digits[1];
714 fSignal = digits[2];
715 for(Int_t i=0; i<10; i++) {
716 fTcharges[i] = charges[i];
717 fTracks[i] = tracks[i];
718 }
fe4da5cc 719}
720
ddae0931 721ClassImp(AliRICHlist)
722
723//____________________________________________________________________________
724AliRICHlist::AliRICHlist(Int_t ich, Int_t *digits):
725 AliRICHdigit(digits)
fe4da5cc 726{
ddae0931 727 //
728 // Creates a RICH digit list object
729 //
730
731 fChamber = ich;
732 fTrackList = new TObjArray;
733
fe4da5cc 734}
ddae0931 735//_____________________________________________________________________________
fe4da5cc 736
fe4da5cc 737
ddae0931 738ClassImp(AliRICHhit)
739
740//___________________________________________
741AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
742 AliHit(shunt, track)
fe4da5cc 743{
ddae0931 744 fChamber=vol[0];
745 fParticle=hits[0];
746 fX=hits[1];
747 fY=hits[2];
748 fZ=hits[3];
749 fTheta=hits[4];
750 fPhi=hits[5];
751 fTlength=hits[6];
752 fEloss=hits[7];
753 fPHfirst=(Int_t) hits[8];
754 fPHlast=(Int_t) hits[9];
fe4da5cc 755}
ddae0931 756ClassImp(AliRICHreccluster)
757
758ClassImp(AliRICHCerenkov)
759//___________________________________________
760AliRICHCerenkov::AliRICHCerenkov(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
761 AliHit(shunt, track)
fe4da5cc 762{
ddae0931 763 fChamber=vol[0];
764 fX=hits[1];
765 fY=hits[2];
766 fZ=hits[3];
767 fTheta=hits[4];
768 fPhi=hits[5];
769 fTlength=hits[6];
770 fPHfirst=(Int_t) hits[8];
771 fPHlast=(Int_t) hits[9];
fe4da5cc 772}
773
ddae0931 774ClassImp(AliRICHRecCluster)
775
776//_____________________________________________________________________
777AliRICHRecCluster::AliRICHRecCluster()
fe4da5cc 778{
ddae0931 779 fDigits=0;
780 fNdigit=-1;
fe4da5cc 781}
ddae0931 782
783AliRICHRecCluster::AliRICHRecCluster(Int_t FirstDigit,Int_t Ichamber, Int_t Icathod)
fe4da5cc 784{
ddae0931 785 fX = 0.;
786 fY = 0.;
787 fDigits = new TArrayI(10);
788 fNdigit=0;
789 AddDigit(FirstDigit);
790 fChamber=Ichamber;
791 fCathod=Icathod;
fe4da5cc 792}
793
ddae0931 794void AliRICHRecCluster::AddDigit(Int_t Digit)
fe4da5cc 795{
ddae0931 796 if (fNdigit==fDigits->GetSize()) {
797 //enlarge the list by hand!
798 Int_t *array= new Int_t[fNdigit*2];
799 for(Int_t i=0;i<fNdigit;i++)
800 array[i] = fDigits->At(i);
801 fDigits->Adopt(fNdigit*2,array);
802 }
803 fDigits->AddAt(Digit,fNdigit);
804 fNdigit++;
fe4da5cc 805}
806
ddae0931 807
808AliRICHRecCluster::~AliRICHRecCluster()
fe4da5cc 809{
ddae0931 810 if (fDigits)
811 delete fDigits;
fe4da5cc 812}
813
ddae0931 814Int_t AliRICHRecCluster::FirstDigitIndex()
fe4da5cc 815{
ddae0931 816 fCurrentDigit=0;
817 return fDigits->At(fCurrentDigit);
fe4da5cc 818}
819
ddae0931 820Int_t AliRICHRecCluster::NextDigitIndex()
821{
822 fCurrentDigit++;
823 if (fCurrentDigit<fNdigit)
824 return fDigits->At(fCurrentDigit);
825 else
826 return InvalidDigitIndex();
827}
fe4da5cc 828
ddae0931 829Int_t AliRICHRecCluster::NDigits()
fe4da5cc 830{
ddae0931 831 return fNdigit;
fe4da5cc 832}
833
ddae0931 834void AliRICHRecCluster::Finish()
fe4da5cc 835{
ddae0931 836 // In order to reconstruct coordinates, one has to
837 // get back to the digits which is not trivial here,
838 // because we don't know where digits are stored!
839 // Center Of Gravity, or other method should be
840 // a property of AliRICH class!
fe4da5cc 841}
842
fe4da5cc 843
844
ddae0931 845void AliRICH::Digitise(Int_t nev,Option_t *option,Text_t *filename)
fe4da5cc 846{
ddae0931 847 // keep galice.root for signal and name differently the file for
848 // background when add! otherwise the track info for signal will be lost !
849
c90dd3e2 850 static Bool_t first=kTRUE;
ddae0931 851 static TTree *TH1;
852 static TFile *File;
c90dd3e2 853 Int_t i;
ddae0931 854 char *Add = strstr(option,"Add");
855
856 AliRICHchamber* iChamber;
857 AliRICHsegmentation* segmentation;
858
859
860 Int_t trk[50];
861 Int_t chtrk[50];
862 TObjArray *list=new TObjArray;
863 Int_t digits[3];
864
865 AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
866 AliRICHHitMap* HitMap[10];
c90dd3e2 867 for (i=0; i<10; i++) {HitMap[i]=0;}
ddae0931 868 if (Add ) {
869 if(first) {
870 fFileName=filename;
871 cout<<"filename"<<fFileName<<endl;
872 File=new TFile(fFileName);
873 cout<<"I have opened "<<fFileName<<" file "<<endl;
874 fHits2 = new TClonesArray("AliRICHhit",1000 );
875 fClusters2 = new TClonesArray("AliRICHcluster",10000);
c90dd3e2 876 first=kFALSE;
ddae0931 877 }
878 File->cd();
879 File->ls();
880 // Get Hits Tree header from file
881 if(fHits2) fHits2->Clear();
882 if(fClusters2) fClusters2->Clear();
883 if(TH1) delete TH1;
884 TH1=0;
885 //
886 char treeName[20];
887 sprintf(treeName,"TreeH%d",nev);
888 TH1 = (TTree*)gDirectory->Get(treeName);
889 if (!TH1) {
890 printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
891 }
892 // Set branch addresses
893 TBranch *branch;
894 char branchname[20];
895 sprintf(branchname,"%s",GetName());
896 if (TH1 && fHits2) {
897 branch = TH1->GetBranch(branchname);
898 if (branch) branch->SetAddress(&fHits2);
899 }
900 if (TH1 && fClusters2) {
901 branch = TH1->GetBranch("RICHCluster");
902 if (branch) branch->SetAddress(&fClusters2);
903 }
904 }
905 //
906 // loop over cathodes
907 //
908 AliRICHHitMap* hm;
909
f91473f6 910 for (int icat=0; icat<1; icat++) {
c90dd3e2 911 for (i=0; i<7; i++) {
ddae0931 912 if (HitMap[i]) {
913 hm=HitMap[i];
914 delete hm;
915 HitMap[i]=0;
916 }
917 }
918 Int_t counter=0;
c90dd3e2 919 for (i =0; i<7; i++) {
ddae0931 920 iChamber=(AliRICHchamber*) (*fChambers)[i];
921 if (iChamber->Nsec()==1 && icat==1) {
922 continue;
923 } else {
924 segmentation=iChamber->GetSegmentationModel(icat+1);
925 }
926 HitMap[i] = new AliRICHHitMapA1(segmentation, list);
927 }
928 printf("Start loop over tracks \n");
929//
930// Loop over tracks
931//
932
933 TTree *TH = gAlice->TreeH();
934 Int_t ntracks =(Int_t) TH->GetEntries();
935 for (Int_t track=0; track<ntracks; track++) {
936 gAlice->ResetHits();
937 TH->GetEvent(track);
938//
939// Loop over hits
940 for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1);
941 mHit;
942 mHit=(AliRICHhit*)RICH->NextHit())
943 {
944 Int_t nch = mHit->fChamber-1; // chamber number
945 if (nch >7) continue;
946 iChamber = &(RICH->Chamber(nch));
947
948//
949// Loop over pad hits
950 for (AliRICHcluster* mPad=
951 (AliRICHcluster*)RICH->FirstPad(mHit,fClusters);
952 mPad;
953 mPad=(AliRICHcluster*)RICH->NextPad(fClusters))
954 {
955 Int_t cathode = mPad->fCathode; // cathode number
956 Int_t ipx = mPad->fPadX; // pad number on X
957 Int_t ipy = mPad->fPadY; // pad number on Y
958 Int_t iqpad = mPad->fQpad; // charge per pad
959//
960//
961
962 if (cathode != (icat+1)) continue;
963 // fill the info array
964 Float_t thex, they;
965 segmentation=iChamber->GetSegmentationModel(cathode);
966 segmentation->GetPadCxy(ipx,ipy,thex,they);
967 TVector *trinfo_p= new TVector(2);
968 TVector &trinfo = *trinfo_p;
969 trinfo(0)=(Float_t)track;
970 trinfo(1)=(Float_t)iqpad;
971
972 digits[0]=ipx;
973 digits[1]=ipy;
974 digits[2]=iqpad;
975
976 AliRICHlist* pdigit;
977 // build the list of fired pads and update the info
978 if (!HitMap[nch]->TestHit(ipx, ipy)) {
979 list->AddAtAndExpand(
980 new AliRICHlist(nch,digits),counter);
981 HitMap[nch]->SetHit(ipx, ipy, counter);
982 counter++;
983 pdigit=(AliRICHlist*)list->At(list->GetLast());
984 // list of tracks
985 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
986 trlist->Add(&trinfo);
987 } else {
988 pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
989 // update charge
990 (*pdigit).fSignal+=iqpad;
991 // update list of tracks
992 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
993 Int_t last_entry=trlist->GetLast();
994 TVector *ptrk_p=(TVector*)trlist->At(last_entry);
995 TVector &ptrk=*ptrk_p;
996 Int_t last_track=Int_t(ptrk(0));
997 Int_t last_charge=Int_t(ptrk(1));
998 if (last_track==track) {
999 last_charge+=iqpad;
1000 trlist->RemoveAt(last_entry);
1001 trinfo(0)=last_track;
1002 trinfo(1)=last_charge;
1003 trlist->AddAt(&trinfo,last_entry);
1004 } else {
1005 trlist->Add(&trinfo);
1006 }
1007 // check the track list
1008 Int_t nptracks=trlist->GetEntriesFast();
1009 if (nptracks > 2) {
1010 printf("Attention - nptracks > 2 %d \n",nptracks);
1011 printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
1012 for (Int_t tr=0;tr<nptracks;tr++) {
1013 TVector *pptrk_p=(TVector*)trlist->At(tr);
1014 TVector &pptrk=*pptrk_p;
1015 trk[tr]=Int_t(pptrk(0));
1016 chtrk[tr]=Int_t(pptrk(1));
1017 }
1018 } // end if nptracks
1019 } // end if pdigit
1020 } //end loop over clusters
1021 } // hit loop
1022 } // track loop
1023
1024 Int_t nentr1=list->GetEntriesFast();
1025 printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
1026
1027 // open the file with background
1028
1029 if (Add ) {
1030 ntracks =(Int_t)TH1->GetEntries();
1031 printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
1032 printf("background - Start loop over tracks \n");
1033//
1034// Loop over tracks
1035//
f91473f6 1036 for (Int_t trak=0; trak<ntracks; trak++) {
ddae0931 1037
1038 if (fHits2) fHits2->Clear();
1039 if (fClusters2) fClusters2->Clear();
1040
f91473f6 1041 TH1->GetEvent(trak);
ddae0931 1042//
1043// Loop over hits
1044 AliRICHhit* mHit;
f91473f6 1045 for(int j=0;j<fHits2->GetEntriesFast();++j)
ddae0931 1046 {
f91473f6 1047 mHit=(AliRICHhit*) (*fHits2)[j];
ddae0931 1048 Int_t nch = mHit->fChamber-1; // chamber number
1049 if (nch >9) continue;
1050 iChamber = &(RICH->Chamber(nch));
1051 Int_t rmin = (Int_t)iChamber->RInner();
1052 Int_t rmax = (Int_t)iChamber->ROuter();
1053//
1054// Loop over pad hits
1055 for (AliRICHcluster* mPad=
1056 (AliRICHcluster*)RICH->FirstPad(mHit,fClusters2);
1057 mPad;
1058 mPad=(AliRICHcluster*)RICH->NextPad(fClusters2))
1059 {
1060 Int_t cathode = mPad->fCathode; // cathode number
1061 Int_t ipx = mPad->fPadX; // pad number on X
1062 Int_t ipy = mPad->fPadY; // pad number on Y
1063 Int_t iqpad = mPad->fQpad; // charge per pad
f91473f6 1064 if (trak==3 && nch==0 && icat==0) printf("bgr - trak,iqpad,ipx,ipy %d %d %d %d\n",trak,iqpad,ipx,ipy);
ddae0931 1065//
1066//
1067 if (cathode != (icat+1)) continue;
1068 // fill the info array
1069 Float_t thex, they;
1070 segmentation=iChamber->GetSegmentationModel(cathode);
1071 segmentation->GetPadCxy(ipx,ipy,thex,they);
1072 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1073 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1074
1075 TVector *trinfo_p;
1076 trinfo_p = new TVector(2);
1077 TVector &trinfo = *trinfo_p;
1078 trinfo(0)=-1; // tag background
1079 trinfo(1)=-1;
1080
1081 digits[0]=ipx;
1082 digits[1]=ipy;
1083 digits[2]=iqpad;
1084
1085
f91473f6 1086 if (trak <4 && icat==0 && nch==0)
1087 printf("bgr - HitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
1088 HitMap[nch]->TestHit(ipx, ipy),trak);
ddae0931 1089 AliRICHlist* pdigit;
1090 // build the list of fired pads and update the info
1091 if (!HitMap[nch]->TestHit(ipx, ipy)) {
1092 list->AddAtAndExpand(new AliRICHlist(nch,digits),counter);
1093
1094 HitMap[nch]->SetHit(ipx, ipy, counter);
1095 counter++;
1096 printf("bgr new elem in list - counter %d\n",counter);
1097
1098 pdigit=(AliRICHlist*)list->At(list->GetLast());
1099 // list of tracks
1100 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
1101 trlist->Add(&trinfo);
1102 } else {
1103 pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
1104 // update charge
1105 (*pdigit).fSignal+=iqpad;
1106 // update list of tracks
1107 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
1108 Int_t last_entry=trlist->GetLast();
1109 TVector *ptrk_p=(TVector*)trlist->At(last_entry);
1110 TVector &ptrk=*ptrk_p;
1111 Int_t last_track=Int_t(ptrk(0));
1112 if (last_track==-1) {
1113 continue;
1114 } else {
1115 trlist->Add(&trinfo);
1116 }
1117 // check the track list
1118 Int_t nptracks=trlist->GetEntriesFast();
1119 if (nptracks > 0) {
1120 for (Int_t tr=0;tr<nptracks;tr++) {
1121 TVector *pptrk_p=(TVector*)trlist->At(tr);
1122 TVector &pptrk=*pptrk_p;
1123 trk[tr]=Int_t(pptrk(0));
1124 chtrk[tr]=Int_t(pptrk(1));
1125 }
1126 } // end if nptracks
1127 } // end if pdigit
1128 } //end loop over clusters
1129 } // hit loop
1130 } // track loop
1131 Int_t nentr2=list->GetEntriesFast();
1132 printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1133 TTree *fAli=gAlice->TreeK();
1134 if (fAli) File =fAli->GetCurrentFile();
1135 File->cd();
1136 } // if Add
1137
1138 Int_t tracks[10];
1139 Int_t charges[10];
1140 cout<<"start filling digits \n "<<endl;
1141 Int_t nentries=list->GetEntriesFast();
1142 printf(" \n \n nentries %d \n",nentries);
1143
1144 // start filling the digits
1145
1146 for (Int_t nent=0;nent<nentries;nent++) {
1147 AliRICHlist *address=(AliRICHlist*)list->At(nent);
1148 if (address==0) continue;
1149 Int_t ich=address->fChamber;
1150 Int_t q=address->fSignal;
1151 iChamber=(AliRICHchamber*) (*fChambers)[ich];
1152 // add white noise and do zero-suppression and signal truncation
1153 Float_t MeanNoise = gRandom->Gaus(1, 0.2);
1154 Float_t ZeroSupp=5*MeanNoise;
1155 Float_t Noise = gRandom->Gaus(0, MeanNoise);
1156 q+=(Int_t)Noise;
1157 if ( q <= ZeroSupp) continue;
1158 digits[0]=address->fPadX;
1159 digits[1]=address->fPadY;
1160 digits[2]=q;
1161
1162 TObjArray* trlist=(TObjArray*)address->TrackList();
1163 Int_t nptracks=trlist->GetEntriesFast();
1164
1165 // this was changed to accomodate the real number of tracks
1166 if (nptracks > 10) {
1167 cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1168 nptracks=10;
1169 }
1170 if (nptracks > 2) {
1171 printf("Attention - nptracks > 2 %d \n",nptracks);
1172 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1173 }
1174 for (Int_t tr=0;tr<nptracks;tr++) {
1175 TVector *pp_p=(TVector*)trlist->At(tr);
1176 TVector &pp =*pp_p;
1177 tracks[tr]=Int_t(pp(0));
1178 charges[tr]=Int_t(pp(1));
1179 } //end loop over list of tracks for one pad
1180 if (nptracks < 10 ) {
f91473f6 1181 for (Int_t t=nptracks; t<10; t++) {
1182 tracks[t]=0;
1183 charges[t]=0;
ddae0931 1184 }
1185 }
1186 // fill digits
1187 RICH->AddDigits(ich,tracks,charges,digits);
1188
1189 delete address;
1190 }
1191 cout<<"I'm out of the loops for digitisation"<<endl;
1192 gAlice->TreeD()->Fill();
1193 TTree *TD=gAlice->TreeD();
1194 Stat_t ndig=TD->GetEntries();
1195 cout<<"number of digits "<<ndig<<endl;
1196 TClonesArray *fDch;
f91473f6 1197 for (int k=0;k<7;k++) {
1198 fDch= RICH->DigitsAddress(k);
c90dd3e2 1199 int ndigit=fDch->GetEntriesFast();
1200 printf (" k, ndigits %d %d \n",k,ndigit);
ddae0931 1201 }
1202 RICH->ResetDigits();
1203
1204 list->Clear();
1205
1206 } //end loop over cathodes
1207 char hname[30];
1208 sprintf(hname,"TreeD%d",nev);
1209 gAlice->TreeD()->Write(hname);
fe4da5cc 1210}
1211
1212
ddae0931 1213
1214
1215
1216
1217
1218
1219