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 #include <TH1F.h>
  25
  26 #include "AliLog.h"
  27 #include "AliSTARTDigitizer.h"
  28 #include "AliSTART.h"
  29 #include "AliSTARThit.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    fSTART(0),
  53    fHits(0),
  54    fdigits(0),
  55    ftimeTDC(0),
  56    fADC(0),
  57    ftimeTDCAmp(0),
  58    fADCAmp(0),
  59    fSumMult(0),
  60    fEff(0)
  61 {
  62 // ctor which should be used
  63
  64   AliDebug(1,"processed");
  65
  66   fSTART = 0;
  67   fHits = 0;
  68   fdigits = 0;
  69
  70   ftimeTDC = new TArrayI(24);
  71   fADC = new TArrayI(24);
  72   ftimeTDCAmp = new TArrayI(24);
  73   fADCAmp = new TArrayI(24);
  74   fSumMult = new TArrayI(6);
  75
  76
  77   TFile* file = TFile::Open("$ALICE_ROOT/START/PMTefficiency.root");
  78   fEff = (TH1F*) file->Get("hEff")->Clone();
  79   fEff->SetDirectory(NULL);
  80   file->Close();
  81   delete file;
  82 }
  83
  84 //------------------------------------------------------------------------
  85 AliSTARTDigitizer::~AliSTARTDigitizer()
  86 {
  87 // Destructor
  88
  89   AliDebug(1,"START");
  90
  91   delete ftimeTDC;
  92   delete fADC;
  93   delete fEff;
  94   delete ftimeTDCAmp;
  95   delete  fADCAmp;
  96   delete fSumMult;
  97 }
  98
  99  //------------------------------------------------------------------------
 100 Bool_t AliSTARTDigitizer::Init()
 101 {
 102 // Initialization
 103   AliDebug(1," Init");
 104  return kTRUE;
 105 }
 106
 107
 108 //---------------------------------------------------------------------
 109
 110 void AliSTARTDigitizer::Exec(Option_t* /*option*/)
 111 {
 112
 113   /*
 114     Produde digits from hits
 115         digits is TObject and includes
 116         We are writing array if left & right  TDC
 117         left & right  ADC (will need for slow simulation)
 118         TOF first particle left & right
 119         mean time and time difference (vertex position)
 120
 121   */
 122
 123   //output loader
 124   AliRunLoader *outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 125   AliLoader * pOutStartLoader = outRL->GetLoader("STARTLoader");
 126
 127   AliDebug(1,"start...");
 128   //input loader
 129   //
 130   // From hits to digits
 131   //
 132   Int_t hit, nhits;
 133   Int_t countE[24], sumMult[3];
 134   Int_t volume,pmt,tr,qt,qtAmp;
 135   Int_t  bestRightTDC,bestLeftTDC;
 136   Float_t time[24],besttime[24], timeGaus[24] ;
 137   Float_t channelWidth=25.; //ps
 138     //Q->T-> coefficients !!!! should be asked!!!
 139   Float_t ph2mV = 150./500.;
 140   Int_t threshold =50; //photoelectrons
 141   Int_t mV2channel=200000/(25*25);  //5V -> 200ns
 142
 143
 144   AliSTARThit  *startHit;
 145   TBranch *brHits=0;
 146
 147   Int_t nFiles=fManager->GetNinputs();
 148   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
 149     if (inputFile < nFiles-1) {
 150       AliWarning(Form("ignoring input stream %d", inputFile));
 151       continue;
 152
 153     }
 154
 155     Float_t besttimeright=99999.;
 156     Float_t besttimeleft=99999.;
 157     Int_t timeDiff, meanTime;
 158
 159     for (Int_t i0=0; i0<24; i0++)
 160       {
 161         time[i0]=besttime[i0]=timeGaus[i0]=99999; countE[i0]=0;
 162       }
 163     for ( Int_t imu=0; imu<3; imu++) sumMult[imu]=0;
 164     AliRunLoader * inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
 165     AliLoader * pInStartLoader = inRL->GetLoader("STARTLoader");
 166     if (!inRL->GetAliRun()) inRL->LoadgAlice();
 167     AliSTART *fSTART  = (AliSTART*)inRL ->GetAliRun()->GetDetector("START");
 168
 169        //read Hits
 170     pInStartLoader->LoadHits("READ");//probably it is necessary to load them before
 171     TClonesArray *fHits = fSTART->Hits ();
 172     TTree *th = pInStartLoader->TreeH();
 173     brHits = th->GetBranch("START");
 174     if (brHits) {
 175       fSTART->SetHitsAddressBranch(brHits);
 176     }else{
 177       AliError("Branch START hit not found");
 178       exit(111);
 179     }
 180     Int_t ntracks    = (Int_t) th->GetEntries();
 181
 182     if (ntracks<=0) return;
 183     // Start loop on tracks in the hits containers
 184     for (Int_t track=0; track<ntracks;track++) {
 185       brHits->GetEntry(track);
 186       nhits = fHits->GetEntriesFast();
 187       for (hit=0;hit<nhits;hit++)
 188         {
 189           startHit   = (AliSTARThit*) fHits->UncheckedAt(hit);
 190           if (!startHit) {
 191             AliError("The unchecked hit doesn't exist");
 192             break;
 193           }
 194           pmt=startHit->Pmt();
 195           Int_t numpmt=pmt-1;
 196                   Double_t e=startHit->Etot();
 197           volume = startHit->Volume();
 198
 199           if(e>0 && RegisterPhotoE(e)) {
 200             countE[numpmt]++;
 201             besttime[numpmt] = startHit->Time();
 202             if(besttime[numpmt]<time[numpmt])
 203               {
 204                 time[numpmt]=besttime[numpmt];
 205               }
 206           }
 207         } //hits loop
 208     } //track loop
 209
 210     //spread time right&left by 25ps   && besttime
 211     for (Int_t ipmt=0; ipmt<12; ipmt++){
 212       if(countE[ipmt] > threshold) {
 213         timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25);
 214         if(timeGaus[ipmt]<besttimeleft) besttimeleft=timeGaus[ipmt]; //timeleft
 215         sumMult[0] +=  countE[ipmt];
 216         sumMult[1] +=  countE[ipmt];
 217       }
 218     }
 219     for ( Int_t ipmt=12; ipmt<24; ipmt++){
 220       if(countE[ipmt] > threshold) {
 221         timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25);
 222         if(timeGaus[ipmt]<besttimeright)  besttimeright=timeGaus[ipmt]; //timeright
 223         sumMult[0] +=  countE[ipmt];
 224         sumMult[2] +=  countE[ipmt];
 225       }
 226     }
 227     //ADC features fill digits
 228    //folding with experimental time distribution
 229     Float_t c = 0.0299792; // cm/ps
 230     Float_t koef=(350.-69.7)/c;
 231     Float_t  besttimeleftR= besttimeleft;
 232     besttimeleft=koef+besttimeleft;
 233     bestLeftTDC=Int_t (besttimeleftR/channelWidth);
 234     bestRightTDC=Int_t (besttimeright/channelWidth);
 235     timeDiff=Int_t ((besttimeright-besttimeleftR)/channelWidth);
 236     meanTime=Int_t (((besttimeright+besttimeleft)/2.)/channelWidth);
 237     AliDebug(2,Form(" time in ns %f ", besttimeleft));
 238     for (Int_t i=0; i<24; i++)
 239       {
 240         //  fill TDC
 241         tr= Int_t (timeGaus[i]/channelWidth);
 242         if(timeGaus[i]>50000) tr=0;
 243
 244         //fill ADC
 245         Int_t al= countE[i];
 246         // QTC procedure:
 247         // phe -> mV 0.3; 1MIP ->500phe -> ln (amp (mV)) = 5;
 248         // max 200ns, HIJING  mean 50000phe -> 15000mv -> ln = 15 (s zapasom)
 249         // channel 25ps
 250         if (al>threshold) {
 251           qt=Int_t (TMath::Log(al*ph2mV) * mV2channel);
 252           qtAmp=Int_t (TMath::Log(al*10*ph2mV) * mV2channel);
 253           fADC->AddAt(qt,i);
 254           ftimeTDC->AddAt(tr,i);
 255           fADCAmp->AddAt(qtAmp,i);
 256           ftimeTDCAmp->AddAt(tr,i); //now is the same as non-amplified but will be change
 257         }
 258       } //pmt loop
 259     for (Int_t im=0; im<3; im++)
 260       {
 261         if (sumMult[im]>threshold){
 262           Int_t sum=Int_t (TMath::Log(sumMult[im]*ph2mV) * mV2channel);
 263           Int_t sumAmp=Int_t (TMath::Log(10*sumMult[im]*ph2mV) * mV2channel);
 264           fSumMult->AddAt(sum,im);
 265           fSumMult->AddAt(sumAmp,im+3);
 266         }
 267       }
 268
 269     fSTART->AddDigit(bestRightTDC,bestLeftTDC,meanTime,timeDiff,fSumMult,
 270                      ftimeTDC,fADC,ftimeTDCAmp,fADCAmp);
 271     pOutStartLoader->UnloadHits();
 272   } //input streams loop
 273
 274     //load digits
 275     pOutStartLoader->LoadDigits("UPDATE");
 276     TTree *treeD  = pOutStartLoader->TreeD();
 277     if (treeD == 0x0) {
 278       pOutStartLoader->MakeTree("D");
 279       treeD = pOutStartLoader->TreeD();
 280     }
 281     AliSTART *fSTART  = (AliSTART*)outRL ->GetAliRun()->GetDetector("START");
 282     fSTART->MakeBranch("D");
 283      treeD->Reset();
 284      treeD->Fill();
 285
 286   pOutStartLoader->WriteDigits("OVERWRITE");
 287
 288   fSTART->ResetDigits();
 289   pOutStartLoader->UnloadDigits();
 290
 291 }
 292
 293
 294 //------------------------------------------------------------------------
 295 Bool_t AliSTARTDigitizer::RegisterPhotoE(Double_t energy)
 296 {
 297
 298
 299   //  Float_t hc=197.326960*1.e6; //mev*nm
 300   Double_t hc=1.973*1.e-6; //gev*nm
 301   Float_t lambda=hc/energy;
 302   Int_t bin=  fEff->GetXaxis()->FindBin(lambda);
 303   Float_t eff=fEff->GetBinContent(bin);
 304   Double_t  p = gRandom->Rndm();
 305   if (p > eff)
 306     return kFALSE;
 307
 308   return kTRUE;
 309 }
 310 //----------------------------------------------------------------------------