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 <TVector.h>
  20 #include <TObjArray.h>
  21 #include <TFile.h>
  22 #include <TDirectory.h>
  23 #include <TRandom.h>
  24 #include <TArrayI.h>
  25 #include <TH1.h>
  26
  27
  28 #include "AliSTARTDigitizer.h"
  29 #include "AliSTART.h"
  30 #include "AliSTARThit.h"
  31 #include "AliSTARThitPhoton.h"
  32 #include "AliSTARTdigit.h"
  33 #include "AliRunDigitizer.h"
  34 #include "AliHeader.h"
  35 #include "AliGenEventHeader.h"
  36 #include "AliRun.h"
  37 #include "AliPDG.h"
  38
  39 #include <stdlib.h>
  40 #include <iostream.h>
  41 #include <fstream.h>
  42
  43 ClassImp(AliSTARTDigitizer)
  44
  45 //___________________________________________
  46   AliSTARTDigitizer::AliSTARTDigitizer()  :AliDigitizer()
  47 {
  48 // Default ctor - don't use it
  49   ;
  50 }
  51
  52 //___________________________________________
  53 AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
  54     :AliDigitizer(manager)
  55 {
  56         cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
  57 // ctor which should be used
  58 //  fDebug =0;
  59  // if (GetDebug()>2)
  60   //  cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
  61    //     <<"(AliRunDigitizer* manager) was processed"<<endl;
  62
  63 }
  64
  65 //------------------------------------------------------------------------
  66 AliSTARTDigitizer::~AliSTARTDigitizer()
  67 {
  68 // Destructor
  69
  70
  71 }
  72
  73  //------------------------------------------------------------------------
  74 Bool_t AliSTARTDigitizer::Init()
  75 {
  76 // Initialization
  77  cout<<"AliSTARTDigitizer::Init"<<endl;
  78  return kTRUE;
  79 }
  80
  81
  82 //---------------------------------------------------------------------
  83
  84 void AliSTARTDigitizer::Exec(Option_t* option)
  85 {
  86
  87   cout<<" AliSTARTDigitizer::Exec"<<endl;
  88 #ifdef DEBUG
  89   cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
  90 #endif
  91   //
  92   // From hits to digits
  93   //
  94   Int_t hit, nhits;
  95   Int_t CountEr[13],CountEl[13];                                                        //!!!
  96   Int_t volume,pmt,tr,tl,sumRight;
  97   char nameDigits[20];
  98   Float_t timediff,timeav;
  99   Float_t besttimeright,besttimeleft,meanTime;
 100   Int_t  bestRightADC,bestLeftADC;
 101   Float_t besttimeleftGaus, besttimerightGaus;
 102   Float_t timeright[13]={13*0};
 103   Float_t timeleft[13]={13*0};
 104   Float_t channelWidth=2.5; //ps
 105   Int_t channelWidthADC=1; //ps
 106   Int_t thresholdAmpl=10;
 107
 108   ftimeRightTDC = new TArrayI(12);
 109   ftimeLeftTDC = new TArrayI(12);
 110   fRightADC = new TArrayI(12);
 111   fLeftADC = new TArrayI(12);
 112
 113
 114   fHits = new TClonesArray ("AliSTARThit", 1000);
 115   fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000);                     //!!!
 116   AliSTART *START  = (AliSTART*) gAlice->GetDetector("START");
 117
 118   AliSTARThit  *startHit;
 119   AliSTARThitPhoton  *startHitPhoton;                                           //!!!
 120   TBranch *brHits=0;
 121   TBranch *brHitPhoton=0;
 122   fdigits= new AliSTARTdigit();
 123
 124   Int_t nFiles=fManager->GetNinputs();
 125   for (Int_t inputFile=0; inputFile<nFiles;  inputFile++) {
 126
 127
 128
 129     besttimeright=9999.;
 130     besttimeleft=9999.;
 131     Int_t timeDiff=0;
 132     Int_t timeAv=0;
 133     sumRight=0;
 134     for (Int_t i0=0; i0<13; i0++)
 135       {
 136         timeright[i0]=0; timeleft[i0]=0;
 137         CountEr[i0]=0;   CountEl[i0]=0;
 138       }
 139     TClonesArray *STARThits = START->Hits ();
 140     TClonesArray *STARThitsPhotons = START->Photons ();
 141
 142    TTree *th = fManager->GetInputTreeH(inputFile);
 143     brHits = th->GetBranch("START");
 144     brHitPhoton = th->GetBranch("STARThitPhoton");
 145     if (brHits) {
 146       START->SetHitsAddressBranch(brHits,brHitPhoton);
 147     }else{
 148       cerr<<"EXEC Branch START hit not found"<<endl;
 149       exit(111);
 150     }
 151     Int_t ntracks    = (Int_t) th->GetEntries();
 152     cout<<" ntracks "<<ntracks<<endl;
 153     if (ntracks<=0) return;
 154     // Start loop on tracks in the photon hits containers
 155     // for amplitude
 156     if(brHitPhoton) {
 157       for (Int_t track=0; track<ntracks;track++) {
 158         brHitPhoton -> GetEntry(track);;
 159         nhits = STARThitsPhotons->GetEntriesFast();
 160         for (hit=0;hit<nhits;hit++) {
 161           startHitPhoton   = (AliSTARThitPhoton*)
 162            STARThitsPhotons ->UncheckedAt(hit);
 163           pmt=startHitPhoton->fPmt;
 164           volume = startHitPhoton->fArray;
 165           if(RegisterPhotoE(startHitPhoton))
 166             {
 167               if (volume == 1) CountEr[pmt]++;
 168               if (volume == 2) CountEl[pmt]++;
 169             }
 170         } //hit photons
 171       } //track photons
 172     } // was photons
 173
 174     // Start loop on tracks in the hits containers
 175     for (Int_t track=0; track<ntracks;track++) {
 176       brHits->GetEntry(track);
 177       nhits = STARThits->GetEntriesFast();
 178       for (hit=0;hit<nhits;hit++) {
 179         startHit   = (AliSTARThit*) STARThits->UncheckedAt(hit);
 180         pmt=startHit->fPmt;
 181         volume = startHit->fVolume;
 182         if(volume==1){
 183           timeright[pmt] = startHit->fTime;
 184           if(timeright[pmt]<besttimeright)
 185             //&&CountEr[pmt-1]>thresholdAmpl)
 186             {
 187             besttimeright=timeright[pmt];
 188           } //timeright
 189         }//time for right shoulder
 190         if(volume==2){
 191           timeleft[pmt] = startHit->fTime;
 192           if(timeleft[pmt]<besttimeleft)
 193             //&&CountEl[pmt-1]>thresholdAmpl)
 194             {
 195             besttimeleft=timeleft[pmt];
 196
 197           } //timeleftbest
 198         }//time for left shoulder
 199       } //hit loop
 200     } //track loop
 201
 202     // z position
 203     cout<<" right time  "<<besttimeright<<
 204       " right distance "<<besttimeright*30<<endl;;
 205     cout<<" left time  "<<besttimeleft<<
 206       " right distance "<<besttimeleft*30<<endl;;
 207
 208
 209     //folding with experimental time distribution
 210
 211     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
 212     //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
 213     bestRightADC=(Int_t) besttimerightGaus*1000/channelWidth;
 214     Float_t koef=69.7/350.;
 215     besttimeleft=koef*besttimeleft;
 216     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
 217
 218     bestLeftADC=(Int_t) besttimeleftGaus*1000/channelWidth;
 219     timediff=besttimerightGaus-besttimeleftGaus;
 220     cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
 221     meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
 222     if ( TMath::Abs(timediff)<TMath::Abs(0.3) )
 223       {
 224         Float_t t1=1000.*besttimeleftGaus;
 225         Float_t t2=1000.*besttimerightGaus;
 226         t1=t1/channelWidth;   //time in ps to channelWidth
 227         t2=t2/channelWidth;   //time in ps to channelWidth
 228         timeav=(t1+t2)/2.;
 229
 230         // Time to TDC signal
 231         // 256 channels for timediff, range 1ns
 232
 233         timediff=512+1000*timediff/channelWidth; // time in ps
 234
 235         timeAv = (Int_t)(timeav);   // time  channel numbres
 236         timeDiff = (Int_t)(timediff); // time  channel numbres
 237         //       fill digits
 238         fdigits->SetTimeBestLeft(bestLeftADC);
 239         fdigits->SetTimeBestRight(bestRightADC);
 240         fdigits->SetMeanTime(timeAv);
 241         fdigits->SetTimeDiff(timeDiff);
 242         for (Int_t i=0; i<12; i++)
 243           {
 244             //  fill TDC
 245             timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
 246             timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
 247             tr= Int_t (timeright[i+1]*1000/channelWidth);
 248             if(tr<200) tr=0;
 249             tl= Int_t (timeleft[i+1]*1000/channelWidth);
 250             if(tl<1000) tl=0;
 251
 252             ftimeRightTDC->AddAt(tr,i);
 253             ftimeLeftTDC->AddAt(tl,i);
 254             //fill ADC
 255             Int_t al=( Int_t ) CountEl[i+1]/ channelWidthADC;
 256             Int_t ar=( Int_t ) CountEr[i+1]/ channelWidthADC;
 257             fRightADC->AddAt(ar,i);
 258             fLeftADC ->AddAt(al,i);
 259             sumRight+=CountEr[i+1];
 260           }
 261         fdigits->SetTimeRight(*ftimeRightTDC);
 262         fdigits->SetTimeLeft(*ftimeLeftTDC);
 263         fdigits->SetADCRight(*fRightADC);
 264         fdigits->SetADCLeft(*fLeftADC);
 265         // cout<<" before sum"<<endl;
 266         fdigits->SetSumADCRight(sumRight);
 267       }
 268     else
 269       {timeAv=999999; timeDiff=99999;}
 270
 271     // trick to find out output dir:
 272     TTree *outTree = fManager->GetTreeD();
 273     if (!outTree) {
 274       cerr<<"something wrong with output...."<<endl;
 275       exit(111);
 276     }
 277
 278     TDirectory *wd = gDirectory;
 279     outTree->GetDirectory()->cd();
 280     sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
 281     fdigits->Write(nameDigits);
 282     cout<<nameDigits<<endl;
 283     wd->cd();
 284   }
 285 }
 286
 287
 288 //------------------------------------------------------------------------
 289 Bool_t AliSTARTDigitizer::RegisterPhotoE(AliSTARThitPhoton *hit)
 290 {
 291     Double_t    P = 0.2;
 292     Double_t    p;
 293
 294     p = gRandom->Rndm();
 295     if (p > P)
 296       return kFALSE;
 297
 298     return kTRUE;
 299 }