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