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