From 62741b41ca9a389bc9e3f6f650f6bf8df5bfb125 Mon Sep 17 00:00:00 2001 From: schutz Date: Wed, 8 Jan 2003 17:24:55 +0000 Subject: [PATCH] added the HCAL section and removed obsolete method --- EMCAL/AliEMCALSDigitizer.cxx | 308 ++++++++++++++++------------------- EMCAL/AliEMCALSDigitizer.h | 19 +-- 2 files changed, 147 insertions(+), 180 deletions(-) diff --git a/EMCAL/AliEMCALSDigitizer.cxx b/EMCAL/AliEMCALSDigitizer.cxx index d8df1253cd2..1639b06bc22 100644 --- a/EMCAL/AliEMCALSDigitizer.cxx +++ b/EMCAL/AliEMCALSDigitizer.cxx @@ -76,6 +76,8 @@ ClassImp(AliEMCALSDigitizer) { // ctor InitParameters() ; + fSplitFile = 0 ; + fToSplit = kFALSE ; fDefaultInit = kTRUE ; } @@ -84,9 +86,9 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigi TTask(sDigitsTitle, headerFile) { // ctor - InitParameters() ; fToSplit = toSplit ; Init(); + InitParameters() ; fDefaultInit = kFALSE ; } @@ -118,6 +120,7 @@ void AliEMCALSDigitizer::Init(){ gime->PostSDigits( GetName(), GetTitle() ) ; fSplitFile = 0 ; + if(fToSplit){ // construct the name of the file as /path/EMCAL.SDigits.root // First - extract full path if necessary @@ -147,6 +150,7 @@ void AliEMCALSDigitizer::Init(){ sdname.Append(GetTitle() ) ; SetName(sdname) ; gime->PostSDigitizer(this) ; + } //____________________________________________________________________________ @@ -154,17 +158,27 @@ void AliEMCALSDigitizer::InitParameters() { fA = 0; fB = 10000000.; - fTowerPrimThreshold = 0.01 ; - fPreShowerPrimThreshold = 0.0001 ; - fPhotonElectronFactor = 5000. ; // photoelectrons per GeV - fSplitFile = 0 ; - fToSplit = kFALSE ; + + AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; + const AliEMCALGeometry * geom = gime->EMCALGeometry() ; + if (geom->GetSampling() == 0.) { + Error("InitParameters", "Sampling factor not set !") ; + abort() ; + } + else + Info("InitParameters", "Sampling factor set to %f\n", geom->GetSampling()) ; + + // this threshold corresponds approximately to 100 MeV + fECPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNECLayers() ; + fPREPrimThreshold = 100E-3 / ( geom->GetSampling() * ( geom->GetNPRLayers() + geom->GetNECLayers()) ) * geom->GetNPRLayers() ; + fHCPrimThreshold = fECPrimThreshold/5. ; // 5 is totally arbitrary + } //____________________________________________________________________________ void AliEMCALSDigitizer::Exec(Option_t *option) { - // Collects all hits in the same active volume into digit + // Collects all hits in the section (PRE/ECAL/HCAL) of the same tower into digit if( strcmp(GetName(), "") == 0 ) Init() ; @@ -186,163 +200,137 @@ void AliEMCALSDigitizer::Exec(Option_t *option) sdname.Remove(sdname.Index(GetTitle())-1) ; Int_t nevents = gime->MaxEvent() ; - Int_t ievent ; - - for(ievent = 0; ievent < nevents; ievent++){ - gime->Event(ievent,"H") ; - - const TClonesArray * hits = gime->Hits() ; - + Int_t ievent ; + for(ievent = 0; ievent < nevents; ievent++){ + gime->Event(ievent,"H") ; + const TClonesArray * hits = gime->Hits() ; TClonesArray * sdigits = gime->SDigits(sdname.Data()) ; sdigits->Clear(); Int_t nSdigits = 0 ; - - //Collects all hits in the same active volume into digit //Now make SDigits from hits, for EMCAL it is the same, so just copy Int_t nPrim = static_cast((gAlice->TreeH())->GetEntries()) ; // Attention nPrim is the number of primaries tracked by Geant // and this number could be different to the number of Primaries in TreeK; Int_t iprim ; - for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { - //=========== Get the EMCAL branch from Hits Tree for the Primary iprim - gime->Track(iprim) ; - Int_t i; - for ( i = 0 ; i < hits->GetEntries() ; i++ ) { - AliEMCALHit * hit = dynamic_cast(hits->At(i)) ; - AliEMCALDigit * curSDigit = 0 ; - AliEMCALDigit * sdigit = 0 ; - Bool_t newsdigit = kTRUE; - // Assign primary number only if deposited energy is significant - - if( (!hit->IsInPreShower() && hit->GetEnergy() > fTowerPrimThreshold) || - (hit->IsInPreShower() && hit->GetEnergy() > fPreShowerPrimThreshold)) + for ( iprim = 0 ; iprim < nPrim ; iprim++ ) { + //=========== Get the EMCAL branch from Hits Tree for the Primary iprim + gime->Track(iprim) ; + Int_t i; + for ( i = 0 ; i < hits->GetEntries() ; i++ ) { + AliEMCALHit * hit = dynamic_cast(hits->At(i)) ; + AliEMCALDigit * curSDigit = 0 ; + AliEMCALDigit * sdigit = 0 ; + Bool_t newsdigit = kTRUE; + + // Assign primary number only if deposited energy is significant + + const AliEMCALGeometry * geom = gime->EMCALGeometry() ; + + if (gDebug) + Info("Exec", "id = %d energy = %f thresholdPRE = %f thresholdEC = %f thresholdHC = %f \n", + hit->GetId(), hit->GetEnergy(), fPREPrimThreshold, fECPrimThreshold, fHCPrimThreshold) ; + + if( geom->IsInPRE(hit->GetId()) ) + if( hit->GetEnergy() > fPREPrimThreshold ) curSDigit = new AliEMCALDigit( hit->GetPrimary(), - hit->GetIparent(),Layer2TowerID(hit->GetId(),hit->IsInPreShower()), + hit->GetIparent(), hit->GetId(), Digitize(hit->GetEnergy()), hit->GetTime() ) ; - else + else curSDigit = new AliEMCALDigit( -1 , -1 , - Layer2TowerID(hit->GetId(),hit->IsInPreShower()), - Digitize(hit->GetEnergy()), hit->GetTime() ) ; - Int_t check = 0 ; - for(check= 0; check < nSdigits ; check++) { - sdigit = dynamic_cast(sdigits->At(check)) ; - if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same tower or the same preshower ? - *sdigit = *sdigit + *curSDigit; - newsdigit = kFALSE; - } - } - if (newsdigit) { - new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit); - nSdigits++ ; - } - delete curSDigit ; - } // loop over all hits (hit = deposited energy/layer/entering particle) - } // loop over iprim - - sdigits->Sort() ; - - nSdigits = sdigits->GetEntriesFast() ; - fSDigitsInRun += nSdigits ; - sdigits->Expand(nSdigits) ; - - const AliEMCALGeometry * geom = gime->EMCALGeometry() ; - - if (nSdigits != 0 ) { - Int_t lastPreShowerIndex = nSdigits - 1 ; - - - if (!(dynamic_cast(sdigits->At(lastPreShowerIndex))->IsInPreShower())) - - lastPreShowerIndex = -2; - - Int_t firstPreShowerIndex = 100000 ; - Int_t index ; - AliEMCALDigit * sdigit = 0 ; - - for ( index = 0; index < nSdigits ; index++) { - sdigit = dynamic_cast(sdigits->At(index) ) ; - if (sdigit->IsInPreShower() ){ - firstPreShowerIndex = index ; - break ; - } - } - - AliEMCALDigit * preshower ; - AliEMCALDigit * tower ; - Int_t lastIndex = lastPreShowerIndex +1 ; + hit->GetId(), + Digitize(hit->GetEnergy()), hit->GetTime() ) ; + else if( geom->IsInECAL(hit->GetId()) ) + if( hit->GetEnergy() > fECPrimThreshold ) + curSDigit = new AliEMCALDigit( hit->GetPrimary(), + hit->GetIparent(), hit->GetId(), + Digitize(hit->GetEnergy()), hit->GetTime() ) ; + else + curSDigit = new AliEMCALDigit( -1 , + -1 , + hit->GetId(), + Digitize(hit->GetEnergy()), hit->GetTime() ) ; + else if( geom->IsInHCAL(hit->GetId()) ) + if( hit->GetEnergy() > fHCPrimThreshold ) + + curSDigit = new AliEMCALDigit( hit->GetPrimary(), + hit->GetIparent(), hit->GetId(), + Digitize(hit->GetEnergy()), hit->GetTime() ) ; + else + curSDigit = new AliEMCALDigit( -1 , + -1 , + hit->GetId(), + Digitize(hit->GetEnergy()), hit->GetTime() ) ; - for (index = firstPreShowerIndex ; index <= lastPreShowerIndex; index++) { - preshower = dynamic_cast(sdigits->At(index) ); - Bool_t towerFound = kFALSE ; - Int_t jndex ; - for (jndex = 0; jndex < firstPreShowerIndex; jndex++) { - tower = dynamic_cast(sdigits->At(jndex) ); - if ( (preshower->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) == tower->GetId() ) { - Float_t towerEnergy = static_cast(tower->GetAmp()) ; - Float_t preshoEnergy = static_cast(preshower->GetAmp()) ; - towerEnergy +=preshoEnergy ; - *tower = *tower + *preshower ; // and add preshower multiplied by layer ratio to tower - tower->SetAmp(static_cast(TMath::Ceil(towerEnergy))) ; - towerFound = kTRUE ; - } + Int_t check = 0 ; + for(check= 0; check < nSdigits ; check++) { + sdigit = dynamic_cast(sdigits->At(check)) ; + if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL/HCAL/preshower tower ? + *sdigit = *sdigit + *curSDigit; + newsdigit = kFALSE; } - if ( !towerFound ) { - new((*sdigits)[lastIndex]) AliEMCALDigit(*preshower); - AliEMCALDigit * temp = dynamic_cast(sdigits->At(lastIndex)) ; - temp->SetId(temp->GetId() - (geom->GetNZ() * geom->GetNPhi()) ) ; - lastIndex++ ; - } - } - sdigits->Sort() ; - Int_t NPrimarymax = -1 ; - Int_t i ; - for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { - sdigit = dynamic_cast(sdigits->At(i)) ; - sdigit->SetIndexInList(i) ; } - - for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { - if (((dynamic_cast(sdigits->At(i)))->GetNprimary()) > NPrimarymax) - NPrimarymax = ((dynamic_cast(sdigits->At(i)))->GetNprimary()) ; + if (newsdigit) { + new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit); + nSdigits++ ; } - } - - //Now write SDigits - - if(gAlice->TreeS() == 0 || (fSplitFile)) //<--- To be checked: we should not create TreeS if it is already here - gAlice->MakeTree("S",fSplitFile); - - if(fSplitFile) - fSplitFile->cd() ; - - //First list of sdigits - Int_t bufferSize = 32000 ; - TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize); - sdigitsBranch->SetTitle(sdname); - - //NEXT - SDigitizer - Int_t splitlevel = 0 ; - AliEMCALSDigitizer * sd = this ; - TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer", - &sd,bufferSize,splitlevel); - sdigitizerBranch->SetTitle(sdname); - - sdigitsBranch->Fill() ; - sdigitizerBranch->Fill() ; - gAlice->TreeS()->AutoSave() ; - - if(strstr(option,"deb")) - PrintSDigits(option) ; - + delete curSDigit ; + } // loop over all hits (hit = deposited energy/layer/entering particle) + } // loop over iprim + + sdigits->Sort() ; + + nSdigits = sdigits->GetEntriesFast() ; + fSDigitsInRun += nSdigits ; + sdigits->Expand(nSdigits) ; + + Int_t NPrimarymax = -1 ; + Int_t i ; + for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { + AliEMCALDigit * sdigit = dynamic_cast(sdigits->At(i)) ; + sdigit->SetIndexInList(i) ; + } + + for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { + if (((dynamic_cast(sdigits->At(i)))->GetNprimary()) > NPrimarymax) + NPrimarymax = ((dynamic_cast(sdigits->At(i)))->GetNprimary()) ; + } + + // Now write SDigits + + if(gAlice->TreeS() == 0 || (fSplitFile)) //<--- To be checked: we should not create TreeS if it is already here + gAlice->MakeTree("S",fSplitFile); + + if(fSplitFile) + fSplitFile->cd() ; + + //First list of sdigits + Int_t bufferSize = 32000 ; + TBranch * sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&sdigits,bufferSize); + sdigitsBranch->SetTitle(sdname); + + //NEXT - SDigitizer + Int_t splitlevel = 0 ; + AliEMCALSDigitizer * sd = this ; + TBranch * sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer", + &sd,bufferSize,splitlevel); + sdigitizerBranch->SetTitle(sdname); + + sdigitsBranch->Fill() ; + sdigitizerBranch->Fill() ; + + gAlice->TreeS()->AutoSave() ; + + if(strstr(option,"deb")) + PrintSDigits(option) ; } - + if(strstr(option,"tim")){ gBenchmark->Stop("EMCALSDigitizer"); Info("Exec", "took %f seconds for SDigitizing %f seconds per event", - gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ; - } + gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer") ) ; + } } //__________________________________________________________________ @@ -389,10 +377,12 @@ void AliEMCALSDigitizer::Print(Option_t* option)const message += fA ; message += "\n B = " ; message += fB ; - message += "\n Threshold for Primary assignment in Tower = " ; - message += fTowerPrimThreshold ; message += "\n Threshold for Primary assignment in PreShower = " ; - message += fPreShowerPrimThreshold ; + message += fPREPrimThreshold ; + message += "\n Threshold for Primary assignment in EC section= " ; + message += fECPrimThreshold ; + message += "\n Threshold for Primary assignment in HC section= " ; + message += fHCPrimThreshold ; message += "\n---------------------------------------------------" ; Info("Print", message.Data() ) ; @@ -405,8 +395,9 @@ Bool_t AliEMCALSDigitizer::operator==( AliEMCALSDigitizer const &sd )const // SDititizers are equal if their pedestal, slope and threshold are equal if( (fA==sd.fA)&&(fB==sd.fB)&& - (fTowerPrimThreshold==sd.fTowerPrimThreshold) && - (fPreShowerPrimThreshold==sd.fPreShowerPrimThreshold)) + (fECPrimThreshold==sd.fECPrimThreshold) && + (fHCPrimThreshold==sd.fHCPrimThreshold) && + (fPREPrimThreshold==sd.fPREPrimThreshold)) return kTRUE ; else return kFALSE ; @@ -450,27 +441,6 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option){ Info("PrintSDigits", message.Data() ) ; } -//________________________________________________________________________ -const Int_t AliEMCALSDigitizer::Layer2TowerID(Int_t ihit, Bool_t preshower) -{ - // Method to Transform from Hit Id to Digit Id - // This function should be one to one - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - const AliEMCALGeometry * geom = gime->EMCALGeometry(); - Int_t ieta = ((ihit-1)/geom->GetNPhi())%geom->GetNZ(); // eta Tower Index - Int_t iphi = (ihit-1)%(geom->GetNPhi())+1; //phi Tower Index - Int_t it = -10; - Int_t ipre = 0; - - if (preshower)ipre = 1; - if (iphi > 0 && ieta >= 0){ - it = iphi+ieta*geom->GetNPhi() + ipre*geom->GetNPhi()*geom->GetNZ(); - return it; - }else{ - Error("Layer2TowerID", "there is an error: Eta number = %f Phi number = %f", ieta, iphi) ; - return it; - } // end if iphi>0 && ieta>0 -} //_______________________________________________________________________________________ // void AliEMCALSDigitizer::TestTowerID(void) // { diff --git a/EMCAL/AliEMCALSDigitizer.h b/EMCAL/AliEMCALSDigitizer.h index 4598365f439..24acce0f824 100644 --- a/EMCAL/AliEMCALSDigitizer.h +++ b/EMCAL/AliEMCALSDigitizer.h @@ -41,30 +41,27 @@ public: void SetSlopeParameter(Float_t B){fB = B ;} void UseHitsFrom(const char * filename) ; Bool_t operator == (const AliEMCALSDigitizer & sd) const ; - const Int_t Segment2TowerID(Int_t SegmentID){ - return Layer2TowerID(SegmentID,kFALSE); -} +// const Int_t Segment2TowerID(Int_t SegmentID){ +// return Layer2TowerID(SegmentID,0); +// } private: void Init() ; void InitParameters() ; void PrintSDigits(Option_t * option) ; - const Int_t Layer2TowerID(Int_t,Bool_t) ; private: Float_t fA ; // Pedestal parameter Float_t fB ; // Slope Digitizition parameters - Float_t fPhotonElectronFactor ; // number of photon electrons per GeV - // should be calculated independently for each layer as : - // LightYield*LightCollectionEfficiency*LightAttenuation*APDPhotoElectronEfficiency*APDGain - Float_t fTowerPrimThreshold ; // To store primary in Tower if Elos > threshold - Float_t fPreShowerPrimThreshold ;// To store primary if Pre Shower Elos > threshold - Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) + Float_t fPREPrimThreshold ; // To store primary if Pre Shower Elos > threshold + Float_t fECPrimThreshold ; // To store primary if EC Shower Elos > threshold + Float_t fHCPrimThreshold ; // To store primary if HC Shower Elos > threshold + Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) Int_t fSDigitsInRun ; //! Total number of sdigits in one run TFile * fSplitFile ; //! file in which SDigits will eventually be stored Bool_t fToSplit ; //! Says that sigits should be written into splip file - ClassDef(AliEMCALSDigitizer,3) // description + ClassDef(AliEMCALSDigitizer,4) // description }; -- 2.39.3