RICH/AliRICHDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2000, 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 //Piotr.Skowronski@cern.ch :
  18 //Corrections applied in order to compile (only) with new I/O and folder structure
  19 //To be implemented correctly by responsible
  20
  21 #include <Riostream.h>
  22
  23 #include <TTree.h>
  24 #include <TObjArray.h>
  25 #include <TFile.h>
  26 #include <TDirectory.h>
  27 #include <TParticle.h>
  28
  29 #include "AliRunLoader.h"
  30 #include "AliLoader.h"
  31
  32 #include "AliRICHDigitizer.h"
  33 #include "AliRICHChamber.h"
  34 #include "AliHitMap.h"
  35 #include "AliRICHHitMapA1.h"
  36 #include "AliRICH.h"
  37 #include "AliRICHSDigit.h"
  38 #include "AliRICHDigit.h"
  39 #include "AliRICHTransientDigit.h"
  40 #include "AliRun.h"
  41 #include "AliPDG.h"
  42 #include "AliRunDigitizer.h"
  43
  44 ClassImp(AliRICHDigitizer)
  45
  46 //__________________________________________________________________________________________________
  47 AliRICHDigitizer::AliRICHDigitizer()
  48 {//default ctor - don't use it
  49   fHits = 0;
  50   fSDigits = 0;
  51   fHitMap = 0;
  52   fTDList = 0;
  53 }//default ctor
  54 //__________________________________________________________________________________________________
  55 AliRICHDigitizer::AliRICHDigitizer(AliRunDigitizer *pManager)
  56                  :AliDigitizer(pManager)
  57 {//main ctor which should be used
  58   if(fManager->GetDebug())Info("main ctor","Start.");
  59   fHits = 0;
  60   fSDigits = 0;
  61   fHitMap = 0;
  62   fTDList = 0;
  63   fDebug  = pManager->GetDebug();
  64 }//main ctor
  65 //__________________________________________________________________________________________________
  66 AliRICHDigitizer::~AliRICHDigitizer()
  67 {//dtor
  68   if(fManager->GetDebug())Info("dtor","Start.");
  69
  70   if(fHits)    {fHits->Delete();   delete fHits;}
  71   if(fSDigits) {fSDigits->Delete();delete fSDigits;}
  72   if(fTDList)  {fTDList->Delete(); delete fTDList;}
  73   for(Int_t i=0; i<kNCH; i++) delete fHitMap[i];    delete [] fHitMap;
  74 }//dtor
  75 //__________________________________________________________________________________________________
  76 Bool_t AliRICHDigitizer::Exists(const AliRICHSDigit *p)
  77 {
  78   return (fHitMap[fNch]->TestHit(p->PadX(),p->PadY()));
  79 }//Exists
  80 //__________________________________________________________________________________________________
  81 void AliRICHDigitizer::Update(AliRICHSDigit *padhit)
  82 {
  83   AliRICHTransientDigit *pdigit =
  84     static_cast<AliRICHTransientDigit*>(
  85       fHitMap[fNch]->GetHit(padhit->PadX(),padhit->PadY()));
  86
  87   // update charge
  88   //
  89   Int_t iqpad    = Int_t(padhit->QPad()); // charge per pad
  90   pdigit->AddSignal(iqpad);
  91   pdigit->AddPhysicsSignal(iqpad);
  92
  93   // update list of tracks
  94   //
  95   Int_t track, charge;
  96   track = fTrack+fMask;
  97   if (fSignal) {
  98     charge = iqpad;
  99   } else {
 100     charge = kBgTag;
 101   }
 102   pdigit->UpdateTrackList(track,charge);
 103 }//Update()
 104 //__________________________________________________________________________________________________
 105 void AliRICHDigitizer::CreateNew(AliRICHSDigit *padhit)
 106 {
 107   fTDList->AddAtAndExpand(
 108     new AliRICHTransientDigit(fNch,fDigits),fCounter);
 109   fHitMap[fNch]->SetHit(padhit->PadX(),padhit->PadY(),fCounter);
 110
 111   AliRICHTransientDigit* pdigit =
 112     static_cast<AliRICHTransientDigit*>(fTDList->Last());
 113   // list of tracks
 114   Int_t track, charge;
 115   if (fSignal) {
 116     track = fTrack;
 117     charge = padhit->QPad();
 118   } else {
 119     track = kBgTag;
 120     charge = kBgTag;
 121   }
 122   pdigit->AddToTrackList(track,charge);
 123   fCounter++;
 124 }//CreateNew()
 125 //__________________________________________________________________________________________________
 126 Bool_t AliRICHDigitizer::Init()
 127 {// Initialisation
 128   if(GetDebug())Info("Init","Start.");
 129   fHits     = new TClonesArray("AliRICHhit",1000);
 130   fSDigits  = new TClonesArray("AliRICHSDigit",1000);
 131   return kTRUE;
 132 }//Init()
 133 //__________________________________________________________________________________________________
 134 void AliRICHDigitizer::Exec(Option_t* option)
 135 {
 136   if(GetDebug())Info("Exec","Start with option=%s",option);
 137
 138   AliRICHChamber*       iChamber;
 139   AliSegmentation*  segmentation;
 140
 141   AliRunLoader *pInAL, *pOutAL;//in and out Run loaders
 142   AliLoader    *pInRL, *pOutRL;//in and out RICH loaders
 143
 144   pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 145   pOutRL = pOutAL->GetLoader("RICHLoader");
 146
 147   fTDList = new TObjArray;
 148
 149
 150   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
 151
 152   if(!pOutRL->TreeD()) pOutRL->MakeTree("D");  pRICH->MakeBranch("D");
 153
 154
 155   fHitMap= new AliHitMap* [kNCH];
 156
 157   for (Int_t i =0; i<kNCH; i++) {
 158     iChamber= &(pRICH->Chamber(i));
 159     segmentation=iChamber->GetSegmentationModel();
 160     fHitMap[i] = new AliRICHHitMapA1(segmentation, fTDList);
 161   }
 162
 163 // Loop over files to digitize
 164   fSignal = kTRUE;
 165   fCounter = 0;
 166   for (Int_t inputFile=0;inputFile<fManager->GetNinputs();inputFile++){//files loop
 167     pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
 168     pInRL = pInAL->GetLoader("RICHLoader");
 169
 170 // Connect RICH branches
 171
 172     if (inputFile > 0 ) fSignal = kFALSE;
 173     TBranch *branchHits = 0;
 174     TBranch *branchSDigits = 0;
 175
 176     pInRL->LoadHits("READ");
 177
 178     TTree *treeH = pInRL->TreeH();
 179     if(pInRL->TreeH()) {
 180       branchHits = treeH->GetBranch("RICH");
 181       if(branchHits){
 182         fHits->Clear();
 183         branchHits->SetAddress(&fHits);}
 184       else
 185         Error("Exec","branch RICH was not found");
 186     }
 187
 188     if(treeH) {
 189       branchSDigits = treeH->GetBranch("RICHSpecials");
 190       if(branchSDigits)
 191         branchSDigits->SetAddress(&fSDigits);
 192       else
 193         Error("exec","branch RICHSpecials was not found");
 194     }
 195
 196     for (fTrack=0; fTrack<treeH->GetEntries(); fTrack++) {//prims loop
 197       fHits->Clear();                   fSDigits->Clear();
 198       fHits->Print();
 199
 200       branchHits->GetEntry(fTrack);     branchSDigits->GetEntry(fTrack);
 201       for(Int_t i=0;i<fHits->GetEntriesFast();++i){//hits loop
 202         AliRICHhit* pHit = static_cast<AliRICHhit*>(fHits->At(i));
 203         fNch = pHit->Chamber()-1;  // chamber number
 204         if (fNch >= kNCH) {
 205           cerr<<"AliRICHDigitizer: chamber nr. fNch out of range: "<<fNch<<endl;
 206           cerr<<"               track: "<<fTrack<<endl;
 207           continue;
 208         }
 209         iChamber = &(pRICH->Chamber(fNch));
 210
 211         for (AliRICHSDigit* mPad=FirstPad(pHit,fSDigits);mPad;mPad=NextPad(fSDigits)){//specials loop
 212           Int_t iqpad    = mPad->QPad();       // charge per pad
 213           fDigits[0]=mPad->PadX();
 214           fDigits[1]=mPad->PadY();
 215           fDigits[2]=iqpad;
 216           fDigits[3]=iqpad;
 217           fDigits[4]=mPad->HitNumber();
 218           if (Exists(mPad))       // build the list of fired pads and update the info
 219             Update(mPad);
 220           else
 221             CreateNew(mPad);
 222         }//specials loop
 223       }//hits loop
 224     }//prims loop
 225   }//files loop
 226   if(GetDebug()) Info("Exec","END OF FILE LOOP");
 227
 228   Int_t tracks[kMAXTRACKSPERRICHDIGIT];
 229   Int_t charges[kMAXTRACKSPERRICHDIGIT];
 230   Int_t nentries=fTDList->GetEntriesFast();
 231
 232   for (Int_t nent=0;nent<nentries;nent++){//transient digits loop
 233     AliRICHTransientDigit *transDigit=(AliRICHTransientDigit*)fTDList->At(nent);
 234     if (transDigit==0) continue;
 235     Int_t ich=transDigit->GetChamber();
 236     Int_t q=transDigit->Signal();
 237     iChamber=&(pRICH->Chamber(ich));
 238     AliRICHResponse * response=iChamber->GetResponseModel();
 239     Int_t adcmax= (Int_t) response->MaxAdc();
 240
 241
 242     // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
 243     //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY()));
 244
 245 // tmp change
 246 //    Int_t pedestal = iChamber->fTresh->GetHitIndex(transDigit->PadX(),transDigit->PadY());
 247     Int_t pedestal = 0;
 248
 249     Float_t treshold = (pedestal + 4*2.2);
 250
 251     Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
 252     Float_t noise     = gRandom->Gaus(0, meanNoise);
 253     q+=(Int_t)(noise + pedestal);
 254     //q+=(Int_t)(noise);
 255     //          magic number to be parametrised !!!
 256     if ( q <= treshold)
 257     {
 258       q = q - pedestal;
 259       continue;
 260     }
 261     q = q - pedestal;
 262     if ( q >= adcmax) q=adcmax;
 263     fDigits[0]=transDigit->PadX();
 264     fDigits[1]=transDigit->PadY();
 265     fDigits[2]=q;
 266     fDigits[3]=transDigit->Physics();
 267     fDigits[4]=transDigit->Hit();
 268
 269     Int_t nptracks = transDigit->GetNTracks();
 270
 271
 272     if(nptracks>kMAXTRACKSPERRICHDIGIT){// this was changed to accomodate the real number of tracks
 273       printf("Attention - tracks > 10 %d\n",nptracks);
 274       nptracks=kMAXTRACKSPERRICHDIGIT;
 275     }
 276     if (nptracks > 2) {
 277       printf("Attention - tracks > 2  %d \n",nptracks);
 278     }
 279     for (Int_t tr=0;tr<nptracks;tr++) {
 280       tracks[tr]=transDigit->GetTrack(tr);
 281       charges[tr]=transDigit->GetCharge(tr);
 282     }      //end loop over list of tracks for one pad
 283     if (nptracks < kMAXTRACKSPERRICHDIGIT ) {
 284       for (Int_t t=nptracks; t<kMAXTRACKSPERRICHDIGIT; t++) {
 285       tracks[t]=0;
 286       charges[t]=0;
 287       }
 288     }
 289     pRICH->AddDigitOld(ich+1,tracks,charges,fDigits);//OLD
 290   }//transient digits loop
 291   pOutRL->TreeD()->Fill();
 292
 293   //pRICH->ResetDigits();
 294   fTDList->Delete();    // or fTDList->Clear(); ???
 295   for(Int_t ii=0;ii<kNCH;++ii) {
 296     if (fHitMap[ii]) {
 297       delete fHitMap[ii];
 298       fHitMap[ii]=0;
 299     }
 300   }
 301
 302   TClonesArray *richDigits;
 303   for (Int_t k=0;k<kNCH;k++) {
 304     richDigits = pRICH->DigitsAddress(k);
 305     Int_t ndigit=richDigits->GetEntriesFast();
 306     printf ("Chamber %d digits %d \n",k+1,ndigit);
 307   }
 308   pRICH->ResetDigits(); /// ??? should it be here???
 309
 310   pOutRL->WriteDigits("OVERWRITE");
 311
 312   delete [] fHitMap;
 313   delete fTDList;
 314   pOutRL->UnloadHits();
 315   pOutRL->UnloadDigits();
 316 //  if (fHits)    fHits->Delete();
 317 //  if (fSDigits) fSDigits->Delete();
 318   if(GetDebug())Info("Exec","Stop.");
 319 }//Exec()
 320 //__________________________________________________________________________________________________
 321
 322 static Int_t sMaxIterPad=0;    // Static variables for the pad-hit iterator routines
 323 static Int_t sCurIterPad=0;
 324
 325 //__________________________________________________________________________________________________
 326 AliRICHSDigit* AliRICHDigitizer::FirstPad(AliRICHhit*  hit,TClonesArray *clusters )
 327 {// Initialise the pad iterator Return the address of the first sdigit for hit
 328     TClonesArray *theClusters = clusters;
 329     Int_t nclust = theClusters->GetEntriesFast();
 330     if (nclust && hit->PHlast() > 0) {
 331         sMaxIterPad=Int_t(hit->PHlast());
 332         sCurIterPad=Int_t(hit->PHfirst());
 333         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
 334     } else {
 335         return 0;
 336     }
 337
 338 }//FirstPad
 339 //__________________________________________________________________________________________________
 340 AliRICHSDigit* AliRICHDigitizer::NextPad(TClonesArray *clusters)
 341 {// Iterates over pads
 342
 343     sCurIterPad++;
 344     if (sCurIterPad <= sMaxIterPad) {
 345         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
 346     } else {
 347         return 0;
 348     }
 349 }//NextPad
 350 //__________________________________________________________________________________________________