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   cout<<" outgime "<<outgime<<endl;
  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   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     //   AliSTARTLoader *pSTARTloader = (AliSTARTLoader*)fLoader;
 150     //  pSTARTLoader->LoadHits("READ");
 151     ingime = inRL->GetLoader("STARTLoader");
 152     cout<<" ingime "<<ingime<<endl;
 153     //    AliSTARTLoader *pSTARTLoader ;
 154     ingime->LoadHits("READ");//probably it is necessary to load them before
 155     ingime->LoadDigits("UPDATE");//probably it is necessary to load them before
 156     TClonesArray *STARThitsPhotons = START->Photons ();
 157     TClonesArray *fHits = START->Hits ();
 158     //    cout<<" Load  "<<AliSTARTLoader::LoadDigits()<<endl;
 159
 160     TTree *th = ingime->TreeH();
 161     brHits = th->GetBranch("START");
 162     brHitPhoton = th->GetBranch("STARThitPhoton");
 163     if (brHits) {
 164       cout<<" brHits "<<endl;
 165       START->SetHitsAddressBranch(brHits,brHitPhoton);
 166     }else{
 167       cerr<<"EXEC Branch START hit not found"<<endl;
 168       exit(111);
 169     }
 170     Int_t ntracks    = (Int_t) th->GetEntries();
 171     cout<<" ntracks "<<ntracks<<endl;
 172     if (ntracks<=0) return;
 173     // Start loop on tracks in the photon hits containers
 174     // for amplitude
 175     /*
 176     if(brHitPhoton) {
 177       cout<<"brHitPhoton "<<endl;
 178       for (Int_t track=0; track<ntracks;track++) {
 179         brHitPhoton -> GetEntry(track);;
 180         nhits = STARThitsPhotons->GetEntriesFast();
 181         for (hit=0;hit<nhits;hit++) {
 182           startHitPhoton   = (AliSTARThitPhoton*)
 183            STARThitsPhotons ->UncheckedAt(hit);
 184           pmt=startHitPhoton->fPmt;
 185           volume = startHitPhoton->fArray;
 186           if(RegisterPhotoE(startHitPhoton))
 187             {
 188               if (volume == 1) CountEr[pmt]++;
 189               if (volume == 2) CountEl[pmt]++;
 190             }
 191         } //hit photons
 192       } //track photons
 193     } // was photons
 194     */
 195     // Start loop on tracks in the hits containers
 196     cout<<" fHits "<<fHits<<endl;;
 197     for (Int_t track=0; track<ntracks;track++) {
 198       brHits->GetEntry(track);
 199       nhits = fHits->GetEntriesFast();
 200       //  cout<<" brHits hits "<<nhits<<endl;
 201       for (hit=0;hit<nhits;hit++) {
 202         startHit   = (AliSTARThit*) fHits->UncheckedAt(hit);
 203         pmt=startHit->fPmt;
 204         volume = startHit->fVolume;
 205         if(volume==1){
 206           timeright[pmt] = startHit->fTime;
 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->fTime;
 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       " right distance "<<besttimeleft*30<<endl;;
 230
 231
 232     //folding with experimental time distribution
 233
 234     besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
 235     //    cout<<" besttimerightGaus "<<besttimerightGaus<<endl;
 236     bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
 237     Float_t koef=69.7/350.;
 238     besttimeleft=koef*besttimeleft;
 239     besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
 240
 241     bestLeftADC=Int_t (besttimeleftGaus*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
 297 /*
 298    // trick to find out output dir:
 299     TTree *outTree = fManager->GetTreeD();
 300     if (!outTree) {
 301       cerr<<"something wrong with output...."<<endl;
 302       exit(111);
 303     }
 304
 305     Char_t nameDigits[20];
 306     TDirectory *wd = gDirectory;
 307     outTree->GetDirectory()->cd();
 308     fdigits->Write(nameDigits);
 309     cout<<nameDigits<<endl;
 310     wd->cd();
 311 */
 312     cout<<" outgime v konce "<<outgime<<endl;
 313     outgime->Dump();
 314     outgime->Print();
 315     Char_t nameDigits[20];
 316     sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
 317     fdigits->Write(nameDigits);
 318
 319     //    outgime->WriteDigits("OVERWRITE");
 320   }
 321 }
 322
 323
 324 //------------------------------------------------------------------------
 325 Bool_t AliSTARTDigitizer::RegisterPhotoE(AliSTARThitPhoton *hit)
 326 {
 327     Double_t    P = 0.2;
 328     Double_t    p;
 329
 330     p = gRandom->Rndm();
 331     if (p > P)
 332       return kFALSE;
 333
 334     return kTRUE;
 335 }
 336 //----------------------------------------------------------------------------