#include <TTree.h>
-#include <TVector.h>
-#include <TObjArray.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TRandom.h>
+#include <TArrayI.h>
+#include <TError.h>
+#include <TH1F.h>
+#include <TGraph.h>
-
+#include "AliLog.h"
#include "AliSTARTDigitizer.h"
#include "AliSTART.h"
#include "AliSTARThit.h"
#include "AliSTARTdigit.h"
#include "AliRunDigitizer.h"
-
+#include <AliDetector.h>
#include "AliRun.h"
-#include "AliPDG.h"
-
+#include <AliLoader.h>
+#include <AliRunLoader.h>
#include <stdlib.h>
-#include <iostream.h>
-#include <fstream.h>
+#include <Riostream.h>
+#include <Riostream.h>
+#include "AliSTARTParameters.h"
+#include "AliCDBLocal.h"
+#include "AliCDBStorage.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
ClassImp(AliSTARTDigitizer)
//___________________________________________
AliSTARTDigitizer::AliSTARTDigitizer(AliRunDigitizer* manager)
- :AliDigitizer(manager)
+ :AliDigitizer(manager),
+ fSTART(0),
+ fHits(0),
+ fdigits(0),
+ ftimeCFD(0),
+ ftimeLED(0),
+ fADC(0),
+ fADC0(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;
+
+ ftimeCFD = new TArrayI(24);
+ fADC = new TArrayI(24);
+ ftimeLED = new TArrayI(24);
+ fADC0 = new TArrayI(24);
+
}
{
// Destructor
+ AliDebug(1,"START");
+ delete ftimeCFD;
+ delete fADC;
+ delete ftimeLED;
+ delete fADC0;
}
//------------------------------------------------------------------------
Bool_t AliSTARTDigitizer::Init()
{
// Initialization
- cout<<"AliSTARTDigitizer::Init"<<endl;
+ AliDebug(1," Init");
return kTRUE;
}
//---------------------------------------------------------------------
-void AliSTARTDigitizer::Exec(Option_t* option)
+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)
+
+ */
+ //output loader
+ AliRunLoader *outRL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+ AliLoader * pOutStartLoader = outRL->GetLoader("STARTLoader");
-#ifdef DEBUG
- cout<<"AliSTARTDigitizer::>SDigits2Digits start...\n";
-#endif
+ AliDebug(1,"start...");
+ //input loader
//
// From hits to digits
//
- Int_t hit;
- Int_t nhits;
- Int_t volume,pmt;
- char nameDigits[20];
- Float_t timediff,timeright,timeleft,timeav;
- Float_t besttimeright,besttimeleft,meanTime;
- Int_t channelWidth=10;
- fHits = new TClonesArray ("AliSTARThit", 1000);
- AliSTART *START = (AliSTART*) gAlice->GetDetector("START");
+ Int_t hit, nhits;
+ Int_t countE[24];
+ Int_t volume, pmt, trCFD, trLED;
+ Float_t sl, qt;
+ Int_t bestRightTDC, bestLeftTDC, qtCh;
+ Float_t time[24], besttime[24], timeGaus[24] ;
+ //Q->T-> coefficients !!!! should be asked!!!
+ Float_t gain[24],timeDelayCFD[24], timeDelayLED[24];
+ Int_t threshold =50; //photoelectrons
+ Float_t zdetA, zdetC;
+ Int_t sumMultCoeff = 100;
+ TObjArray slewingLED;
+ TObjArray slewingRec;
+ AliSTARTParameters* param = AliSTARTParameters::Instance();
+ param->Init();
+
+ Int_t ph2Mip = param->GetPh2Mip();
+ Int_t channelWidth = param->GetChannelWidth() ;
+ Float_t delayVertex = param->GetTimeDelayTVD();
+ for (Int_t i=0; i<24; i++){
+ timeDelayCFD[i] = param->GetTimeDelayCFD(i);
+ timeDelayLED[i] = param->GetTimeDelayLED(i);
+ gain[i] = param->GetGain(i);
+ TGraph* gr = param ->GetSlew(i);
+ slewingLED.AddAtAndExpand(gr,i);
+
+ TGraph* gr1 = param ->GetSlewRec(i);
+ slewingRec.AddAtAndExpand(gr1,i);
+ TGraph* grEff = param ->GetPMTeff(i);
+ fEffPMT.AddAtAndExpand(grEff,i);
+ }
+ zdetC = param->GetZposition(0);
+ zdetA = param->GetZposition(1);
+
AliSTARThit *startHit;
TBranch *brHits=0;
- fdigits= new AliSTARTdigit();
-
+
Int_t nFiles=fManager->GetNinputs();
for (Int_t inputFile=0; inputFile<nFiles; inputFile++) {
- sprintf(nameDigits,"START_D_%d",fManager->GetOutputEventNr());
+ if (inputFile < nFiles-1) {
+ AliWarning(Form("ignoring input stream %d", inputFile));
+ continue;
+
+ }
+
+ Float_t besttimeright=99999.;
+ Float_t besttimeleft=99999.;
+ Int_t pmtBestRight=9999;
+ Int_t pmtBestLeft=9999;
+ Int_t timeDiff=999, meanTime=0;
+ Int_t sumMult =0; fSumMult=0;
+ bestRightTDC = 99999; bestLeftTDC = 99999;
- besttimeright=9999.;
- besttimeleft=9999.;
- Int_t timeDiff=0;
- Int_t timeAv=0;
- TClonesArray *STARThits = START->Hits ();
+ ftimeCFD -> Reset();
+ fADC -> Reset();
+ fADC0 -> Reset();
+ ftimeLED ->Reset();
+ for (Int_t i0=0; i0<24; i0++)
+ {
+ time[i0]=besttime[i0]=timeGaus[i0]=99999; countE[i0]=0;
+ }
+ AliRunLoader * inRL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+ AliLoader * pInStartLoader = inRL->GetLoader("STARTLoader");
+ if (!inRL->GetAliRun()) inRL->LoadgAlice();
+ AliSTART *fSTART = (AliSTART*)inRL ->GetAliRun()->GetDetector("START");
- TTree *th = fManager->GetInputTreeH(inputFile);
+ //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");
if (brHits) {
- START->SetHitsAddressBranch(brHits);
+ fSTART->SetHitsAddressBranch(brHits);
}else{
- cerr<<"EXEC Branch START hit not found"<<exit;
+ AliError("Branch START hit not found");
+ exit(111);
}
Int_t ntracks = (Int_t) th->GetEntries();
+
if (ntracks<=0) return;
// Start loop on tracks in the hits containers
for (Int_t track=0; track<ntracks;track++) {
brHits->GetEntry(track);
- nhits = STARThits->GetEntriesFast();
- for (hit=0;hit<nhits;hit++) {
- startHit = (AliSTARThit*) STARThits->UncheckedAt(hit);
- pmt=startHit->fPmt;
- volume = startHit->fVolume;
- if(volume==1){
- timeright = startHit->fTime;
- if(timeright<besttimeright) {
- besttimeright=timeright;
- } //timeright
- }//time for right shoulder
- if(volume==2){
- timeleft = startHit->fTime;
- if(timeleft<besttimeleft) {
- besttimeleft=timeleft;
- } //timeleftbest
- }//time for left shoulder
- } //hit loop
+ nhits = fHits->GetEntriesFast();
+ 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(numpmt,e)) {
+ if(e>0 ) {
+ countE[numpmt]++;
+ besttime[numpmt] = startHit->Time();
+ if(besttime[numpmt]<time[numpmt])
+ {
+ time[numpmt]=besttime[numpmt];
+ }
+ } //photoelectron accept
+ } //hits loop
} //track loop
-
- //folding with experimental time distribution
- Float_t besttimerightGaus=gRandom->Gaus(besttimeright,0.05);
- Float_t koef=69.7/350.;
- besttimeleft=koef*besttimeleft;
- Float_t besttimeleftGaus=gRandom->Gaus(besttimeleft,0.05);
- timediff=besttimerightGaus-besttimeleftGaus;
- meanTime=(besttimerightGaus+besttimeleftGaus)/2.;
- if ( TMath::Abs(timediff)<TMath::Abs(3.) && meanTime<TMath::Abs(5.))
- {
- //we assume centre of bunch is 5ns after TTS signal
- //TOF values are relative of the end of bunch
- Float_t ppBunch=25;
- ppBunch=ppBunch-10/2;
- Float_t t1=1000.*besttimeleftGaus;
- Float_t t2=1000.*besttimerightGaus;
- t1=t1/channelWidth+ppBunch; //time in ps to channelWidth
- t2=t2/channelWidth+ppBunch; //time in ps to channelWidth
- timeav=(t1+t2)/2.;
-
- // Time to TDC signal
- // 256 channels for timediff, range 1ns
-
- timediff=128+1000*timediff/channelWidth; // time in ps
-
- timeAv = (Int_t)(timeav); // time (ps) channel numbres
- timeDiff = (Int_t)(timediff); // time ( ps) channel numbres
- fdigits->Set(timeAv,timeDiff);
- fdigits->Print();
+ //spread time right&left by 25ps && besttime
+ Float_t c = 0.0299792; // cm/ps
+
+ Float_t koef=(zdetA-zdetC)/c; //correction position difference by cable
+ for (Int_t ipmt=0; ipmt<12; ipmt++){
+ if(countE[ipmt] > threshold) {
+ timeGaus[ipmt]=gRandom->Gaus(time[ipmt],25)+koef;
+ if(timeGaus[ipmt]<besttimeleft){
+ besttimeleft=timeGaus[ipmt]; //timeleft
+ pmtBestLeft=ipmt;}
}
- else
- {timeAv=999999; timeDiff=99999;}
-
-// trick to find out output dir:
- TTree *outTree = fManager->GetTreeD();
- if (!outTree) {
- cerr<<"something wrong with output...."<<exit;
- return;
}
- TDirectory *wd = gDirectory;
- outTree->GetDirectory()->cd();
- fdigits->Write(nameDigits);
- wd->cd();
- }
+ 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
+ pmtBestRight=ipmt;}
+ }
+ }
+ //folding with alignmentz position distribution
+ if( besttimeleft > 10000. && besttimeleft <15000)
+ bestLeftTDC=Int_t ((besttimeleft+1000*timeDelayCFD[pmtBestLeft])
+ /channelWidth);
+
+ if( besttimeright > 10000. && besttimeright <15000)
+ bestRightTDC=Int_t ((besttimeright+1000*timeDelayCFD[pmtBestRight])
+ /channelWidth);
+
+ if (bestRightTDC < 99999 && bestLeftTDC < 99999)
+ {
+ timeDiff=Int_t (((besttimeleft-besttimeright)+1000*delayVertex)
+ /channelWidth);
+ meanTime=Int_t (((besttimeright+1000*timeDelayCFD[pmtBestLeft]+
+ besttimeleft+1000*timeDelayCFD[pmtBestLeft])/2.)
+ /channelWidth);
+ }
+ AliDebug(10,Form(" time right& left %i %i time diff && mean time in channels %i %i",bestRightTDC,bestLeftTDC, timeDiff, meanTime));
+ for (Int_t i=0; i<24; i++)
+ {
+ Float_t al = countE[i];
+ if (al>threshold && timeGaus[i]<50000 ) {
+ //fill ADC
+ // QTC procedure:
+ // phe -> mV 0.3; 1MIP ->500phe -> ln (amp (mV)) = 5;
+ // max 200ns, HIJING mean 50000phe -> 15000mv -> ln = 15 (s zapasom)
+ // channel 25ps
+ qt= 50.*al*gain[i]/ph2Mip; // 50mv/Mip amp in mV
+ // fill TDC
+ trCFD = Int_t (timeGaus[i] + 1000.*timeDelayCFD[i])/channelWidth;
+ trLED= Int_t (timeGaus[i] + 1000.*timeDelayLED[i]);
+ sl = ((TGraph*)slewingLED.At(i))->Eval(qt);
+ trLED = Int_t(( trLED + 1000*sl )/channelWidth);
+ qtCh=Int_t (1000.*TMath::Log(qt)) / channelWidth;
+ fADC0->AddAt(0,i);
+ fADC->AddAt(qtCh,i);
+ ftimeCFD->AddAt(Int_t (trCFD),i);
+ ftimeLED->AddAt(trLED,i);
+ // sumMult += Int_t ((al*gain[i]/ph2Mip)*50) ;
+ sumMult += Int_t (qt/sumMultCoeff) ;
+
+ AliDebug(10,Form(" pmt %i : time in ns %f time in channels %i ",
+ i, timeGaus[i],trCFD ));
+ AliDebug(10,Form(" qt in mV %f qt in ns %f qt in channels %i ",qt,
+ TMath::Log(qt), qtCh));
+ }
+ } //pmt loop
+ if (sumMult > threshold){
+ fSumMult = Int_t (1000.* TMath::Log(Double_t(sumMult) / Double_t(sumMultCoeff))
+ /channelWidth);
+ AliDebug(10,Form("summult mv %i mult in chammens %i in ps %i ",
+ sumMult, fSumMult, fSumMult*channelWidth));
+ }
+ // if ( besttimeright<99999 || besttimeleft < 99999) {
+
+ fSTART->AddDigit(bestRightTDC,bestLeftTDC,meanTime,timeDiff,fSumMult,
+ ftimeCFD,fADC,ftimeLED,fADC0);
+ // }
+
+ AliDebug(10,Form(" Digits wrote bestRightTDC %i bestLeftTDC %i meanTime %i timeDiff %i fSumMult %i ", bestRightTDC,bestLeftTDC,meanTime,timeDiff,fSumMult ));
+ pOutStartLoader->UnloadHits();
+ } //input streams loop
+
+ //load digits
+ pOutStartLoader->LoadDigits("UPDATE");
+ TTree *treeD = pOutStartLoader->TreeD();
+ if (treeD == 0x0) {
+ pOutStartLoader->MakeTree("D");
+ treeD = pOutStartLoader->TreeD();
+ }
+ treeD->Reset();
+ fSTART = (AliSTART*)outRL ->GetAliRun()->GetDetector("START");
+ // Make a branch in the tree
+ fSTART->MakeBranch("D");
+ treeD->Fill();
+
+ pOutStartLoader->WriteDigits("OVERWRITE");
+
+ fSTART->ResetDigits();
+ pOutStartLoader->UnloadDigits();
+
}
+//------------------------------------------------------------------------
+Bool_t AliSTARTDigitizer::RegisterPhotoE(Int_t ipmt,Double_t energy)
+{
+
+
+ // Float_t hc=197.326960*1.e6; //mev*nm
+ Double_t hc=1.973*1.e-6; //gev*nm
+ Float_t lambda=hc/energy;
+ Float_t eff = ((TGraph*) fEffPMT.At(ipmt))->Eval(lambda);
+ Double_t p = gRandom->Rndm();
+
+ if (p > eff)
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//----------------------------------------------------------------------------
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff