#include <TTree.h>
-#include <TVector.h>
-#include <TObjArray.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TRandom.h>
#include <TArrayI.h>
-#include <TH1.h>
#include <TError.h>
+#include <TH1F.h>
-
+#include "AliLog.h"
#include "AliSTARTDigitizer.h"
#include "AliSTART.h"
#include "AliSTARThit.h"
-#include "AliSTARThitPhoton.h"
#include "AliSTARTdigit.h"
#include "AliRunDigitizer.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
+#include <AliDetector.h>
#include "AliRun.h"
-#include "AliPDG.h"
-#include "AliLoader.h"
-#include "AliRunLoader.h"
-#include "AliSTARTLoader.h"
-
+#include <AliLoader.h>
+#include <AliRunLoader.h>
#include <stdlib.h>
#include <Riostream.h>
#include <Riostream.h>
//___________________________________________
AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
- :AliDigitizer(manager)
+ :AliDigitizer(manager),
+ fSTART(0),
+ fHits(0),
+ fdigits(0),
+ ftimeTDC(0),
+ fADC(0),
+ ftimeTDCAmp(0),
+ fADCAmp(0),
+ fSumMult(0),
+ fEff(0)
{
// cout<<"AliSTARTDigitizer::AliSTARTDigitizer"<<endl;
// ctor which should be used
-// fDebug =0;
- // if (GetDebug()>2)
- // cerr<<"AliSTARTDigitizer::AliSTARTDigitizer"
- // <<"(AliRunDigitizer* manager) was processed"<<endl;
+ AliDebug(1,"processed");
+
+ fSTART = 0;
+ fHits = 0;
+ fdigits = 0;
+
+ ftimeTDC = new TArrayI(24);
+ fADC = new TArrayI(24);
+ ftimeTDCAmp = new TArrayI(24);
+ fADCAmp = new TArrayI(24);
+ fSumMult = new TArrayI(6);
+
+
+ TFile* file = TFile::Open("$ALICE_ROOT/START/PMTefficiency.root");
+ fEff = (TH1F*) file->Get("hEff")->Clone();
+ fEff->SetDirectory(NULL);
+ file->Close();
+ delete file;
}
//------------------------------------------------------------------------
{
// Destructor
+ AliDebug(1,"START");
+ delete ftimeTDC;
+ delete fADC;
+ delete fEff;
+ delete ftimeTDCAmp;
+ delete fADCAmp;
+ delete fSumMult;
}
//------------------------------------------------------------------------
Bool_t AliSTARTDigitizer::Init()
{
// Initialization
-// cout<<"AliSTARTDigitizer::Init"<<endl;
+ AliDebug(1," Init");
return kTRUE;
}
void AliSTARTDigitizer::Exec(Option_t* /*option*/)
{
+ /*
+ Produde digits from hits
+ digits is TObject and includes
+ We are writing array if left & right TDC
+ left & right ADC (will need for slow simulation)
+ TOF first particle left & right
+ mean time and time difference (vertex position)
+
+ */
- AliRunLoader *inRL, *outRL;//in and out Run Loaders
- AliLoader *ingime, *outgime;// in and out ITSLoaders
-
- outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
- outgime = outRL->GetLoader("STARTLoader");
+ //output loader
+ AliRunLoader *outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+ AliLoader * pOutStartLoader = outRL->GetLoader("STARTLoader");
+ AliSTART *fSTART = (AliSTART*)outRL ->GetAliRun()->GetDetector("START");
-#ifdef DEBUG
- cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
-#endif
+ AliDebug(1,"start...");
+ cout<<" AliSTARTDigitizer::Exec "<<endl;
+ //input loader
//
// From hits to digits
//
Int_t hit, nhits;
- Int_t CountEr[13],CountEl[13]; //!!!
- Int_t volume,pmt,tr,tl,sumRight;
- Float_t timediff,timeav;
- Float_t besttimeright,besttimeleft,meanTime;
- Int_t bestRightADC,bestLeftADC;
- Float_t besttimeleftGaus, besttimerightGaus;
- Float_t timeright[13]={13*0};
- Float_t timeleft[13]={13*0};
- Float_t channelWidth=2.5; //ps
- Int_t channelWidthADC=1; //ps
- // Int_t thresholdAmpl=10;
+ Int_t countE[24], sumMult[3];
+ Int_t volume,pmt,tr,qt,qtAmp;
+ Int_t bestRightTDC,bestLeftTDC;
+ Float_t time[24],besttime[24], timeGaus[24] ;
+ Float_t channelWidth=25.; //ps
+ //Q->T-> coefficients !!!! should be asked!!!
+ Float_t ph2mV = 150./500.;
+ Int_t threshold =50; //photoelectrons
+ Int_t mV2channel=200000/(25*25); //5V -> 200ns
- ftimeRightTDC = new TArrayI(12);
- ftimeLeftTDC = new TArrayI(12);
- fRightADC = new TArrayI(12);
- fLeftADC = new TArrayI(12);
- inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
- inRL->LoadgAlice();
-
- // fHits = new TClonesArray ("AliSTARThit", 1000);
- fPhotons = new TClonesArray ("AliSTARThitPhoton", 10000); //!!!
- AliSTART *START = (AliSTART*) gAlice->GetDetector("START");
AliSTARThit *startHit;
- //use if Cherenkov photons
- // AliSTARThitPhoton *startHitPhoton; //!!!
TBranch *brHits=0;
- TBranch *brHitPhoton=0;
- fdigits= new AliSTARTdigit();
Int_t nFiles=fManager->GetNinputs();
for (Int_t inputFile=0; inputFile<nFiles; inputFile++) {
-
- besttimeright=9999.;
- besttimeleft=9999.;
- Int_t timeDiff=0;
- Int_t timeAv=0;
- sumRight=0;
- for (Int_t i0=0; i0<13; i0++)
+ if (inputFile < nFiles-1) {
+ AliWarning(Form("ignoring input stream %d", inputFile));
+ continue;
+
+ }
+
+ Float_t besttimeright=99999.;
+ Float_t besttimeleft=99999.;
+ Int_t timeDiff, meanTime;
+
+ for (Int_t i0=0; i0<24; i0++)
{
- timeright[i0]=0; timeleft[i0]=0;
- CountEr[i0]=0; CountEl[i0]=0;
+ time[i0]=99999; besttime[i0]=99999; countE[i0]=0;
}
+ for ( Int_t imu=0; imu<3; imu++) sumMult[imu]=0;
+ AliRunLoader * inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+ AliLoader * pInStartLoader = inRL->GetLoader("STARTLoader");
+ if (!inRL->GetAliRun()) inRL->LoadgAlice();
- inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
- ingime = inRL->GetLoader("STARTLoader");
- ingime->LoadHits("READ");//probably it is necessary to load them before
- outgime->LoadDigits("UPDATE");//probably it is necessary to load them before
- //use if Cherenkov photons
- // TClonesArray *STARThitsPhotons = START->Photons ();
- TClonesArray *fHits = START->Hits ();
- // cout<<" Load "<<AliSTARTLoader::LoadDigits()<<endl;
-
- TTree *th = ingime->TreeH();
+ //read Hits
+ pInStartLoader->LoadHits("READ");//probably it is necessary to load them before
+ TClonesArray *fHits = fSTART->Hits ();
+ TTree *th = pInStartLoader->TreeH();
brHits = th->GetBranch("START");
- brHitPhoton = th->GetBranch("STARThitPhoton");
if (brHits) {
- START->SetHitsAddressBranch(brHits,brHitPhoton);
+ fSTART->SetHitsAddressBranch(brHits);
}else{
- cerr<<"EXEC Branch START hit not found"<<endl;
+ AliError("Branch START hit not found");
exit(111);
}
Int_t ntracks = (Int_t) th->GetEntries();
- cout<<" ntracks "<<ntracks<<endl;
+
if (ntracks<=0) return;
- // Start loop on tracks in the photon hits containers
- // for amplitude
- /*
- if(brHitPhoton) {
- cout<<"brHitPhoton "<<endl;
- for (Int_t track=0; track<ntracks;track++) {
- brHitPhoton -> GetEntry(track);;
- nhits = STARThitsPhotons->GetEntriesFast();
- for (hit=0;hit<nhits;hit++) {
- startHitPhoton = (AliSTARThitPhoton*)
- STARThitsPhotons ->UncheckedAt(hit);
- pmt=startHitPhoton->fPmt;
- volume = startHitPhoton->fArray;
- if(RegisterPhotoE(startHitPhoton))
- {
- if (volume == 1) CountEr[pmt]++;
- if (volume == 2) CountEl[pmt]++;
- }
- } //hit photons
- } //track photons
- } // was photons
- */
// Start loop on tracks in the hits containers
for (Int_t track=0; track<ntracks;track++) {
brHits->GetEntry(track);
nhits = fHits->GetEntriesFast();
- // cout<<" brHits hits "<<nhits<<endl;
- for (hit=0;hit<nhits;hit++) {
- startHit = (AliSTARThit*) fHits->UncheckedAt(hit);
- if (!startHit) {
- ::Error("Exec","The unchecked hit doesn't exist");
- break;
- }
- pmt=startHit->Pmt();
- volume = startHit->Volume();
- if(volume==1){
- timeright[pmt] = startHit->Time();
- if(timeright[pmt]<besttimeright)
- //&&CountEr[pmt-1]>thresholdAmpl)
- {
- besttimeright=timeright[pmt];
- } //timeright
- }//time for right shoulder
- if(volume==2){
- timeleft[pmt] = startHit->Time();
- if(timeleft[pmt]<besttimeleft)
- //&&CountEl[pmt-1]>thresholdAmpl)
- {
- besttimeleft=timeleft[pmt];
-
- } //timeleftbest
- }//time for left shoulder
- } //hit loop
+ for (hit=0;hit<nhits;hit++)
+ {
+ startHit = (AliSTARThit*) fHits->UncheckedAt(hit);
+ if (!startHit) {
+ AliError("The unchecked hit doesn't exist");
+ break;
+ }
+ pmt=startHit->Pmt();
+ Int_t numpmt=pmt-1;
+ Double_t e=startHit->Etot();
+ volume = startHit->Volume();
+
+ if(e>0 && RegisterPhotoE(e)) {
+ countE[numpmt]++;
+ besttime[numpmt] = startHit->Time();
+ if(besttime[numpmt]<time[numpmt])
+ {
+ time[numpmt]=besttime[numpmt];
+ }
+ }
+ } //hits loop
} //track loop
-
- // z position
- cout<<" right time "<<besttimeright<<
- " right distance "<<besttimeright*30<<endl;;
- cout<<" left time "<<besttimeleft<<
- " left distance "<<besttimeleft*30<<endl;;
-
-
- //folding with experimental time distribution
- besttimeleftGaus=gRandom->Gaus(besttimeright,0.05);
- cout<<" besttimeleftGaus "<<besttimeleftGaus<<endl;
- bestLeftADC=Int_t (besttimeleftGaus*1000/channelWidth);
- Float_t koef=69.7/350.;
- besttimeright=koef*besttimeleft;
- besttimerightGaus=gRandom->Gaus(besttimeleft,0.05);
-
- bestRightADC=Int_t (besttimerightGaus*1000/channelWidth);
- timediff=besttimerightGaus-besttimeleftGaus;
- cout<<" timediff in ns "<<timediff<<" z= "<<timediff*30<<endl;
- meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
- if ( TMath::Abs(timediff)<TMath::Abs(0.3) )
+ //spread time right&left by 25ps && besttime
+ for (Int_t ipmt=0; ipmt<12; ipmt++){
+ if(countE[ipmt] > threshold) {
+ timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25);
+ if(timeGaus[ipmt]<besttimeleft) besttimeleft=timeGaus[ipmt]; //timeleft
+ sumMult[0] += countE[ipmt];
+ sumMult[1] += countE[ipmt];
+ }
+ }
+ for ( Int_t ipmt=12; ipmt<24; ipmt++){
+ if(countE[ipmt] > threshold) {
+ timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25);
+ if(timeGaus[ipmt]<besttimeright) besttimeright=timeGaus[ipmt]; //timeright
+ sumMult[0] += countE[ipmt];
+ sumMult[2] += countE[ipmt];
+ }
+ }
+ //ADC features fill digits
+ //folding with experimental time distribution
+ Float_t c = 0.0299792; // cm/ps
+ Float_t koef=(350.-69.7)/c;
+ Float_t besttimeleftR= besttimeleft;
+ besttimeleft=koef+besttimeleft;
+ bestLeftTDC=Int_t (besttimeleftR/channelWidth);
+ bestRightTDC=Int_t (besttimeright/channelWidth);
+ timeDiff=Int_t ((besttimeright-besttimeleftR)/channelWidth);
+ meanTime=Int_t (((besttimeright+besttimeleft)/2.)/channelWidth);
+ AliDebug(2,Form(" time in ns %f ", besttimeleft));
+ for (Int_t i=0; i<24; i++)
{
- Float_t t1=1000.*besttimeleftGaus;
- Float_t t2=1000.*besttimerightGaus;
- t1=t1/channelWidth; //time in ps to channelWidth
- t2=t2/channelWidth; //time in ps to channelWidth
- timeav=(t1+t2)/2.;
-
- // Time to TDC signal
- // 256 channels for timediff, range 1ns
+ // fill TDC
+ tr= Int_t (timeGaus[i]/channelWidth);
+ if(timeGaus[i]>50000) tr=0;
- timediff=512+1000*timediff/channelWidth; // time in ps
-
- timeAv = (Int_t)(timeav); // time channel numbres
- timeDiff = (Int_t)(timediff); // time channel numbres
- // fill digits
- fdigits->SetTimeBestLeft(bestLeftADC);
- fdigits->SetTimeBestRight(bestRightADC);
- fdigits->SetMeanTime(timeAv);
- fdigits->SetTimeDiff(timeDiff);
- for (Int_t i=0; i<12; i++)
- {
- // fill TDC
- timeright[i+1]=gRandom->Gaus(timeright[i+1],0.05);
- timeleft[i+1]=gRandom->Gaus(timeleft[i+1],0.05);
- tr= Int_t (timeright[i+1]*1000/channelWidth);
- if(tr<200) tr=0;
- tl= Int_t (timeleft[i+1]*1000/channelWidth);
- if(tl<1000) tl=0;
-
- ftimeRightTDC->AddAt(tr,i);
- ftimeLeftTDC->AddAt(tl,i);
- //fill ADC
- Int_t al=( Int_t ) CountEl[i+1]/ channelWidthADC;
- Int_t ar=( Int_t ) CountEr[i+1]/ channelWidthADC;
- fRightADC->AddAt(ar,i);
- fLeftADC ->AddAt(al,i);
- sumRight+=CountEr[i+1];
- }
- fdigits->SetTimeRight(*ftimeRightTDC);
- fdigits->SetTimeLeft(*ftimeLeftTDC);
- fdigits->SetADCRight(*fRightADC);
- fdigits->SetADCLeft(*fLeftADC);
- // cout<<" before sum"<<endl;
- fdigits->SetSumADCRight(sumRight);
- }
- else
- {timeAv=999999; timeDiff=99999;}
-
-// trick to find out output dir:
-
- Char_t nameDigits[20];
- sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
- TDirectory *wd = gDirectory;
- outgime->GetDigitsDataLoader()->GetDirectory()->cd();
- fdigits->Write(nameDigits);
- wd->cd();
+ //fill ADC
+ Int_t al= countE[i];
+ // QTC procedure:
+ // phe -> mV 0.3; 1MIP ->500phe -> ln (amp (mV)) = 5;
+ // max 200ns, HIJING mean 50000phe -> 15000mv -> ln = 15 (s zapasom)
+ // channel 25ps
+ if (al>threshold) {
+ qt=Int_t (TMath::Log(al*ph2mV) * mV2channel);
+ qtAmp=Int_t (TMath::Log(al*10*ph2mV) * mV2channel);
+ cout<<i<<" "<<qt<<" "<<qtAmp<<endl;
+ fADC->AddAt(qt,i);
+ ftimeTDC->AddAt(tr,i);
+ fADCAmp->AddAt(qtAmp,i);
+ ftimeTDCAmp->AddAt(tr,i); //now is the same as non-amplified but will be change
+ }
+ } //pmt loop
+ for (Int_t im=0; im<3; im++)
+ {
+ if (sumMult[im]>threshold){
+ Int_t sum=Int_t (TMath::Log(sumMult[im]*ph2mV) * mV2channel);
+ Int_t sumAmp=Int_t (TMath::Log(10*sumMult[im]*ph2mV) * mV2channel);
+ fSumMult->AddAt(sum,im);
+ fSumMult->AddAt(sumAmp,im+3);
+ }
+ }
- // outgime->WriteDigits("OVERWRITE");
- }
+ fSTART->AddDigit(bestRightTDC,bestLeftTDC,meanTime,timeDiff,fSumMult,
+ ftimeTDC,fADC,ftimeTDCAmp,fADCAmp);
+ pOutStartLoader->UnloadHits();
+ } //input streams loop
+
+ //load digits
+ pOutStartLoader->LoadDigits("UPDATE");
+ TTree *treeD = pOutStartLoader->TreeD();
+ if (treeD == 0x0) {
+ pOutStartLoader->MakeTree("D");
+ treeD = pOutStartLoader->TreeD();
+ }
+ fSTART->MakeBranch("D");
+ treeD->Reset();
+ treeD->Fill();
+
+ pOutStartLoader->WriteDigits("OVERWRITE");
+
+ fSTART->ResetDigits();
+ pOutStartLoader->UnloadDigits();
+
}
//------------------------------------------------------------------------
-Bool_t AliSTARTDigitizer::RegisterPhotoE(/*AliSTARThitPhoton *hit*/)
+Bool_t AliSTARTDigitizer::RegisterPhotoE(Double_t energy)
{
- Double_t P = 0.2;
- Double_t p;
-
- p = gRandom->Rndm();
- if (p > P)
- return kFALSE;
-
- return kTRUE;
+
+
+ // Float_t hc=197.326960*1.e6; //mev*nm
+ Double_t hc=1.973*1.e-6; //gev*nm
+ // cout<<"AliSTARTDigitizer::RegisterPhotoE >> energy "<<energy<<endl;
+ Float_t lambda=hc/energy;
+ Int_t bin= fEff->GetXaxis()->FindBin(lambda);
+ Float_t eff=fEff->GetBinContent(bin);
+ Double_t p = gRandom->Rndm();
+ if (p > eff)
+ return kFALSE;
+
+ return kTRUE;
}
//----------------------------------------------------------------------------
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff