a4451541474aa57c0944736ed0575fa00a050d08
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
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   Revision 1.17  2000/06/09 14:58:37  jbarbosa
19   New digitisation per particle type
20
21   Revision 1.16  2000/04/19 12:55:43  morsch
22   Newly structured and updated version (JB, AM)
23
24 */
25
26
27 ////////////////////////////////////////////////
28 //  Manager and hits classes for set:RICH     //
29 ////////////////////////////////////////////////
30
31 #include <TBRIK.h>
32 #include <TTUBE.h>
33 #include <TNode.h> 
34 #include <TRandom.h> 
35 #include <TObject.h>
36 #include <TVector.h>
37 #include <TObjArray.h>
38 #include <TArrayF.h>
39 #include <TFile.h>
40 #include <TParticle.h>
41 #include <iostream.h>
42
43 #include "AliRICH.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"
54 #include "AliRun.h"
55 #include "AliMC.h"
56 #include "AliPoints.h"
57 #include "AliCallf77.h" 
58
59
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;
65 static TTree *TrH1;
66  
67 ClassImp(AliRICH)
68     
69 //___________________________________________
70 AliRICH::AliRICH()
71 {
72 // Default constructor for RICH manager class
73
74     fIshunt     = 0;
75     fHits       = 0;
76     fPadHits    = 0;
77     fNPadHits   = 0;
78     fNcerenkovs = 0;
79     fDchambers  = 0;
80     fCerenkovs  = 0;
81     fNdch       = 0;
82 }
83
84 //___________________________________________
85 AliRICH::AliRICH(const char *name, const char *title)
86     : AliDetector(name,title)
87 {
88 //Begin_Html
89 /*
90   <img src="gif/alirich.gif">
91 */
92 //End_Html
93     
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);
100     fNPadHits   = 0;
101     fNcerenkovs = 0;
102     fIshunt     = 0;
103     
104     fNdch      = new Int_t[kNCH];
105     
106     fDchambers = new TObjArray(kNCH);
107
108     fRecHits = new TObjArray(kNCH);
109     
110     Int_t i;
111    
112     for (i=0; i<kNCH ;i++) {
113         (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000); 
114         fNdch[i]=0;
115     }
116
117     fNrawch      = new Int_t[kNCH];
118     
119     fRawClusters = new TObjArray(kNCH);
120     //printf("Created fRwClusters with adress:%p",fRawClusters);
121
122     for (i=0; i<kNCH ;i++) {
123       (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000); 
124       fNrawch[i]=0;
125     }
126
127     fNrechits      = new Int_t[kNCH];
128     
129     for (i=0; i<kNCH ;i++) {
130         (*fRecHits)[i] = new TClonesArray("AliRICHRecHit",1000); 
131     }
132     //printf("Created fRecHits with adress:%p",fRecHits);
133
134         
135     SetMarkerColor(kRed);
136 }
137
138 AliRICH::AliRICH(const AliRICH& RICH)
139 {
140 // Copy Constructor
141 }
142
143
144 //___________________________________________
145 AliRICH::~AliRICH()
146 {
147
148 // Destructor of RICH manager class
149
150     fIshunt  = 0;
151     delete fHits;
152     delete fPadHits;
153     delete fCerenkovs;
154 }
155
156 //___________________________________________
157 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
158 {
159
160 //  
161 // Adds a hit to the Hits list
162 //
163
164     TClonesArray &lhits = *fHits;
165     new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
166 }
167 //_____________________________________________________________________________
168 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
169 {
170
171 //
172 // Adds a RICH cerenkov hit to the Cerenkov Hits list
173 //
174
175     TClonesArray &lcerenkovs = *fCerenkovs;
176     new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
177     //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
178 }
179 //___________________________________________
180 void AliRICH::AddPadHit(Int_t *clhits)
181 {
182
183 //
184 // Add a RICH pad hit to the list
185 //
186
187     TClonesArray &lPadHits = *fPadHits;
188     new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
189 }
190 //_____________________________________________________________________________
191 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
192 {
193
194   //
195   // Add a RICH digit to the list
196   //
197
198     TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
199     new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
200 }
201
202 //_____________________________________________________________________________
203 void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
204 {
205     //
206     // Add a RICH digit to the list
207     //
208
209     TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
210     new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
211 }
212
213 //_____________________________________________________________________________
214 void AliRICH::AddRecHit(Int_t id, Float_t *rechit)
215 {
216   
217   //
218   // Add a RICH reconstructed hit to the list
219   //
220
221     TClonesArray &lrec = *((TClonesArray*)(*fRecHits)[id]);
222     new(lrec[fNrechits[id]++]) AliRICHRecHit(id,rechit);
223 }
224
225 //___________________________________________
226 void AliRICH::BuildGeometry()
227     
228 {
229   
230   //
231   // Builds a TNode geometry for event display
232   //
233     TNode *node, *top;
234     
235     const int kColorRICH = kGreen;
236     //
237     top=gAlice->GetGeometry()->GetNode("alice");
238     
239     
240     new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
241     
242     top->cd();
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");
247     
248
249     node->SetLineColor(kColorRICH);
250     fNodes->Add(node);
251     top->cd();
252     
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");
257     
258     
259     node->SetLineColor(kColorRICH);
260     fNodes->Add(node);
261     top->cd();
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");
266     
267
268     node->SetLineColor(kColorRICH);
269     fNodes->Add(node);
270     top->cd();
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");
275     
276
277     node->SetLineColor(kColorRICH);
278     fNodes->Add(node);
279     top->cd();
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");
284     
285     node->SetLineColor(kColorRICH);
286     fNodes->Add(node);
287     top->cd();
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");
292     
293     
294     node->SetLineColor(kColorRICH);
295     fNodes->Add(node);
296     top->cd();
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);
302     fNodes->Add(node); 
303     
304 }
305
306 //___________________________________________
307 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
308 {
309
310 // Default value
311
312     return 9999;
313 }
314
315 //___________________________________________
316 void AliRICH::MakeBranch(Option_t* option)
317 {
318   // Create Tree branches for the RICH.
319     
320     const Int_t kBufferSize = 4000;
321     char branchname[20];
322     
323     
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);
329     }
330     
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);
335     }
336     
337 // one branch for digits per chamber
338     Int_t i;
339     
340     for (i=0; i<kNCH ;i++) {
341         sprintf(branchname,"%sDigits%d",GetName(),i+1);
342         
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);
346         }       
347     }
348
349 // one branch for raw clusters per chamber
350   for (i=0; i<kNCH ;i++) {
351       sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
352       
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);
356       } 
357   }
358
359   // one branch for rec hits per chamber
360   for (i=0; i<kNCH ;i++) {
361     sprintf(branchname,"%sRecHits%d",GetName(),i+1);
362     
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);
366     }   
367   }
368 }
369
370 //___________________________________________
371 void AliRICH::SetTreeAddress()
372 {
373   // Set branch address for the Hits and Digits Tree.
374   char branchname[20];
375   Int_t i;
376
377     AliDetector::SetTreeAddress();
378     
379     TBranch *branch;
380     TTree *treeH = gAlice->TreeH();
381     TTree *treeD = gAlice->TreeD();
382     TTree *treeR = gAlice->TreeR();
383     
384     if (treeH) {
385         if (fPadHits) {
386             branch = treeH->GetBranch("RICHPadHits");
387             if (branch) branch->SetAddress(&fPadHits);
388         }
389         if (fCerenkovs) {
390             branch = treeH->GetBranch("RICHCerenkov");
391             if (branch) branch->SetAddress(&fCerenkovs);
392         }
393     }
394     
395     if (treeD) {
396         for (int i=0; i<kNCH; i++) {
397             sprintf(branchname,"%sDigits%d",GetName(),i+1);
398             if (fDchambers) {
399                 branch = treeD->GetBranch(branchname);
400                 if (branch) branch->SetAddress(&((*fDchambers)[i]));
401             }
402         }
403     }
404   if (treeR) {
405       for (i=0; i<kNCH; i++) {
406           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
407           if (fRawClusters) {
408               branch = treeR->GetBranch(branchname);
409               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
410           }
411       }
412       
413       for (i=0; i<kNCH; i++) {
414         sprintf(branchname,"%sRecHits%d",GetName(),i+1);
415         if (fRecHits) {
416           branch = treeR->GetBranch(branchname);
417           if (branch) branch->SetAddress(&((*fRecHits)[i]));
418           }
419       }
420       
421   }
422 }
423 //___________________________________________
424 void AliRICH::ResetHits()
425 {
426   // Reset number of clusters and the cluster array for this detector
427     AliDetector::ResetHits();
428     fNPadHits   = 0;
429     fNcerenkovs = 0;
430     if (fPadHits)  fPadHits->Clear();
431     if (fCerenkovs) fCerenkovs->Clear();
432 }
433
434
435 //____________________________________________
436 void AliRICH::ResetDigits()
437 {
438   //
439   // Reset number of digits and the digits array for this detector
440   //
441     for ( int i=0;i<kNCH;i++ ) {
442         if ((*fDchambers)[i])   (*fDchambers)[i]->Clear();
443         if (fNdch)  fNdch[i]=0;
444     }
445 }
446
447 //____________________________________________
448 void AliRICH::ResetRawClusters()
449 {
450   //
451   // Reset number of raw clusters and the raw clust array for this detector
452   //
453     for ( int i=0;i<kNCH;i++ ) {
454         if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
455         if (fNrawch)  fNrawch[i]=0;
456     }
457 }
458
459 //____________________________________________
460 void AliRICH::ResetRecHits()
461 {
462   //
463   // Reset number of raw clusters and the raw clust array for this detector
464   //
465   
466   for ( int i=0;i<kNCH;i++ ) {
467         if ((*fRecHits)[i])    ((TClonesArray*)(*fRecHits)[i])->Clear();
468         if (fNrechits)  fNrechits[i]=0;
469     }
470 }
471
472 //___________________________________________
473 void   AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
474 {
475
476 //
477 // Setter for the RICH geometry model
478 //
479
480
481     ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
482 }
483
484 //___________________________________________
485 void   AliRICH::SetSegmentationModel(Int_t id, AliRICHSegmentation *segmentation)
486 {
487
488 //
489 // Setter for the RICH segmentation model
490 //
491
492     ((AliRICHChamber*) (*fChambers)[id])->SegmentationModel(segmentation);
493 }
494
495 //___________________________________________
496 void   AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
497 {
498
499 //
500 // Setter for the RICH response model
501 //
502
503     ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
504 }
505
506 void   AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
507 {
508
509 //
510 // Setter for the RICH reconstruction model (clusters)
511 //
512
513     ((AliRICHChamber*) (*fChambers)[id])->ReconstructionModel(reconst);
514 }
515
516 void   AliRICH::SetNsec(Int_t id, Int_t nsec)
517 {
518
519 //
520 // Sets the number of padplanes
521 //
522
523     ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
524 }
525
526
527 //___________________________________________
528
529 void AliRICH::StepManager()
530 {
531
532 // Dummy step manager (should never be called)
533
534 }
535
536 void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
537 {
538
539 //
540 // Loop on chambers and on cathode planes
541 //
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)       
549               continue;
550           //
551           // Get ready the current chamber stuff
552           //
553           AliRICHResponse* response = iChamber->GetResponseModel();
554           AliRICHSegmentation*  seg = iChamber->GetSegmentationModel();
555           AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
556           if (seg) {      
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();
563           }  
564           TClonesArray *fRch;
565           fRch=RawClustAddress(ich);
566           fRch->Sort();
567         } // for ich
568
569         gAlice->TreeR()->Fill();
570         TClonesArray *fRch;
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);
575         }
576         
577         ResetRawClusters();
578         
579     } // for icat
580     
581     char hname[30];
582     sprintf(hname,"TreeR%d",nev);
583     gAlice->TreeR()->Write(hname);
584     gAlice->TreeR()->Reset();
585     
586     //gObjectTable->Print();
587 }
588
589
590 //______________________________________________________________________________
591 void AliRICH::Streamer(TBuffer &R__b)
592 {
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;
600       
601     if (R__b.IsReading()) {
602         Version_t R__v = R__b.ReadVersion(); if (R__v) { }
603         AliDetector::Streamer(R__b);
604         R__b >> fNPadHits;
605         R__b >> fPadHits;   // diff
606         R__b >> fNcerenkovs;
607         R__b >> fCerenkovs; // diff
608         R__b >> fDchambers;
609         R__b >> fRawClusters;
610         R__b >> fRecHits;  //diff       
611         R__b.ReadArray(fNdch);
612         R__b.ReadArray(fNrawch);
613         R__b.ReadArray(fNrechits);
614 //
615         R__b >> fChambers;
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);
630         }
631       
632     } else {
633         R__b.WriteVersion(AliRICH::IsA());
634         AliDetector::Streamer(R__b);
635         R__b << fNPadHits;
636         R__b << fPadHits; // diff
637         R__b << fNcerenkovs;
638         R__b << fCerenkovs; // diff
639         R__b << fDchambers;
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);
645         //
646         R__b << fChambers;
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);
661         }
662     }
663 }
664 AliRICHPadHit* AliRICH::FirstPad(AliRICHHit*  hit,TClonesArray *clusters ) 
665 {
666 //
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);
675     } else {
676         return 0;
677     }
678     
679 }
680
681 AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters) 
682 {
683
684   // Iterates over pads
685   
686     sCurIterPad++;
687     if (sCurIterPad <= sMaxIterPad) {
688         return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
689     } else {
690         return 0;
691     }
692 }
693
694
695 void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
696 {
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 !
699
700     static Bool_t first=kTRUE;
701     static TFile *pFile;
702     char *addBackground = strstr(option,"Add");
703
704     FILE* points; //these will be the digits...
705
706     points=fopen("points.dat","w");
707
708     AliRICHChamber*       iChamber;
709     AliRICHSegmentation*  segmentation;
710
711     Int_t digitse=0;
712     Int_t trk[50];
713     Int_t chtrk[50];  
714     TObjArray *list=new TObjArray;
715     static TClonesArray *pAddress=0;
716     if(!pAddress) pAddress=new TClonesArray("TVector",1000);
717     Int_t digits[5]; 
718     
719     AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
720     AliRICHHitMap* pHitMap[10];
721     Int_t i;
722     for (i=0; i<10; i++) {pHitMap[i]=0;}
723     if (addBackground ) {
724         if(first) {
725             fFileName=filename;
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);
731             first=kFALSE;
732         }
733         pFile->cd();
734         // Get Hits Tree header from file
735         if(fHits2) fHits2->Clear();
736         if(fClusters2) fClusters2->Clear();
737         if(TrH1) delete TrH1;
738         TrH1=0;
739
740         char treeName[20];
741         sprintf(treeName,"TreeH%d",nev);
742         TrH1 = (TTree*)gDirectory->Get(treeName);
743         if (!TrH1) {
744             printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
745         }
746         // Set branch addresses
747         TBranch *branch;
748         char branchname[20];
749         sprintf(branchname,"%s",GetName());
750         if (TrH1 && fHits2) {
751             branch = TrH1->GetBranch(branchname);
752             if (branch) branch->SetAddress(&fHits2);
753         }
754         if (TrH1 && fClusters2) {
755             branch = TrH1->GetBranch("RICHCluster");
756             if (branch) branch->SetAddress(&fClusters2);
757         }
758     }
759     //
760     // loop over cathodes
761     //
762     AliRICHHitMap* hm;
763     Int_t countadr=0;
764     for (int icat=0; icat<1; icat++) { 
765         Int_t counter=0;
766         for (i =0; i<kNCH; i++) {
767             iChamber=(AliRICHChamber*) (*fChambers)[i];
768             if (iChamber->Nsec()==1 && icat==1) {
769                 continue;
770             } else {
771                 segmentation=iChamber->GetSegmentationModel(icat+1);
772             }
773             pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
774         }
775 //
776 //   Loop over tracks
777 //
778
779         TTree *treeH = gAlice->TreeH();
780         Int_t ntracks =(Int_t) treeH->GetEntries();
781         for (Int_t track=0; track<ntracks; track++) {
782             gAlice->ResetHits();
783             treeH->GetEvent(track);
784 //
785 //   Loop over hits
786             for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1); 
787                 mHit;
788                 mHit=(AliRICHHit*)pRICH->NextHit()) 
789             {
790               
791               digitse=0;
792               
793               Int_t   nch   = mHit->fChamber-1;  // chamber number
794               if (nch >kNCH) continue;
795               iChamber = &(pRICH->Chamber(nch));
796               
797               TParticle *current = (TParticle*)(*gAlice->Particles())[track];
798               
799               Int_t particle = current->GetPdgCode();
800               
801               //printf("Flag:%d\n",flag);
802               //printf("Track:%d\n",track);
803               //printf("Particle:%d\n",particle);
804               
805               if (flag == 0)
806                 digitse=1;
807               
808               if (flag == 1) 
809                 if(TMath::Abs(particle) == 211 || TMath::Abs(particle) == 111)
810                   digitse=1;
811               
812               if (flag == 2)
813                 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
814                    || TMath::Abs(particle)==311)
815                   digitse=1;
816               
817               if (flag == 3 && TMath::Abs(particle)==2212)
818                 digitse=1;
819               
820               if (flag == 4 && TMath::Abs(particle)==13)
821                 digitse=1;
822
823               if (flag == 5 && TMath::Abs(particle)==11)
824                 digitse=1;
825               
826               if (flag == 6 && TMath::Abs(particle)==2112)
827                 digitse=1;
828                   
829
830               //printf ("Particle: %d, Flag: %d, Digitse: %d\n",particle,flag,digitse); 
831
832
833               if (digitse)
834                 {
835                 
836 //
837 // Loop over pad hits
838                   for (AliRICHPadHit* mPad=
839                          (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
840                        mPad;
841                        mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
842                     {
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
847                       //
848                       //
849                       //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
850                       
851                       Float_t thex, they;
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;
858                       
859                       digits[0]=ipx;
860                       digits[1]=ipy;
861                       digits[2]=iqpad;
862                       
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);
868                         counter++;
869                         pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
870                         // list of tracks
871                         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
872                         trlist->Add(&trinfo);
873                       } else {
874                         pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
875                         // update charge
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) {
885                             lastCharge+=iqpad;
886                             trlist->RemoveAt(lastEntry);
887                             trinfo(0)=lastTrack;
888                             trinfo(1)=lastCharge;
889                             trlist->AddAt(&trinfo,lastEntry);
890                         } else {
891                             trlist->Add(&trinfo);
892                         }
893                         // check the track list
894                         Int_t nptracks=trlist->GetEntriesFast();
895                         if (nptracks > 2) {
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));
903                             }
904                         } // end if nptracks
905                       } //  end if pdigit
906                     } //end loop over clusters
907                 }// track type condition
908             } // hit loop
909         } // track loop
910         
911         // open the file with background
912         
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");     
917 //
918 //   Loop over tracks
919 //
920             for (Int_t trak=0; trak<ntracks; trak++) {
921                 if (fHits2)       fHits2->Clear();
922                 if (fClusters2)   fClusters2->Clear();
923                 TrH1->GetEvent(trak);
924 //
925 //   Loop over hits
926                 AliRICHHit* mHit;
927                 for(int j=0;j<fHits2->GetEntriesFast();++j) 
928                 {
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();
935 //
936 // Loop over pad hits
937                     for (AliRICHPadHit* mPad=
938                              (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
939                          mPad;
940                          mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2))
941                     {
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);
947 //
948 //
949                         Float_t thex, they;
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
957                         trinfo(1)=-1;
958                         digits[0]=ipx;
959                         digits[1]=ipy;
960                         digits[2]=iqpad;
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);
968                             
969                             pHitMap[nch]->SetHit(ipx, ipy, counter);
970                             counter++;
971                             printf("bgr new elem in list - counter %d\n",counter);
972                             
973                             pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
974                             // list of tracks
975                             TObjArray *trlist=(TObjArray*)pdigit->TrackList();
976                             trlist->Add(&trinfo);
977                         } else {
978                             pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
979                             // update charge
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));
987                             if (lastTrack==-1) {
988                                 continue;
989                             } else {
990                                 trlist->Add(&trinfo);
991                             }
992                             // check the track list
993                             Int_t nptracks=trlist->GetEntriesFast();
994                             if (nptracks > 0) {
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));
1000                                 }
1001                             } // end if nptracks
1002                         } //  end if pdigit
1003                     } //end loop over clusters
1004                 } // hit loop
1005             } // track loop
1006             TTree *fAli=gAlice->TreeK();
1007             if (fAli) pFile =fAli->GetCurrentFile();
1008             pFile->cd();
1009         } // if Add     
1010         
1011         Int_t tracks[10];
1012         Int_t charges[10];
1013         //cout<<"Start filling digits \n "<<endl;
1014         Int_t nentries=list->GetEntriesFast();
1015         //printf(" \n \n nentries %d \n",nentries);
1016
1017         // start filling the digits
1018         
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);
1030             q+=(Int_t)noise;
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;
1036             digits[2]=q;
1037             
1038             TObjArray* trlist=(TObjArray*)address->TrackList();
1039             Int_t nptracks=trlist->GetEntriesFast();
1040             
1041             // this was changed to accomodate the real number of tracks
1042             if (nptracks > 10) {
1043                 cout<<"Attention - tracks > 10 "<<nptracks<<endl;
1044                 nptracks=10;
1045             }
1046             if (nptracks > 2) {
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);
1050             }
1051             for (Int_t tr=0;tr<nptracks;tr++) {
1052                 TVector *ppP=(TVector*)trlist->At(tr);
1053                 TVector &pp  =*ppP;
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++) {
1059                     tracks[t]=0;
1060                     charges[t]=0;
1061                 }
1062             }
1063             //write file
1064             if (ich==2)
1065               fprintf(points,"%4d,      %4d,      %4d\n",digits[0],digits[1],digits[2]);
1066             
1067             // fill digits
1068             pRICH->AddDigits(ich,tracks,charges,digits);
1069         }
1070         gAlice->TreeD()->Fill();
1071
1072         list->Delete();
1073         for(Int_t ii=0;ii<kNCH;++ii) {
1074             if (pHitMap[ii]) {
1075                 hm=pHitMap[ii];
1076                 delete hm;
1077                 pHitMap[ii]=0;
1078             }
1079         }
1080         
1081         //TTree *TD=gAlice->TreeD();
1082         //Stat_t ndig=TD->GetEntries();
1083         //cout<<"number of digits  "<<ndig<<endl;
1084         TClonesArray *fDch;
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);
1089         }
1090         pRICH->ResetDigits();
1091     } //end loop over cathodes
1092     char hname[30];
1093     sprintf(hname,"TreeD%d",nev);
1094     gAlice->TreeD()->Write(hname);
1095
1096 // reset tree
1097 //    gAlice->TreeD()->Reset();
1098     delete list;
1099     pAddress->Clear();
1100 // gObjectTable->Print();
1101 }
1102
1103 AliRICH& AliRICH::operator=(const AliRICH& rhs)
1104 {
1105 // Assignment operator
1106     return *this;
1107     
1108 }
1109
1110 Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
1111 {
1112 //
1113 //  Calls the charge disintegration method of the current chamber and adds
1114 //  the simulated cluster to the root treee 
1115 //
1116     Int_t clhits[kNCH];
1117     Float_t newclust[6][500];
1118     Int_t nnew;
1119     
1120 //
1121 //  Integrated pulse height on chamber
1122     
1123     clhits[0]=fNhits+1;
1124     
1125     ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
1126     Int_t ic=0;
1127     
1128 //
1129 //  Add new clusters
1130     for (Int_t i=0; i<nnew; i++) {
1131         if (Int_t(newclust[3][i]) > 0) {
1132             ic++;
1133 // Cathode plane
1134             clhits[1] = Int_t(newclust[5][i]);
1135 //  Cluster Charge
1136             clhits[2] = Int_t(newclust[0][i]);
1137 //  Pad: ix
1138             clhits[3] = Int_t(newclust[1][i]);
1139 //  Pad: iy 
1140             clhits[4] = Int_t(newclust[2][i]);
1141 //  Pad: charge
1142             clhits[5] = Int_t(newclust[3][i]);
1143 //  Pad: chamber sector
1144             clhits[6] = Int_t(newclust[4][i]);
1145             
1146             AddPadHit(clhits);
1147         }
1148     }
1149 return nnew;
1150 }
1151