Added kDebugLevel variable to control output size on demand
[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/*
2e5f0f7b 17 $Log$
33984590 18 Revision 1.18 2000/06/12 15:15:46 jbarbosa
19 Cleaned up version.
20
237c933d 21 Revision 1.17 2000/06/09 14:58:37 jbarbosa
22 New digitisation per particle type
23
d28b5fc2 24 Revision 1.16 2000/04/19 12:55:43 morsch
25 Newly structured and updated version (JB, AM)
26
4c039060 27*/
28
2e5f0f7b 29
ddae0931 30////////////////////////////////////////////////
31// Manager and hits classes for set:RICH //
32////////////////////////////////////////////////
fe4da5cc 33
ddae0931 34#include <TBRIK.h>
35#include <TTUBE.h>
fe4da5cc 36#include <TNode.h>
37#include <TRandom.h>
ddae0931 38#include <TObject.h>
39#include <TVector.h>
40#include <TObjArray.h>
2e5f0f7b 41#include <TArrayF.h>
237c933d 42#include <TFile.h>
43#include <TParticle.h>
44#include <iostream.h>
fe4da5cc 45
fe4da5cc 46#include "AliRICH.h"
237c933d 47#include "AliRICHSegmentation.h"
48#include "AliRICHHit.h"
49#include "AliRICHCerenkov.h"
50#include "AliRICHPadHit.h"
51#include "AliRICHDigit.h"
52#include "AliRICHTransientDigit.h"
53#include "AliRICHRawCluster.h"
54#include "AliRICHRecHit.h"
55#include "AliRICHHitMapA1.h"
2e5f0f7b 56#include "AliRICHClusterFinder.h"
fe4da5cc 57#include "AliRun.h"
ddae0931 58#include "AliMC.h"
2e5f0f7b 59#include "AliPoints.h"
ddae0931 60#include "AliCallf77.h"
237c933d 61
ddae0931 62
63// Static variables for the pad-hit iterator routines
64static Int_t sMaxIterPad=0;
65static Int_t sCurIterPad=0;
66static TClonesArray *fClusters2;
67static TClonesArray *fHits2;
2e5f0f7b 68static TTree *TrH1;
ddae0931 69
fe4da5cc 70ClassImp(AliRICH)
ddae0931 71
72//___________________________________________
fe4da5cc 73AliRICH::AliRICH()
74{
237c933d 75// Default constructor for RICH manager class
76
ddae0931 77 fIshunt = 0;
78 fHits = 0;
2e5f0f7b 79 fPadHits = 0;
80 fNPadHits = 0;
81 fNcerenkovs = 0;
ddae0931 82 fDchambers = 0;
ddae0931 83 fCerenkovs = 0;
84 fNdch = 0;
fe4da5cc 85}
ddae0931 86
87//___________________________________________
fe4da5cc 88AliRICH::AliRICH(const char *name, const char *title)
ddae0931 89 : AliDetector(name,title)
fe4da5cc 90{
ddae0931 91//Begin_Html
92/*
93 <img src="gif/alirich.gif">
94*/
95//End_Html
96
2e5f0f7b 97 fHits = new TClonesArray("AliRICHHit",1000 );
1cedd08a 98 gAlice->AddHitList(fHits);
2e5f0f7b 99 fPadHits = new TClonesArray("AliRICHPadHit",100000);
100 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
101 gAlice->AddHitList(fCerenkovs);
102 //gAlice->AddHitList(fHits);
103 fNPadHits = 0;
104 fNcerenkovs = 0;
105 fIshunt = 0;
ddae0931 106
237c933d 107 fNdch = new Int_t[kNCH];
ddae0931 108
237c933d 109 fDchambers = new TObjArray(kNCH);
2e5f0f7b 110
237c933d 111 fRecHits = new TObjArray(kNCH);
ddae0931 112
113 Int_t i;
114
237c933d 115 for (i=0; i<kNCH ;i++) {
2e5f0f7b 116 (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
ddae0931 117 fNdch[i]=0;
118 }
2e5f0f7b 119
237c933d 120 fNrawch = new Int_t[kNCH];
ddae0931 121
237c933d 122 fRawClusters = new TObjArray(kNCH);
d28b5fc2 123 //printf("Created fRwClusters with adress:%p",fRawClusters);
2e5f0f7b 124
237c933d 125 for (i=0; i<kNCH ;i++) {
2e5f0f7b 126 (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
127 fNrawch[i]=0;
128 }
129
237c933d 130 fNrechits = new Int_t[kNCH];
ddae0931 131
237c933d 132 for (i=0; i<kNCH ;i++) {
2e5f0f7b 133 (*fRecHits)[i] = new TClonesArray("AliRICHRecHit",1000);
134 }
d28b5fc2 135 //printf("Created fRecHits with adress:%p",fRecHits);
2e5f0f7b 136
137
ddae0931 138 SetMarkerColor(kRed);
fe4da5cc 139}
140
237c933d 141AliRICH::AliRICH(const AliRICH& RICH)
142{
143// Copy Constructor
144}
145
146
ddae0931 147//___________________________________________
fe4da5cc 148AliRICH::~AliRICH()
149{
237c933d 150
151// Destructor of RICH manager class
152
ddae0931 153 fIshunt = 0;
154 delete fHits;
2e5f0f7b 155 delete fPadHits;
ddae0931 156 delete fCerenkovs;
fe4da5cc 157}
158
ddae0931 159//___________________________________________
fe4da5cc 160void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
161{
237c933d 162
163//
164// Adds a hit to the Hits list
165//
166
ddae0931 167 TClonesArray &lhits = *fHits;
2e5f0f7b 168 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
fe4da5cc 169}
fe4da5cc 170//_____________________________________________________________________________
ddae0931 171void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
172{
237c933d 173
174//
175// Adds a RICH cerenkov hit to the Cerenkov Hits list
176//
177
ddae0931 178 TClonesArray &lcerenkovs = *fCerenkovs;
179 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
2e5f0f7b 180 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
fe4da5cc 181}
ddae0931 182//___________________________________________
2e5f0f7b 183void AliRICH::AddPadHit(Int_t *clhits)
ddae0931 184{
237c933d 185
186//
187// Add a RICH pad hit to the list
188//
189
2e5f0f7b 190 TClonesArray &lPadHits = *fPadHits;
191 new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
33984590 192}
fe4da5cc 193//_____________________________________________________________________________
ddae0931 194void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
195{
237c933d 196
197 //
198 // Add a RICH digit to the list
199 //
ddae0931 200
201 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
2e5f0f7b 202 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
fe4da5cc 203}
204
205//_____________________________________________________________________________
2e5f0f7b 206void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
fe4da5cc 207{
ddae0931 208 //
2e5f0f7b 209 // Add a RICH digit to the list
ddae0931 210 //
2e5f0f7b 211
212 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
213 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
fe4da5cc 214}
215
2e5f0f7b 216//_____________________________________________________________________________
33984590 217void AliRICH::AddRecHit(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
2e5f0f7b 218{
237c933d 219
220 //
221 // Add a RICH reconstructed hit to the list
222 //
fe4da5cc 223
2e5f0f7b 224 TClonesArray &lrec = *((TClonesArray*)(*fRecHits)[id]);
33984590 225 new(lrec[fNrechits[id]++]) AliRICHRecHit(id,rechit,photons,padsx,padsy);
2e5f0f7b 226}
fe4da5cc 227
ddae0931 228//___________________________________________
229void AliRICH::BuildGeometry()
230
fe4da5cc 231{
237c933d 232
233 //
234 // Builds a TNode geometry for event display
235 //
236 TNode *node, *top;
ddae0931 237
238 const int kColorRICH = kGreen;
239 //
237c933d 240 top=gAlice->GetGeometry()->GetNode("alice");
ddae0931 241
242
243 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
244
237c933d 245 top->cd();
ddae0931 246 Float_t pos1[3]={0,471.8999,165.2599};
2e5f0f7b 247 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
248 new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
237c933d 249 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
ddae0931 250
251
237c933d 252 node->SetLineColor(kColorRICH);
253 fNodes->Add(node);
254 top->cd();
ddae0931 255
256 Float_t pos2[3]={171,470,0};
2e5f0f7b 257 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
258 new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
237c933d 259 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
ddae0931 260
261
237c933d 262 node->SetLineColor(kColorRICH);
263 fNodes->Add(node);
264 top->cd();
ddae0931 265 Float_t pos3[3]={0,500,0};
2e5f0f7b 266 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
267 new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
237c933d 268 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
ddae0931 269
270
237c933d 271 node->SetLineColor(kColorRICH);
272 fNodes->Add(node);
273 top->cd();
ddae0931 274 Float_t pos4[3]={-171,470,0};
2e5f0f7b 275 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
276 new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
237c933d 277 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
ddae0931 278
279
237c933d 280 node->SetLineColor(kColorRICH);
281 fNodes->Add(node);
282 top->cd();
ddae0931 283 Float_t pos5[3]={161.3999,443.3999,-165.3};
2e5f0f7b 284 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
285 new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
237c933d 286 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
ddae0931 287
237c933d 288 node->SetLineColor(kColorRICH);
289 fNodes->Add(node);
290 top->cd();
ddae0931 291 Float_t pos6[3]={0., 471.9, -165.3,};
2e5f0f7b 292 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
293 new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
237c933d 294 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
ddae0931 295
296
237c933d 297 node->SetLineColor(kColorRICH);
298 fNodes->Add(node);
299 top->cd();
ddae0931 300 Float_t pos7[3]={-161.399,443.3999,-165.3};
2e5f0f7b 301 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
302 new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
237c933d 303 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
304 node->SetLineColor(kColorRICH);
305 fNodes->Add(node);
ddae0931 306
fe4da5cc 307}
308
ddae0931 309//___________________________________________
310Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
311{
237c933d 312
313// Default value
314
ddae0931 315 return 9999;
316}
fe4da5cc 317
ddae0931 318//___________________________________________
319void AliRICH::MakeBranch(Option_t* option)
fe4da5cc 320{
237c933d 321 // Create Tree branches for the RICH.
fe4da5cc 322
237c933d 323 const Int_t kBufferSize = 4000;
ddae0931 324 char branchname[20];
fe4da5cc 325
ddae0931 326
327 AliDetector::MakeBranch(option);
328 sprintf(branchname,"%sCerenkov",GetName());
329 if (fCerenkovs && gAlice->TreeH()) {
237c933d 330 gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize);
ddae0931 331 printf("Making Branch %s for Cerenkov Hits\n",branchname);
fe4da5cc 332 }
ddae0931 333
2e5f0f7b 334 sprintf(branchname,"%sPadHits",GetName());
335 if (fPadHits && gAlice->TreeH()) {
237c933d 336 gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
2e5f0f7b 337 printf("Making Branch %s for PadHits\n",branchname);
fe4da5cc 338 }
339
ddae0931 340// one branch for digits per chamber
341 Int_t i;
342
237c933d 343 for (i=0; i<kNCH ;i++) {
ddae0931 344 sprintf(branchname,"%sDigits%d",GetName(),i+1);
345
346 if (fDchambers && gAlice->TreeD()) {
237c933d 347 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
ddae0931 348 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
349 }
fe4da5cc 350 }
2e5f0f7b 351
352// one branch for raw clusters per chamber
237c933d 353 for (i=0; i<kNCH ;i++) {
2e5f0f7b 354 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
355
356 if (fRawClusters && gAlice->TreeR()) {
237c933d 357 gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
2e5f0f7b 358 printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
359 }
360 }
361
362 // one branch for rec hits per chamber
237c933d 363 for (i=0; i<kNCH ;i++) {
2e5f0f7b 364 sprintf(branchname,"%sRecHits%d",GetName(),i+1);
365
366 if (fRecHits && gAlice->TreeR()) {
237c933d 367 gAlice->TreeR()->Branch(branchname,&((*fRecHits)[i]), kBufferSize);
2e5f0f7b 368 printf("Making Branch %s for rec. hits in chamber %d\n",branchname,i+1);
369 }
370 }
fe4da5cc 371}
372
ddae0931 373//___________________________________________
374void AliRICH::SetTreeAddress()
fe4da5cc 375{
237c933d 376 // Set branch address for the Hits and Digits Tree.
2e5f0f7b 377 char branchname[20];
378 Int_t i;
379
ddae0931 380 AliDetector::SetTreeAddress();
381
382 TBranch *branch;
383 TTree *treeH = gAlice->TreeH();
384 TTree *treeD = gAlice->TreeD();
2e5f0f7b 385 TTree *treeR = gAlice->TreeR();
ddae0931 386
387 if (treeH) {
2e5f0f7b 388 if (fPadHits) {
389 branch = treeH->GetBranch("RICHPadHits");
390 if (branch) branch->SetAddress(&fPadHits);
fe4da5cc 391 }
ddae0931 392 if (fCerenkovs) {
393 branch = treeH->GetBranch("RICHCerenkov");
394 if (branch) branch->SetAddress(&fCerenkovs);
fe4da5cc 395 }
ddae0931 396 }
397
398 if (treeD) {
237c933d 399 for (int i=0; i<kNCH; i++) {
ddae0931 400 sprintf(branchname,"%sDigits%d",GetName(),i+1);
401 if (fDchambers) {
402 branch = treeD->GetBranch(branchname);
403 if (branch) branch->SetAddress(&((*fDchambers)[i]));
404 }
fe4da5cc 405 }
fe4da5cc 406 }
2e5f0f7b 407 if (treeR) {
237c933d 408 for (i=0; i<kNCH; i++) {
2e5f0f7b 409 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
410 if (fRawClusters) {
411 branch = treeR->GetBranch(branchname);
412 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
413 }
414 }
415
237c933d 416 for (i=0; i<kNCH; i++) {
2e5f0f7b 417 sprintf(branchname,"%sRecHits%d",GetName(),i+1);
418 if (fRecHits) {
419 branch = treeR->GetBranch(branchname);
420 if (branch) branch->SetAddress(&((*fRecHits)[i]));
421 }
422 }
423
424 }
ddae0931 425}
426//___________________________________________
427void AliRICH::ResetHits()
428{
237c933d 429 // Reset number of clusters and the cluster array for this detector
ddae0931 430 AliDetector::ResetHits();
2e5f0f7b 431 fNPadHits = 0;
432 fNcerenkovs = 0;
433 if (fPadHits) fPadHits->Clear();
ddae0931 434 if (fCerenkovs) fCerenkovs->Clear();
fe4da5cc 435}
436
2e5f0f7b 437
ddae0931 438//____________________________________________
439void AliRICH::ResetDigits()
440{
237c933d 441 //
442 // Reset number of digits and the digits array for this detector
443 //
444 for ( int i=0;i<kNCH;i++ ) {
ddae0931 445 if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
446 if (fNdch) fNdch[i]=0;
447 }
448}
2e5f0f7b 449
ddae0931 450//____________________________________________
2e5f0f7b 451void AliRICH::ResetRawClusters()
fe4da5cc 452{
237c933d 453 //
454 // Reset number of raw clusters and the raw clust array for this detector
455 //
456 for ( int i=0;i<kNCH;i++ ) {
2e5f0f7b 457 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
458 if (fNrawch) fNrawch[i]=0;
fe4da5cc 459 }
ddae0931 460}
fe4da5cc 461
2e5f0f7b 462//____________________________________________
463void AliRICH::ResetRecHits()
fe4da5cc 464{
237c933d 465 //
466 // Reset number of raw clusters and the raw clust array for this detector
467 //
2e5f0f7b 468
237c933d 469 for ( int i=0;i<kNCH;i++ ) {
2e5f0f7b 470 if ((*fRecHits)[i]) ((TClonesArray*)(*fRecHits)[i])->Clear();
471 if (fNrechits) fNrechits[i]=0;
472 }
fe4da5cc 473}
474
ddae0931 475//___________________________________________
2e5f0f7b 476void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
fe4da5cc 477{
237c933d 478
479//
480// Setter for the RICH geometry model
481//
482
483
2e5f0f7b 484 ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
fe4da5cc 485}
486
ddae0931 487//___________________________________________
2e5f0f7b 488void AliRICH::SetSegmentationModel(Int_t id, AliRICHSegmentation *segmentation)
ddae0931 489{
237c933d 490
491//
492// Setter for the RICH segmentation model
493//
494
2e5f0f7b 495 ((AliRICHChamber*) (*fChambers)[id])->SegmentationModel(segmentation);
ddae0931 496}
fe4da5cc 497
ddae0931 498//___________________________________________
2e5f0f7b 499void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
fe4da5cc 500{
237c933d 501
502//
503// Setter for the RICH response model
504//
505
2e5f0f7b 506 ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
fe4da5cc 507}
508
2e5f0f7b 509void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
ddae0931 510{
237c933d 511
512//
513// Setter for the RICH reconstruction model (clusters)
514//
515
2e5f0f7b 516 ((AliRICHChamber*) (*fChambers)[id])->ReconstructionModel(reconst);
fe4da5cc 517}
518
ddae0931 519void AliRICH::SetNsec(Int_t id, Int_t nsec)
520{
237c933d 521
522//
523// Sets the number of padplanes
524//
525
2e5f0f7b 526 ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
ddae0931 527}
528
529
530//___________________________________________
531
532void AliRICH::StepManager()
fe4da5cc 533{
237c933d 534
535// Dummy step manager (should never be called)
536
ddae0931 537}
2e5f0f7b 538
237c933d 539void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
ddae0931 540{
2e5f0f7b 541
ddae0931 542//
543// Loop on chambers and on cathode planes
544//
2e5f0f7b 545 for (Int_t icat=1;icat<2;icat++) {
546 gAlice->ResetDigits();
547 gAlice->TreeD()->GetEvent(1); // spurious +1 ...
237c933d 548 for (Int_t ich=0;ich<kNCH;ich++) {
2e5f0f7b 549 AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
237c933d 550 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
551 if (pRICHdigits == 0)
2e5f0f7b 552 continue;
553 //
554 // Get ready the current chamber stuff
555 //
556 AliRICHResponse* response = iChamber->GetResponseModel();
557 AliRICHSegmentation* seg = iChamber->GetSegmentationModel();
558 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
559 if (seg) {
560 rec->SetSegmentation(seg);
561 rec->SetResponse(response);
237c933d 562 rec->SetDigits(pRICHdigits);
2e5f0f7b 563 rec->SetChamber(ich);
564 if (nev==0) rec->CalibrateCOG();
565 rec->FindRawClusters();
566 }
567 TClonesArray *fRch;
568 fRch=RawClustAddress(ich);
569 fRch->Sort();
570 } // for ich
571
572 gAlice->TreeR()->Fill();
573 TClonesArray *fRch;
237c933d 574 for (int i=0;i<kNCH;i++) {
2e5f0f7b 575 fRch=RawClustAddress(i);
576 int nraw=fRch->GetEntriesFast();
577 printf ("Chamber %d, raw clusters %d\n",i,nraw);
578 }
579
580 ResetRawClusters();
581
582 } // for icat
583
584 char hname[30];
585 sprintf(hname,"TreeR%d",nev);
33984590 586 gAlice->TreeR()->Write(hname,kOverwrite,0);
2e5f0f7b 587 gAlice->TreeR()->Reset();
588
589 //gObjectTable->Print();
ddae0931 590}
591
592
593//______________________________________________________________________________
594void AliRICH::Streamer(TBuffer &R__b)
595{
596 // Stream an object of class AliRICH.
2e5f0f7b 597 AliRICHChamber *iChamber;
598 AliRICHSegmentation *segmentation;
599 AliRICHResponse *response;
ddae0931 600 TClonesArray *digitsaddress;
2e5f0f7b 601 TClonesArray *rawcladdress;
602 TClonesArray *rechitaddress;
ddae0931 603
604 if (R__b.IsReading()) {
605 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
606 AliDetector::Streamer(R__b);
2e5f0f7b 607 R__b >> fNPadHits;
608 R__b >> fPadHits; // diff
609 R__b >> fNcerenkovs;
610 R__b >> fCerenkovs; // diff
ddae0931 611 R__b >> fDchambers;
2e5f0f7b 612 R__b >> fRawClusters;
33984590 613 R__b >> fRecHits; //diff
614 R__b >> fDebugLevel; //diff
ddae0931 615 R__b.ReadArray(fNdch);
2e5f0f7b 616 R__b.ReadArray(fNrawch);
617 R__b.ReadArray(fNrechits);
618//
ddae0931 619 R__b >> fChambers;
620// Stream chamber related information
237c933d 621 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 622 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 623 iChamber->Streamer(R__b);
2e5f0f7b 624 segmentation=iChamber->GetSegmentationModel();
625 segmentation->Streamer(R__b);
626 response=iChamber->GetResponseModel();
ddae0931 627 response->Streamer(R__b);
2e5f0f7b 628 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
629 rawcladdress->Streamer(R__b);
630 rechitaddress=(TClonesArray*) (*fRecHits)[i];
631 rechitaddress->Streamer(R__b);
ddae0931 632 digitsaddress=(TClonesArray*) (*fDchambers)[i];
633 digitsaddress->Streamer(R__b);
634 }
635
636 } else {
637 R__b.WriteVersion(AliRICH::IsA());
638 AliDetector::Streamer(R__b);
2e5f0f7b 639 R__b << fNPadHits;
640 R__b << fPadHits; // diff
641 R__b << fNcerenkovs;
642 R__b << fCerenkovs; // diff
ddae0931 643 R__b << fDchambers;
2e5f0f7b 644 R__b << fRawClusters;
645 R__b << fRecHits; //diff
33984590 646 R__b << fDebugLevel; //diff
237c933d 647 R__b.WriteArray(fNdch, kNCH);
648 R__b.WriteArray(fNrawch, kNCH);
649 R__b.WriteArray(fNrechits, kNCH);
ddae0931 650 //
ddae0931 651 R__b << fChambers;
652// Stream chamber related information
237c933d 653 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 654 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 655 iChamber->Streamer(R__b);
2e5f0f7b 656 segmentation=iChamber->GetSegmentationModel();
657 segmentation->Streamer(R__b);
658 response=iChamber->GetResponseModel();
ddae0931 659 response->Streamer(R__b);
2e5f0f7b 660 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
661 rawcladdress->Streamer(R__b);
662 rechitaddress=(TClonesArray*) (*fRecHits)[i];
663 rechitaddress->Streamer(R__b);
ddae0931 664 digitsaddress=(TClonesArray*) (*fDchambers)[i];
665 digitsaddress->Streamer(R__b);
666 }
fe4da5cc 667 }
ddae0931 668}
2e5f0f7b 669AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
ddae0931 670{
671//
672 // Initialise the pad iterator
673 // Return the address of the first padhit for hit
674 TClonesArray *theClusters = clusters;
675 Int_t nclust = theClusters->GetEntriesFast();
676 if (nclust && hit->fPHlast > 0) {
2e5f0f7b 677 sMaxIterPad=Int_t(hit->fPHlast);
678 sCurIterPad=Int_t(hit->fPHfirst);
679 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 680 } else {
681 return 0;
fe4da5cc 682 }
ddae0931 683
fe4da5cc 684}
685
2e5f0f7b 686AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters)
fe4da5cc 687{
237c933d 688
689 // Iterates over pads
690
ddae0931 691 sCurIterPad++;
692 if (sCurIterPad <= sMaxIterPad) {
2e5f0f7b 693 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 694 } else {
695 return 0;
696 }
fe4da5cc 697}
698
fe4da5cc 699
d28b5fc2 700void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
fe4da5cc 701{
ddae0931 702 // keep galice.root for signal and name differently the file for
703 // background when add! otherwise the track info for signal will be lost !
2e5f0f7b 704
c90dd3e2 705 static Bool_t first=kTRUE;
237c933d 706 static TFile *pFile;
707 char *addBackground = strstr(option,"Add");
ddae0931 708
2e5f0f7b 709 FILE* points; //these will be the digits...
710
711 points=fopen("points.dat","w");
712
713 AliRICHChamber* iChamber;
714 AliRICHSegmentation* segmentation;
ddae0931 715
d28b5fc2 716 Int_t digitse=0;
ddae0931 717 Int_t trk[50];
718 Int_t chtrk[50];
719 TObjArray *list=new TObjArray;
237c933d 720 static TClonesArray *pAddress=0;
721 if(!pAddress) pAddress=new TClonesArray("TVector",1000);
2e5f0f7b 722 Int_t digits[5];
ddae0931 723
237c933d 724 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
725 AliRICHHitMap* pHitMap[10];
2e5f0f7b 726 Int_t i;
237c933d 727 for (i=0; i<10; i++) {pHitMap[i]=0;}
728 if (addBackground ) {
ddae0931 729 if(first) {
730 fFileName=filename;
731 cout<<"filename"<<fFileName<<endl;
237c933d 732 pFile=new TFile(fFileName);
ddae0931 733 cout<<"I have opened "<<fFileName<<" file "<<endl;
2e5f0f7b 734 fHits2 = new TClonesArray("AliRICHHit",1000 );
735 fClusters2 = new TClonesArray("AliRICHPadHit",10000);
c90dd3e2 736 first=kFALSE;
ddae0931 737 }
237c933d 738 pFile->cd();
ddae0931 739 // Get Hits Tree header from file
740 if(fHits2) fHits2->Clear();
741 if(fClusters2) fClusters2->Clear();
2e5f0f7b 742 if(TrH1) delete TrH1;
743 TrH1=0;
744
ddae0931 745 char treeName[20];
746 sprintf(treeName,"TreeH%d",nev);
2e5f0f7b 747 TrH1 = (TTree*)gDirectory->Get(treeName);
748 if (!TrH1) {
ddae0931 749 printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
750 }
751 // Set branch addresses
752 TBranch *branch;
753 char branchname[20];
754 sprintf(branchname,"%s",GetName());
2e5f0f7b 755 if (TrH1 && fHits2) {
756 branch = TrH1->GetBranch(branchname);
ddae0931 757 if (branch) branch->SetAddress(&fHits2);
758 }
2e5f0f7b 759 if (TrH1 && fClusters2) {
760 branch = TrH1->GetBranch("RICHCluster");
ddae0931 761 if (branch) branch->SetAddress(&fClusters2);
762 }
763 }
33984590 764
ddae0931 765 AliRICHHitMap* hm;
2e5f0f7b 766 Int_t countadr=0;
33984590 767 Int_t counter=0;
768 for (i =0; i<kNCH; i++) {
769 iChamber=(AliRICHChamber*) (*fChambers)[i];
770 segmentation=iChamber->GetSegmentationModel(1);
771 pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
772 }
773 //
774 // Loop over tracks
775 //
776
777 TTree *treeH = gAlice->TreeH();
778 Int_t ntracks =(Int_t) treeH->GetEntries();
779 for (Int_t track=0; track<ntracks; track++) {
780 gAlice->ResetHits();
781 treeH->GetEvent(track);
782 //
783 // Loop over hits
784 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
785 mHit;
786 mHit=(AliRICHHit*)pRICH->NextHit())
787 {
788
789 digitse=0;
790
791 Int_t nch = mHit->fChamber-1; // chamber number
792 if (nch >kNCH) continue;
793 iChamber = &(pRICH->Chamber(nch));
794
795 TParticle *current = (TParticle*)(*gAlice->Particles())[track];
796
797 Int_t particle = current->GetPdgCode();
798
799 //printf("Flag:%d\n",flag);
800 //printf("Track:%d\n",track);
801 //printf("Particle:%d\n",particle);
802
803 if (flag == 0)
804 digitse=1;
805
806 if (flag == 1)
807 if(TMath::Abs(particle) == 211 || TMath::Abs(particle) == 111)
808 digitse=1;
809
810 if (flag == 2)
811 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
812 || TMath::Abs(particle)==311)
813 digitse=1;
814
815 if (flag == 3 && TMath::Abs(particle)==2212)
816 digitse=1;
817
818 if (flag == 4 && TMath::Abs(particle)==13)
819 digitse=1;
820
821 if (flag == 5 && TMath::Abs(particle)==11)
822 digitse=1;
823
824 if (flag == 6 && TMath::Abs(particle)==2112)
825 digitse=1;
826
827
828 //printf ("Particle: %d, Flag: %d, Digitse: %d\n",particle,flag,digitse);
829
830
831 if (digitse)
ddae0931 832 {
d28b5fc2 833
33984590 834 //
835 // Loop over pad hits
836 for (AliRICHPadHit* mPad=
837 (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
838 mPad;
839 mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
d28b5fc2 840 {
33984590 841 Int_t cathode = mPad->fCathode; // cathode number
842 Int_t ipx = mPad->fPadX; // pad number on X
843 Int_t ipy = mPad->fPadY; // pad number on Y
844 Int_t iqpad = mPad->fQpad; // charge per pad
845 //
846 //
847 //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
848
849 Float_t thex, they;
850 segmentation=iChamber->GetSegmentationModel(cathode);
851 segmentation->GetPadCxy(ipx,ipy,thex,they);
852 new((*pAddress)[countadr++]) TVector(2);
853 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
854 trinfo(0)=(Float_t)track;
855 trinfo(1)=(Float_t)iqpad;
856
857 digits[0]=ipx;
858 digits[1]=ipy;
859 digits[2]=iqpad;
860
861 AliRICHTransientDigit* pdigit;
862 // build the list of fired pads and update the info
863 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
864 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
865 pHitMap[nch]->SetHit(ipx, ipy, counter);
866 counter++;
867 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
868 // list of tracks
869 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
870 trlist->Add(&trinfo);
871 } else {
872 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
873 // update charge
874 (*pdigit).fSignal+=iqpad;
875 // update list of tracks
876 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
877 Int_t lastEntry=trlist->GetLast();
878 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
879 TVector &ptrk=*ptrkP;
880 Int_t lastTrack=Int_t(ptrk(0));
881 Int_t lastCharge=Int_t(ptrk(1));
882 if (lastTrack==track) {
883 lastCharge+=iqpad;
884 trlist->RemoveAt(lastEntry);
885 trinfo(0)=lastTrack;
886 trinfo(1)=lastCharge;
887 trlist->AddAt(&trinfo,lastEntry);
888 } else {
889 trlist->Add(&trinfo);
890 }
891 // check the track list
892 Int_t nptracks=trlist->GetEntriesFast();
893 if (nptracks > 2) {
894 printf("Attention - tracks: %d (>2)\n",nptracks);
895 //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
896 for (Int_t tr=0;tr<nptracks;tr++) {
897 TVector *pptrkP=(TVector*)trlist->At(tr);
898 TVector &pptrk=*pptrkP;
899 trk[tr]=Int_t(pptrk(0));
900 chtrk[tr]=Int_t(pptrk(1));
901 }
902 } // end if nptracks
903 } // end if pdigit
904 } //end loop over clusters
905 }// track type condition
906 } // hit loop
907 } // track loop
908
909 // open the file with background
910
911 if (addBackground ) {
912 ntracks =(Int_t)TrH1->GetEntries();
913 //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
914 //printf("background - Start loop over tracks \n");
915 //
916 // Loop over tracks
917 //
918 for (Int_t trak=0; trak<ntracks; trak++) {
919 if (fHits2) fHits2->Clear();
920 if (fClusters2) fClusters2->Clear();
921 TrH1->GetEvent(trak);
922 //
923 // Loop over hits
924 AliRICHHit* mHit;
925 for(int j=0;j<fHits2->GetEntriesFast();++j)
926 {
927 mHit=(AliRICHHit*) (*fHits2)[j];
928 Int_t nch = mHit->fChamber-1; // chamber number
929 if (nch >6) continue;
930 iChamber = &(pRICH->Chamber(nch));
931 Int_t rmin = (Int_t)iChamber->RInner();
932 Int_t rmax = (Int_t)iChamber->ROuter();
933 //
934 // Loop over pad hits
935 for (AliRICHPadHit* mPad=
936 (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
937 mPad;
938 mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2))
939 {
940 Int_t cathode = mPad->fCathode; // cathode number
941 Int_t ipx = mPad->fPadX; // pad number on X
942 Int_t ipy = mPad->fPadY; // pad number on Y
943 Int_t iqpad = mPad->fQpad; // charge per pad
2e5f0f7b 944
33984590 945 Float_t thex, they;
946 segmentation=iChamber->GetSegmentationModel(cathode);
947 segmentation->GetPadCxy(ipx,ipy,thex,they);
948 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
949 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
950 new((*pAddress)[countadr++]) TVector(2);
951 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
952 trinfo(0)=-1; // tag background
953 trinfo(1)=-1;
954 digits[0]=ipx;
955 digits[1]=ipy;
956 digits[2]=iqpad;
957 if (trak <4 && nch==0)
958 printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
959 pHitMap[nch]->TestHit(ipx, ipy),trak);
960 AliRICHTransientDigit* pdigit;
961 // build the list of fired pads and update the info
962 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
963 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
964
965 pHitMap[nch]->SetHit(ipx, ipy, counter);
966 counter++;
967 printf("bgr new elem in list - counter %d\n",counter);
968
969 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
970 // list of tracks
971 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
972 trlist->Add(&trinfo);
973 } else {
974 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
975 // update charge
976 (*pdigit).fSignal+=iqpad;
977 // update list of tracks
978 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
979 Int_t lastEntry=trlist->GetLast();
980 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
981 TVector &ptrk=*ptrkP;
982 Int_t lastTrack=Int_t(ptrk(0));
983 if (lastTrack==-1) {
984 continue;
985 } else {
986 trlist->Add(&trinfo);
987 }
988 // check the track list
989 Int_t nptracks=trlist->GetEntriesFast();
990 if (nptracks > 0) {
991 for (Int_t tr=0;tr<nptracks;tr++) {
992 TVector *pptrkP=(TVector*)trlist->At(tr);
993 TVector &pptrk=*pptrkP;
994 trk[tr]=Int_t(pptrk(0));
995 chtrk[tr]=Int_t(pptrk(1));
996 }
997 } // end if nptracks
998 } // end if pdigit
999 } //end loop over clusters
1000 } // hit loop
1001 } // track loop
ddae0931 1002 TTree *fAli=gAlice->TreeK();
237c933d 1003 if (fAli) pFile =fAli->GetCurrentFile();
1004 pFile->cd();
33984590 1005 } // if Add
1006
1007 Int_t tracks[10];
1008 Int_t charges[10];
1009 //cout<<"Start filling digits \n "<<endl;
1010 Int_t nentries=list->GetEntriesFast();
1011 //printf(" \n \n nentries %d \n",nentries);
1012
1013 // start filling the digits
1014
1015 for (Int_t nent=0;nent<nentries;nent++) {
1016 AliRICHTransientDigit *address=(AliRICHTransientDigit*)list->At(nent);
1017 if (address==0) continue;
1018
1019 Int_t ich=address->fChamber;
1020 Int_t q=address->fSignal;
1021 iChamber=(AliRICHChamber*) (*fChambers)[ich];
1022 AliRICHResponse * response=iChamber->GetResponseModel();
1023 Int_t adcmax= (Int_t) response->MaxAdc();
1024
1025
1026 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
1027 //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY));
1028 Int_t pedestal = iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY);
ddae0931 1029
33984590 1030 //printf("Pedestal:%d\n",pedestal);
1031 //Int_t pedestal=0;
1032 Float_t treshold = (pedestal + 4*1.7);
1033
1034 Float_t meanNoise = gRandom->Gaus(1.7, 0.25);
1035 Float_t noise = gRandom->Gaus(0, meanNoise);
1036 q+=(Int_t)(noise + pedestal);
1037 //q+=(Int_t)(noise);
1038 // magic number to be parametrised !!!
1039 if ( q <= treshold)
1040 {
1041 q = q - pedestal;
1042 continue;
2e5f0f7b 1043 }
33984590 1044 q = q - pedestal;
1045 if ( q >= adcmax) q=adcmax;
1046 digits[0]=address->fPadX;
1047 digits[1]=address->fPadY;
1048 digits[2]=q;
1049
1050 TObjArray* trlist=(TObjArray*)address->TrackList();
1051 Int_t nptracks=trlist->GetEntriesFast();
1052
1053 // this was changed to accomodate the real number of tracks
1054 if (nptracks > 10) {
1055 cout<<"Attention - tracks > 10 "<<nptracks<<endl;
1056 nptracks=10;
1057 }
1058 if (nptracks > 2) {
1059 printf("Attention - tracks > 2 %d \n",nptracks);
1060 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
1061 //icat,ich,digits[0],digits[1],q);
1062 }
1063 for (Int_t tr=0;tr<nptracks;tr++) {
1064 TVector *ppP=(TVector*)trlist->At(tr);
1065 TVector &pp =*ppP;
1066 tracks[tr]=Int_t(pp(0));
1067 charges[tr]=Int_t(pp(1));
1068 } //end loop over list of tracks for one pad
1069 if (nptracks < 10 ) {
1070 for (Int_t t=nptracks; t<10; t++) {
1071 tracks[t]=0;
1072 charges[t]=0;
ddae0931 1073 }
33984590 1074 }
1075 //write file
1076 if (ich==2)
1077 fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
1078
1079 // fill digits
1080 pRICH->AddDigits(ich,tracks,charges,digits);
1081 }
1082 gAlice->TreeD()->Fill();
1083
1084 list->Delete();
1085 for(Int_t ii=0;ii<kNCH;++ii) {
1086 if (pHitMap[ii]) {
1087 hm=pHitMap[ii];
1088 delete hm;
1089 pHitMap[ii]=0;
1090 }
1091 }
1092
1093 //TTree *TD=gAlice->TreeD();
1094 //Stat_t ndig=TD->GetEntries();
1095 //cout<<"number of digits "<<ndig<<endl;
1096 TClonesArray *fDch;
1097 for (int k=0;k<kNCH;k++) {
1098 fDch= pRICH->DigitsAddress(k);
1099 int ndigit=fDch->GetEntriesFast();
1100 printf ("Chamber %d digits %d \n",k,ndigit);
1101 }
1102 pRICH->ResetDigits();
ddae0931 1103 char hname[30];
1104 sprintf(hname,"TreeD%d",nev);
33984590 1105 gAlice->TreeD()->Write(hname,kOverwrite,0);
1106
1107 // reset tree
1108 // gAlice->TreeD()->Reset();
2e5f0f7b 1109 delete list;
237c933d 1110 pAddress->Clear();
33984590 1111 // gObjectTable->Print();
2e5f0f7b 1112}
ddae0931 1113
237c933d 1114AliRICH& AliRICH::operator=(const AliRICH& rhs)
1115{
1116// Assignment operator
1117 return *this;
1118
1119}
ddae0931 1120
237c933d 1121Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
1122{
1123//
1124// Calls the charge disintegration method of the current chamber and adds
1125// the simulated cluster to the root treee
1126//
1127 Int_t clhits[kNCH];
1128 Float_t newclust[6][500];
1129 Int_t nnew;
1130
1131//
1132// Integrated pulse height on chamber
1133
1134 clhits[0]=fNhits+1;
1135
1136 ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
1137 Int_t ic=0;
1138
1139//
1140// Add new clusters
1141 for (Int_t i=0; i<nnew; i++) {
1142 if (Int_t(newclust[3][i]) > 0) {
1143 ic++;
1144// Cathode plane
1145 clhits[1] = Int_t(newclust[5][i]);
1146// Cluster Charge
1147 clhits[2] = Int_t(newclust[0][i]);
1148// Pad: ix
1149 clhits[3] = Int_t(newclust[1][i]);
1150// Pad: iy
1151 clhits[4] = Int_t(newclust[2][i]);
1152// Pad: charge
1153 clhits[5] = Int_t(newclust[3][i]);
1154// Pad: chamber sector
1155 clhits[6] = Int_t(newclust[4][i]);
1156
1157 AddPadHit(clhits);
1158 }
1159 }
1160return nnew;
1161}
ddae0931 1162