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