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