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