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