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