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=true;
 832     static TTree *TH1;
 833     static TFile *File;
 834     char *Add = strstr(option,"Add");
 835
 836     AliRICHchamber*  iChamber;
 837     AliRICHsegmentation*  segmentation;
 838
 839
 840     Int_t trk[50];
 841     Int_t chtrk[50];
 842     TObjArray *list=new TObjArray;
 843     Int_t digits[3];
 844
 845     AliRICH *RICH  = (AliRICH *) gAlice->GetDetector("RICH");
 846     AliRICHHitMap* HitMap[10];
 847     for (Int_t i=0; i<10; i++) {HitMap[i]=0;}
 848     if (Add ) {
 849         if(first) {
 850             fFileName=filename;
 851             cout<<"filename"<<fFileName<<endl;
 852             File=new TFile(fFileName);
 853             cout<<"I have opened "<<fFileName<<" file "<<endl;
 854             fHits2     = new TClonesArray("AliRICHhit",1000  );
 855             fClusters2 = new TClonesArray("AliRICHcluster",10000);
 856             first=false;
 857         }
 858         File->cd();
 859         File->ls();
 860         // Get Hits Tree header from file
 861         if(fHits2) fHits2->Clear();
 862         if(fClusters2) fClusters2->Clear();
 863         if(TH1) delete TH1;
 864         TH1=0;
 865         //
 866         char treeName[20];
 867         sprintf(treeName,"TreeH%d",nev);
 868         TH1 = (TTree*)gDirectory->Get(treeName);
 869         if (!TH1) {
 870             printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
 871         }
 872         // Set branch addresses
 873         TBranch *branch;
 874         char branchname[20];
 875         sprintf(branchname,"%s",GetName());
 876         if (TH1 && fHits2) {
 877             branch = TH1->GetBranch(branchname);
 878             if (branch) branch->SetAddress(&fHits2);
 879         }
 880         if (TH1 && fClusters2) {
 881             branch = TH1->GetBranch("RICHCluster");
 882             if (branch) branch->SetAddress(&fClusters2);
 883         }
 884     }
 885     //
 886     // loop over cathodes
 887     //
 888     AliRICHHitMap* hm;
 889
 890     for (int icat=0; icat<1; icat++) {
 891         for (Int_t i=0; i<7; i++) {
 892             if (HitMap[i]) {
 893                 hm=HitMap[i];
 894                 delete hm;
 895                 HitMap[i]=0;
 896             }
 897         }
 898         Int_t counter=0;
 899         for (Int_t i =0; i<7; i++) {
 900             iChamber=(AliRICHchamber*) (*fChambers)[i];
 901             if (iChamber->Nsec()==1 && icat==1) {
 902                 continue;
 903             } else {
 904                 segmentation=iChamber->GetSegmentationModel(icat+1);
 905             }
 906             HitMap[i] = new AliRICHHitMapA1(segmentation, list);
 907         }
 908         printf("Start loop over tracks \n");
 909 //
 910 //   Loop over tracks
 911 //
 912
 913         TTree *TH = gAlice->TreeH();
 914         Int_t ntracks =(Int_t) TH->GetEntries();
 915         for (Int_t track=0; track<ntracks; track++) {
 916             gAlice->ResetHits();
 917             TH->GetEvent(track);
 918 //
 919 //   Loop over hits
 920             for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1);
 921                 mHit;
 922                 mHit=(AliRICHhit*)RICH->NextHit())
 923             {
 924                 Int_t   nch   = mHit->fChamber-1;  // chamber number
 925                 if (nch >7) continue;
 926                 iChamber = &(RICH->Chamber(nch));
 927
 928 //
 929 // Loop over pad hits
 930                 for (AliRICHcluster* mPad=
 931                          (AliRICHcluster*)RICH->FirstPad(mHit,fClusters);
 932                      mPad;
 933                      mPad=(AliRICHcluster*)RICH->NextPad(fClusters))
 934                 {
 935                     Int_t cathode  = mPad->fCathode;    // cathode number
 936                     Int_t ipx      = mPad->fPadX;       // pad number on X
 937                     Int_t ipy      = mPad->fPadY;       // pad number on Y
 938                     Int_t iqpad    = mPad->fQpad;       // charge per pad
 939 //
 940 //
 941
 942                     if (cathode != (icat+1)) continue;
 943                     // fill the info array
 944                     Float_t thex, they;
 945                     segmentation=iChamber->GetSegmentationModel(cathode);
 946                     segmentation->GetPadCxy(ipx,ipy,thex,they);
 947                     TVector *trinfo_p= new TVector(2);
 948                     TVector &trinfo = *trinfo_p;
 949                     trinfo(0)=(Float_t)track;
 950                     trinfo(1)=(Float_t)iqpad;
 951
 952                     digits[0]=ipx;
 953                     digits[1]=ipy;
 954                     digits[2]=iqpad;
 955
 956                     AliRICHlist* pdigit;
 957                     // build the list of fired pads and update the info
 958                     if (!HitMap[nch]->TestHit(ipx, ipy)) {
 959                         list->AddAtAndExpand(
 960                             new AliRICHlist(nch,digits),counter);
 961                         HitMap[nch]->SetHit(ipx, ipy, counter);
 962                         counter++;
 963                         pdigit=(AliRICHlist*)list->At(list->GetLast());
 964                         // list of tracks
 965                         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
 966                         trlist->Add(&trinfo);
 967                     } else {
 968                         pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
 969                         // update charge
 970                         (*pdigit).fSignal+=iqpad;
 971                         // update list of tracks
 972                         TObjArray* trlist=(TObjArray*)pdigit->TrackList();
 973                         Int_t last_entry=trlist->GetLast();
 974                         TVector *ptrk_p=(TVector*)trlist->At(last_entry);
 975                         TVector &ptrk=*ptrk_p;
 976                         Int_t last_track=Int_t(ptrk(0));
 977                         Int_t last_charge=Int_t(ptrk(1));
 978                         if (last_track==track) {
 979                             last_charge+=iqpad;
 980                             trlist->RemoveAt(last_entry);
 981                             trinfo(0)=last_track;
 982                             trinfo(1)=last_charge;
 983                             trlist->AddAt(&trinfo,last_entry);
 984                         } else {
 985                             trlist->Add(&trinfo);
 986                         }
 987                         // check the track list
 988                         Int_t nptracks=trlist->GetEntriesFast();
 989                         if (nptracks > 2) {
 990                             printf("Attention - nptracks > 2  %d \n",nptracks);
 991                             printf("cat,nch,ix,iy %d %d %d %d  \n",icat+1,nch,ipx,ipy);
 992                             for (Int_t tr=0;tr<nptracks;tr++) {
 993                                 TVector *pptrk_p=(TVector*)trlist->At(tr);
 994                                 TVector &pptrk=*pptrk_p;
 995                                 trk[tr]=Int_t(pptrk(0));
 996                                 chtrk[tr]=Int_t(pptrk(1));
 997                             }
 998                         } // end if nptracks
 999                     } //  end if pdigit
1000                 } //end loop over clusters
1001             } // hit loop
1002         } // track loop
1003
1004         Int_t nentr1=list->GetEntriesFast();
1005         printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
1006
1007         // open the file with background
1008
1009         if (Add ) {
1010             ntracks =(Int_t)TH1->GetEntries();
1011             printf("background - icat,ntracks1  %d %d\n",icat,ntracks);
1012             printf("background - Start loop over tracks \n");
1013 //
1014 //   Loop over tracks
1015 //
1016             for (Int_t track=0; track<ntracks; track++) {
1017
1018                 if (fHits2)       fHits2->Clear();
1019                 if (fClusters2)   fClusters2->Clear();
1020
1021                 TH1->GetEvent(track);
1022 //
1023 //   Loop over hits
1024                 AliRICHhit* mHit;
1025                 for(int i=0;i<fHits2->GetEntriesFast();++i)
1026                 {
1027                     mHit=(AliRICHhit*) (*fHits2)[i];
1028                     Int_t   nch   = mHit->fChamber-1;  // chamber number
1029                     if (nch >9) continue;
1030                     iChamber = &(RICH->Chamber(nch));
1031                     Int_t rmin = (Int_t)iChamber->RInner();
1032                     Int_t rmax = (Int_t)iChamber->ROuter();
1033 //
1034 // Loop over pad hits
1035                     for (AliRICHcluster* mPad=
1036                              (AliRICHcluster*)RICH->FirstPad(mHit,fClusters2);
1037                          mPad;
1038                          mPad=(AliRICHcluster*)RICH->NextPad(fClusters2))
1039                     {
1040                         Int_t cathode  = mPad->fCathode;    // cathode number
1041                         Int_t ipx      = mPad->fPadX;       // pad number on X
1042                         Int_t ipy      = mPad->fPadY;       // pad number on Y
1043                         Int_t iqpad    = mPad->fQpad;       // charge per pad
1044                         if (track==3 && nch==0 && icat==0) printf("bgr - track,iqpad,ipx,ipy %d %d %d %d\n",track,iqpad,ipx,ipy);
1045 //
1046 //
1047                         if (cathode != (icat+1)) continue;
1048                         // fill the info array
1049                         Float_t thex, they;
1050                         segmentation=iChamber->GetSegmentationModel(cathode);
1051                         segmentation->GetPadCxy(ipx,ipy,thex,they);
1052                         Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1053                         if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1054
1055                         TVector *trinfo_p;
1056                         trinfo_p = new TVector(2);
1057                         TVector &trinfo = *trinfo_p;
1058                         trinfo(0)=-1;  // tag background
1059                         trinfo(1)=-1;
1060
1061                         digits[0]=ipx;
1062                         digits[1]=ipy;
1063                         digits[2]=iqpad;
1064
1065
1066                         if (track <4 && icat==0 && nch==0)
1067                             printf("bgr - HitMap[nch]->TestHit(ipx, ipy),track %d %d\n",
1068                                    HitMap[nch]->TestHit(ipx, ipy),track);
1069                         AliRICHlist* pdigit;
1070                         // build the list of fired pads and update the info
1071                         if (!HitMap[nch]->TestHit(ipx, ipy)) {
1072                             list->AddAtAndExpand(new AliRICHlist(nch,digits),counter);
1073
1074                             HitMap[nch]->SetHit(ipx, ipy, counter);
1075                             counter++;
1076                             printf("bgr new elem in list - counter %d\n",counter);
1077
1078                             pdigit=(AliRICHlist*)list->At(list->GetLast());
1079                             // list of tracks
1080                             TObjArray *trlist=(TObjArray*)pdigit->TrackList();
1081                             trlist->Add(&trinfo);
1082                         } else {
1083                             pdigit=(AliRICHlist*) HitMap[nch]->GetHit(ipx, ipy);
1084                             // update charge
1085                             (*pdigit).fSignal+=iqpad;
1086                             // update list of tracks
1087                             TObjArray* trlist=(TObjArray*)pdigit->TrackList();
1088                             Int_t last_entry=trlist->GetLast();
1089                             TVector *ptrk_p=(TVector*)trlist->At(last_entry);
1090                             TVector &ptrk=*ptrk_p;
1091                             Int_t last_track=Int_t(ptrk(0));
1092                             if (last_track==-1) {
1093                                 continue;
1094                             } else {
1095                                 trlist->Add(&trinfo);
1096                             }
1097                             // check the track list
1098                             Int_t nptracks=trlist->GetEntriesFast();
1099                             if (nptracks > 0) {
1100                                 for (Int_t tr=0;tr<nptracks;tr++) {
1101                                     TVector *pptrk_p=(TVector*)trlist->At(tr);
1102                                     TVector &pptrk=*pptrk_p;
1103                                     trk[tr]=Int_t(pptrk(0));
1104                                     chtrk[tr]=Int_t(pptrk(1));
1105                                 }
1106                             } // end if nptracks
1107                         } //  end if pdigit
1108                     } //end loop over clusters
1109                 } // hit loop
1110             } // track loop
1111             Int_t nentr2=list->GetEntriesFast();
1112             printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1113             TTree *fAli=gAlice->TreeK();
1114             if (fAli) File =fAli->GetCurrentFile();
1115             File->cd();
1116         } // if Add
1117
1118         Int_t tracks[10];
1119         Int_t charges[10];
1120         cout<<"start filling digits \n "<<endl;
1121         Int_t nentries=list->GetEntriesFast();
1122         printf(" \n \n nentries %d \n",nentries);
1123
1124         // start filling the digits
1125
1126         for (Int_t nent=0;nent<nentries;nent++) {
1127             AliRICHlist *address=(AliRICHlist*)list->At(nent);
1128             if (address==0) continue;
1129             Int_t ich=address->fChamber;
1130             Int_t q=address->fSignal;
1131             iChamber=(AliRICHchamber*) (*fChambers)[ich];
1132             // add white noise and do zero-suppression and signal truncation
1133             Float_t MeanNoise = gRandom->Gaus(1, 0.2);
1134             Float_t ZeroSupp=5*MeanNoise;
1135             Float_t Noise     = gRandom->Gaus(0, MeanNoise);
1136             q+=(Int_t)Noise;
1137             if ( q <= ZeroSupp) continue;
1138             digits[0]=address->fPadX;
1139             digits[1]=address->fPadY;
1140             digits[2]=q;
1141
1142             TObjArray* trlist=(TObjArray*)address->TrackList();
1143             Int_t nptracks=trlist->GetEntriesFast();
1144
1145             // this was changed to accomodate the real number of tracks
1146             if (nptracks > 10) {
1147                 cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1148                 nptracks=10;
1149             }
1150             if (nptracks > 2) {
1151                 printf("Attention - nptracks > 2  %d \n",nptracks);
1152                 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1153             }
1154             for (Int_t tr=0;tr<nptracks;tr++) {
1155                 TVector *pp_p=(TVector*)trlist->At(tr);
1156                 TVector &pp  =*pp_p;
1157                 tracks[tr]=Int_t(pp(0));
1158                 charges[tr]=Int_t(pp(1));
1159             }      //end loop over list of tracks for one pad
1160             if (nptracks < 10 ) {
1161                 for (Int_t i=nptracks; i<10; i++) {
1162                     tracks[i]=0;
1163                     charges[i]=0;
1164                 }
1165             }
1166             // fill digits
1167             RICH->AddDigits(ich,tracks,charges,digits);
1168
1169             delete address;
1170         }
1171         cout<<"I'm out of the loops for digitisation"<<endl;
1172         gAlice->TreeD()->Fill();
1173         TTree *TD=gAlice->TreeD();
1174         Stat_t ndig=TD->GetEntries();
1175         cout<<"number of digits  "<<ndig<<endl;
1176         TClonesArray *fDch;
1177         for (int i=0;i<7;i++) {
1178             fDch= RICH->DigitsAddress(i);
1179             int ndig=fDch->GetEntriesFast();
1180             printf (" i, ndig %d %d \n",i,ndig);
1181         }
1182         RICH->ResetDigits();
1183
1184         list->Clear();
1185
1186     } //end loop over cathodes
1187     char hname[30];
1188     sprintf(hname,"TreeD%d",nev);
1189     gAlice->TreeD()->Write(hname);
1190 }
1191
1192
1193
1194
1195
1196
1197
1198
1199