START/AliSTARTDigitizer.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17
  18 #include <TTree.h>
  19 #include <TFile.h>
  20 #include <TDirectory.h>
  21 #include <TRandom.h>
  22 #include <TArrayI.h>
  23 #include <TError.h>
  24
  25
  26 #include "AliSTARTDigitizer.h"
  27 #include "AliSTART.h"
  28 #include "AliSTARThit.h"
  29 #include "AliSTARThitPhoton.h"
  30 #include "AliSTARTdigit.h"
  31 #include "AliRunDigitizer.h"
  32 #include <AliDetector.h>
  33 #include "AliRun.h"
  34 #include <AliLoader.h>
  35 #include <AliRunLoader.h>
  36 #include <stdlib.h>
  37 #include <Riostream.h>
  38 #include <Riostream.h>
  39
  40 ClassImp(AliSTARTDigitizer)
  41
  42 //___________________________________________
  43   AliSTARTDigitizer::AliSTARTDigitizer()  :AliDigitizer()
  44 {
  45 // Default ctor - don't use it
  46   ;
  47 }
  48
  49 //___________________________________________
  50 AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
  51     :AliDigitizer(manager)
  52 {
  53   //    cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
  54 // ctor which should be used
  55 //  fDebug =0;
  56   if (GetDebug()>2)
  57     cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
  58         <<"(AliRunDigitizer* manager) was processed"<<endl;
  59
  60   ftimeRightTDC = new TArrayI(12);
  61   ftimeLeftTDC = new TArrayI(12);
  62   fRightADC = new TArrayI(12);
  63   fLeftADC = new TArrayI(12);
  64 }
  65
  66 //------------------------------------------------------------------------
  67 AliSTARTDigitizer::~AliSTARTDigitizer()
  68 {
  69 // Destructor
  70   if(GetDebug()) Info("dtor","START");
  71   delete ftimeRightTDC;
  72   delete ftimeLeftTDC;
  73   delete fRightADC;
  74   delete fLeftADC;
  75 }
  76
  77  //------------------------------------------------------------------------
  78 Bool_t AliSTARTDigitizer::Init()
  79 {
  80 // Initialization
  81 // cout<<"AliSTARTDigitizer::Init"<<endl;
  82  return kTRUE;
  83 }
  84
  85
  86 //---------------------------------------------------------------------
  87
  88 void AliSTARTDigitizer::Exec(Option_t* /*option*/)
  89 {
  90
  91   /*
  92     Produde digits from hits
  93         digits is TObject and includes
  94         We are writing array if left & right  TDC
  95         left & right  ADC (will need for slow simulation)
  96         TOF first particle left & right
  97         mean time and time difference (vertex position)
  98
  99   */
 100
 101   AliRunLoader *inRL, *outRL;//in and out Run Loaders
 102   AliLoader *pInStartLoader, *pOutStartLoader;// in and out STARTLoaders
 103
 104   outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 105   pOutStartLoader = outRL->GetLoader("STARTLoader");
 106
 107 #ifdef DEBUG
 108   cout<<"AliSTARTDigitizer::>Hits2Digits start...\n";
 109 #endif
 110   //
 111   // From hits to digits
 112   //
 113   Int_t hit, nhits;
 114   Float_t meanTime;
 115   Int_t countEr[13],countEl[13];
 116   Int_t volume,pmt,tr,tl,sumRight;
 117   Int_t  bestRightADC,bestLeftADC;
 118   Float_t besttimeleftGaus, besttimerightGaus;
 119   Float_t timeright[13]={13*0};
 120   Float_t timeleft[13]={13*0};
 121   Float_t channelWidth=2.5; //ps
 122   Int_t channelWidthADC=1; //ps
 123   //  Int_t thresholdAmpl=10;
 124
 125
 126   inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
 127
 128
 129   AliSTART *fSTART  = (AliSTART*) inRL->GetAliRun()->GetDetector("START");
 130   AliSTARThit  *startHit;
 131   TBranch *brHits=0;
 132   TBranch *brHitPhoton=0;
 133   fdigits= new AliSTARTdigit();
 134
 135   Int_t nFiles=fManager->GetNinputs();
 136   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
 137
 138     Float_t besttimeright=9999.;
 139     Float_t besttimeleft=9999.;
 140     Int_t iTimeDiff=0;
 141     Int_t iTimeAv=0;
 142     Float_t timeDiff,timeAv;
 143     sumRight=0;
 144     for (Int_t i0=0; i0<13; i0++)
 145       {
 146         timeright[i0]=0; timeleft[i0]=0;
 147         countEr[i0]=0;   countEl[i0]=0;
 148       }
 149
 150     inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
 151     pInStartLoader = inRL->GetLoader("STARTLoader");
 152     pInStartLoader->LoadHits("READ");//probably it is necessary to load them before
 153     pOutStartLoader->LoadDigits("UPDATE");//probably it is necessary to load them before
 154     TClonesArray *fHits = fSTART->Hits ();
 155
 156     TTree *th = pInStartLoader->TreeH();
 157     brHits = th->GetBranch("START");
 158     brHitPhoton = th->GetBranch("STARThitPhoton");
 159     if (brHits) {
 160       fSTART->SetHitsAddressBranch(brHits,brHitPhoton);
 161     }else{
 162       cerr<<"EXEC Branch START hit not found"<<endl;
 163       exit(111);
 164     }
 165     Int_t ntracks    = (Int_t) th->GetEntries();
 166 #ifdef DEBUG
 167    Info("Digitizer",ntracks);
 168 #endif
 169      if (ntracks<=0) return;
 170     // Start loop on tracks in the hits containers
 171     for (Int_t track=0; track<ntracks;track++) {
 172       brHits->GetEntry(track);
 173       nhits = fHits->GetEntriesFast();
 174       for (hit=0;hit<nhits;hit++) {
 175         startHit   = (AliSTARThit*) fHits->UncheckedAt(hit);
 176         if (!startHit) {
 177           ::Error("Exec","The unchecked hit doesn't exist");
 178           break;
 179         }
 180         pmt=startHit->Pmt();
 181         volume = startHit->Volume();
 182         if(volume==1){
 183           timeright[pmt] = startHit->Time();
 184           if(timeright[pmt]<besttimeright)
 185             {
 186               besttimeright=timeright[pmt];
 187           } //timeright
 188         }//time for right shoulder
 189         if(volume==2){
 190           timeleft[pmt] = startHit->Time();
 191           if(timeleft[pmt]<besttimeleft)
 192             {
 193               besttimeleft=timeleft[pmt];
 194
 195           } //timeleftbest
 196         }//time for left shoulder
 197       } //hit loop
 198     } //track loop
 199
 200     // z position
 201
 202     //folding with experimental time distribution
 203
 204     Float_t koef=69.7/350.;
 205     besttimeright=koef*besttimeright;
 206     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
 207     bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
 208     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
 209     bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
 210     timeDiff=besttimerightGaus-besttimeleftGaus;
 211 #ifdef DEBUG
 212     cout<<" timediff in ns "<<timeDiff<<" z= "<<timeDiff*30<<endl;
 213 #endif
 214     meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
 215     if ( TMath::Abs(timeDiff)<TMath::Abs(0.3) )
 216       {
 217         Float_t t1=1000.*besttimeleftGaus;
 218         Float_t t2=1000.*besttimerightGaus;
 219         t1=t1/channelWidth;   //time in ps to channelWidth
 220         t2=t2/channelWidth;   //time in ps to channelWidth
 221         timeAv=(t1+t2)/2.;// time  channel numbres
 222
 223         // Time to TDC signal
 224         // 256 channels for timediff, range 1ns
 225         iTimeAv=(Int_t)timeAv;
 226         timeDiff= 512+1000*timeDiff/channelWidth; // time  channel numbres
 227         iTimeDiff=(Int_t)timeDiff;
 228         //       fill digits
 229         fdigits->SetTimeBestLeft(bestLeftADC);
 230         fdigits->SetTimeBestRight(bestRightADC);
 231         fdigits->SetMeanTime(iTimeAv);
 232         fdigits->SetTimeDiff(iTimeDiff);
 233         for (Int_t i=0; i<12; i++)
 234           {
 235             //  fill TDC
 236             timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
 237             timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
 238             tr= Int_t (timeright[i+1]*1000/channelWidth);
 239             if(tr<200) tr=0;
 240             tl= Int_t (timeleft[i+1]*1000/channelWidth);
 241             if(tl<1000) tl=0;
 242
 243             ftimeRightTDC->AddAt(tr,i);
 244             ftimeLeftTDC->AddAt(tl,i);
 245             //fill ADC
 246             Int_t al=( Int_t ) countEl[i+1]/ channelWidthADC;
 247             Int_t ar=( Int_t ) countEr[i+1]/ channelWidthADC;
 248             fRightADC->AddAt(ar,i);
 249             fLeftADC ->AddAt(al,i);
 250             sumRight+=countEr[i+1];
 251           }
 252         fdigits->SetTimeRight(*ftimeRightTDC);
 253         fdigits->SetTimeLeft(*ftimeLeftTDC);
 254         fdigits->SetADCRight(*fRightADC);
 255         fdigits->SetADCLeft(*fLeftADC);
 256         fdigits->SetSumADCRight(sumRight);
 257       }
 258     else
 259       {timeAv=999999; timeDiff=99999;}
 260
 261 // trick to find out output dir:
 262
 263     TDirectory *wd = gDirectory;
 264     pOutStartLoader->GetDigitsDataLoader()->GetDirectory()->cd();
 265     fdigits->Write("START_D");
 266     wd->cd();
 267     pInStartLoader->UnloadHits();
 268     pOutStartLoader->UnloadDigits();
 269   } //event loop
 270 }
 271
 272
 273 //------------------------------------------------------------------------
 274 Bool_t AliSTARTDigitizer::RegisterPhotoE(/*AliSTARThitPhoton *hit*/)
 275 {
 276     Double_t    pP = 0.2;
 277     Double_t    p;
 278
 279     p = gRandom->Rndm();
 280     if (p > pP)
 281       return kFALSE;
 282
 283     return kTRUE;
 284 }
 285 //----------------------------------------------------------------------------