X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALSDigitizer.cxx;h=a483a5bae41d6f643ee9eaf666f02e01cb607e0c;hb=452656d15558dace1996058a954ab3af884d3482;hp=2bf60b37c7568ec680191c852cb8aa97727a663e;hpb=ffa6d63bef8d3ba10bb8d254a05d7db3cea6f9a1;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALSDigitizer.cxx b/EMCAL/AliEMCALSDigitizer.cxx index 2bf60b37c75..a483a5bae41 100644 --- a/EMCAL/AliEMCALSDigitizer.cxx +++ b/EMCAL/AliEMCALSDigitizer.cxx @@ -16,10 +16,17 @@ /* $Id$ */ //_________________________________________________________________________ -// This is a TTask that makes SDigits out of Hits +// This is a Class that makes SDigits out of Hits // A Summable Digits is the sum of all hits originating // from one in one tower of the EMCAL // A threshold for assignment of the primary to SDigit is applied +// +// JLK 26-Jun-2008 Added explanation: +// SDigits need to hold the energy sum of the hits, but AliEMCALDigit +// can (should) only store amplitude. Therefore, the SDigit energy is +// "digitized" before being stored and must be "calibrated" back to an +// energy before SDigits are summed to form true Digits +// // SDigits are written to TreeS, branch "EMCAL" // AliEMCALSDigitizer with all current parameters is written // to TreeS branch "AliEMCALSDigitizer". @@ -29,13 +36,13 @@ // User case: // root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root") // Warning in : object already instantiated -// root [1] s->ExecuteTask() +// root [1] s->Digitize() // // Makes SDigitis for all events stored in galice.root // root [2] s->SetPedestalParameter(0.001) // // One can change parameters of digitization // root [3] s->SetSDigitsBranch("Redestal 0.001") // // and write them into the new branch -// root [4] s->ExeciteTask("deb all tim") +// root [4] s->Digitize("deb all tim") // // available parameters: // deb - print # of produced SDigitis // deb all - print # and list of produced SDigits @@ -45,363 +52,483 @@ // based on : AliPHOSSDigitzer ////////////////////////////////////////////////////////////////////////////// - // --- ROOT system --- -#include "TFile.h" -#include "TTask.h" -#include "TTree.h" -#include "TSystem.h" -#include "TROOT.h" -#include "TFolder.h" -#include "TBenchmark.h" +#include +#include +//#include +#include +#include + // --- Standard library --- -#include +#include "stdlib.h" // --- AliRoot header files --- -#include "AliRun.h" +#include "AliLog.h" +#include "AliRunLoader.h" +#include "AliStack.h" #include "AliEMCALDigit.h" +#include "AliEMCALLoader.h" #include "AliEMCALHit.h" #include "AliEMCALSDigitizer.h" #include "AliEMCALGeometry.h" -#include "AliEMCALv1.h" +#include "AliEMCALSimParam.h" ClassImp(AliEMCALSDigitizer) - //____________________________________________________________________________ - AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") +AliEMCALSDigitizer::AliEMCALSDigitizer() + : TNamed("",""), + fA(0.),fB(0.),fECPrimThreshold(0.), + fDefaultInit(kTRUE), + fEventFolderName(0), + fInit(0), + fSDigitsInRun(0), + fFirstEvent(0), + fLastEvent(0), + fSampling(0.), + fHits(0) { // ctor - fA = 0; - fB = 10000000. ; - fPrimThreshold = 0.01 ; - fNevents = 0 ; - fSDigits = 0 ; - fHits = 0 ; - fIsInitialized = kFALSE ; - + InitParameters(); } //____________________________________________________________________________ -AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","") +AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, + const char * eventFolderName) + : TNamed("EMCALSDigitizer", alirunFileName), + fA(0.),fB(0.),fECPrimThreshold(0.), + fDefaultInit(kFALSE), + fEventFolderName(eventFolderName), + fInit(0), + fSDigitsInRun(0), + fFirstEvent(0), + fLastEvent(0), + fSampling(0.), + fHits(0) { // ctor - fA = 0; - fB = 10000000.; - fPrimThreshold = 0.01 ; - fNevents = 0 ; - fSDigitsTitle = sDigitsTitle ; - fHeadersFile = headerFile ; - fSDigits = new TClonesArray("AliEMCALDigit",1000); - fHits = new TClonesArray("AliEMCALHit",1000); + Init(); + InitParameters() ; - TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ; - - //File was not opened yet - if(file == 0){ - if(fHeadersFile.Contains("rfio")) - file = TFile::Open(fHeadersFile,"update") ; - else - file = new TFile(fHeadersFile.Data(),"update") ; - gAlice = (AliRun *) file->Get("gAlice") ; - } - - //add Task to //root/Tasks folder - TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; - roottasks->Add(this) ; - - fIsInitialized = kTRUE ; } + //____________________________________________________________________________ -AliEMCALSDigitizer::~AliEMCALSDigitizer() +AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) + : TNamed(sd.GetName(),sd.GetTitle()), + fA(sd.fA), + fB(sd.fB), + fECPrimThreshold(sd.fECPrimThreshold), + fDefaultInit(sd.fDefaultInit), + fEventFolderName(sd.fEventFolderName), + fInit(sd.fInit), + fSDigitsInRun(sd.fSDigitsInRun), + fFirstEvent(sd.fFirstEvent), + fLastEvent(sd.fLastEvent), + fSampling(sd.fSampling), + fHits(0) { - // dtor - if(fSDigits) - delete fSDigits ; - if(fHits) - delete fHits ; + //cpy ctor } -//____________________________________________________________________________ -void AliEMCALSDigitizer::Init(){ - //Initialization can not be done in the default constructor - if(!fIsInitialized){ +//_____________________________________________________________________ +AliEMCALSDigitizer& AliEMCALSDigitizer::operator = (const AliEMCALSDigitizer& source) +{ // assignment operator; use copy ctor + if (&source == this) return *this; - if(fHeadersFile.IsNull()) - fHeadersFile="galice.root" ; + new (this) AliEMCALSDigitizer(source); + return *this; +} - TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data() ) ; - - //if file was not opened yet, read gAlice - if(file == 0){ - file = new TFile(fHeadersFile.Data(),"update") ; - gAlice = (AliRun *) file->Get("gAlice") ; - } - - fHits = new TClonesArray("AliEMCALHit",1000); - fSDigits = new TClonesArray("AliEMCALDigit",1000); - - // add Task to //root/Tasks folder - TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; - roottasks->Add(this) ; - - fIsInitialized = kTRUE ; +//____________________________________________________________________________ +AliEMCALSDigitizer::~AliEMCALSDigitizer() { + //dtor + + if(fHits){ + fHits->Clear(); + delete fHits; } } -//____________________________________________________________________________ -void AliEMCALSDigitizer::Exec(Option_t *option) { - //Collects all hits in the same active volume into digit - - if(!fIsInitialized) - Init() ; - if(strstr(option,"tim")) - gBenchmark->Start("EMCALSDigitizer"); - - fNevents = (Int_t) gAlice->TreeE()->GetEntries() ; - - Int_t ievent ; - for(ievent = 0; ievent < fNevents; ievent++){ - gAlice->GetEvent(ievent) ; - gAlice->SetEvent(ievent) ; - - if(gAlice->TreeH()==0){ - cout << "AliEMCALSDigitizer: There is no Hit Tree" << endl; - return ; - } - - //set address of the hits - TBranch * branch = gAlice->TreeH()->GetBranch("EMCAL"); - if (branch) - branch->SetAddress(&fHits); - else{ - cout << "ERROR in AliEMCALSDigitizer: "<< endl ; - cout << " no branch EMCAL in TreeH"<< endl ; - cout << " do nothing " << endl ; - return ; - } - - fSDigits->Clear(); - Int_t nSdigits = 0 ; - - - //Now made SDigits from hits, for EMCAL it is the same, so just copy - Int_t itrack ; - for (itrack=0; itrack < gAlice->GetNtrack(); itrack++){ - - //=========== Get the EMCAL branch from Hits Tree for the Primary track itrack - branch->GetEntry(itrack,0); - AliEMCALv0 * EMCAL = (AliEMCALv0 *) gAlice->GetDetector("EMCAL") ; - AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() ); - - Int_t i; - for ( i = 0 ; i < fHits->GetEntries() ; i++ ) { - AliEMCALHit * hit = (AliEMCALHit*)fHits->At(i) ; - AliEMCALDigit * curSDigit = 0 ; - AliEMCALDigit * sdigit ; - Bool_t newsdigit = kTRUE; -// Assign primary number only if contribution is significant - if( hit->GetEnergy() > fPrimThreshold) - curSDigit = new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1 ) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; - else - curSDigit = new AliEMCALDigit( -1 , -1 ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1 ) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; - cout << "Hit ID = " <GetId() << endl ; - cout << "ID for detector = " << curSDigit->GetId() << endl ; - cout << hit->GetEnergy() << " - hit energy - Digit Energy - " << curSDigit->GetAmp() << endl; - for(Int_t check= 0; check < nSdigits ; check++) { - sdigit = (AliEMCALDigit *)fSDigits->At(check); - if( sdigit->GetId() == curSDigit->GetId()) - { cout << "SDigit - Get Amp " << sdigit->GetAmp() << endl ; - *sdigit = *sdigit + *curSDigit ; - newsdigit = kFALSE; - cout << " and after addition " << sdigit->GetAmp() << endl ; - } - } - if (newsdigit) - { new((*fSDigits)[nSdigits]) AliEMCALDigit(*curSDigit); - nSdigits++ ; - cout << "Detector nsdigits = " << nSdigits << endl ; } - newsdigit = kTRUE; - +//____________________________________________________________________________ +void AliEMCALSDigitizer::Init(){ + // Initialization: open root-file, allocate arrays for hits and sdigits + // + // Initialization can not be done in the default constructor + //============================================================= YS + // The initialisation is now done by the getter - if( hit->GetEnergy() > fPrimThreshold) - curSDigit = new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), ((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ() + 1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; - else - curSDigit = new AliEMCALDigit( -1 , -1 ,((geom->GetNZ() * geom->GetNPhi()) + ((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; - - if((hit->GetId()/geom->GetNPhi()) < (2*geom->GetNZ())) - { cout << "ID for Preshower = " << curSDigit->GetId() << endl ; - for(Int_t check= 0; check < nSdigits; check++) { - sdigit = (AliEMCALDigit *)fSDigits->At(check); - if( sdigit->GetId() == curSDigit->GetId()) - { - *sdigit = *sdigit + *curSDigit ; - newsdigit = kFALSE ; - } - } - if (newsdigit) - { new((*fSDigits)[nSdigits]) AliEMCALDigit(*curSDigit); - nSdigits++ ; - cout << "Preshower nsdigits = " << nSdigits << endl ;} - newsdigit=kTRUE; - } - } - } // loop over tracks - - fSDigits->Sort() ; - - nSdigits = fSDigits->GetEntriesFast() ; - fSDigits->Expand(nSdigits) ; - - Int_t i ; - for (i = 0 ; i < nSdigits ; i++) { - AliEMCALDigit * digit = (AliEMCALDigit *) fSDigits->At(i) ; - digit->SetIndexInList(i) ; - } + fInit = kTRUE; + + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::Instance()->GetDetectorLoader("EMCAL")); - if(gAlice->TreeS() == 0) - gAlice->MakeTree("S") ; - - //check, if this branch already exits? - TBranch * sdigitsBranch = 0; - TBranch * sdigitizerBranch = 0; - - TObjArray * branches = gAlice->TreeS()->GetListOfBranches() ; - Int_t ibranch; - Bool_t emcalNotFound = kTRUE ; - Bool_t sdigitizerNotFound = kTRUE ; - - for(ibranch = 0;ibranch GetEntries();ibranch++){ - - if(emcalNotFound){ - sdigitsBranch=(TBranch *) branches->At(ibranch) ; - if( (strcmp("EMCAL",sdigitsBranch->GetName())==0 ) && - (fSDigitsTitle.CompareTo(sdigitsBranch->GetTitle()) == 0) ) - emcalNotFound = kFALSE ; - } - if(sdigitizerNotFound){ - sdigitizerBranch = (TBranch *) branches->At(ibranch) ; - if( (strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0)&& - (fSDigitsTitle.CompareTo(sdigitizerBranch->GetTitle()) == 0) ) - sdigitizerNotFound = kFALSE ; - } - } + if ( emcalLoader == 0 ) { + Fatal("Init", "Could not obtain the AliEMCALLoader"); + return ; + } + +} - if(!(sdigitizerNotFound && emcalNotFound)){ - cout << "AliEMCALSdigitizer error:" << endl ; - cout << "Can not overwrite existing branches: do not write" << endl ; - return ; - } - - //Make (if necessary) branches - char * file =0; - if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name - file = new char[strlen(gAlice->GetBaseFile())+20] ; - sprintf(file,"%s/EMCAL.SDigits.root",gAlice->GetBaseFile()) ; - } - - TDirectory *cwd = gDirectory; - - //First list of sdigits - Int_t bufferSize = 32000 ; - sdigitsBranch = gAlice->TreeS()->Branch("EMCAL",&fSDigits,bufferSize); - sdigitsBranch->SetTitle(fSDigitsTitle.Data()); - if (file) { - sdigitsBranch->SetFile(file); - TIter next( sdigitsBranch->GetListOfBranches()); - TBranch * subbr; - while ((subbr=(TBranch*)next())) { - subbr->SetFile(file); - } - cwd->cd(); - } - - //second - SDigitizer - Int_t splitlevel = 0 ; - AliEMCALSDigitizer * sd = this ; - sdigitizerBranch = gAlice->TreeS()->Branch("AliEMCALSDigitizer","AliEMCALSDigitizer", - &sd,bufferSize,splitlevel); - sdigitizerBranch->SetTitle(fSDigitsTitle.Data()); - if (file) { - sdigitizerBranch->SetFile(file); - TIter next( sdigitizerBranch->GetListOfBranches()); - TBranch * subbr ; - while ((subbr=(TBranch*)next())) { - subbr->SetFile(file); - } - cwd->cd(); - delete file; - } +//____________________________________________________________________________ +void AliEMCALSDigitizer::InitParameters() +{ + //initialize parameters for sdigitization - sdigitsBranch->Fill() ; - sdigitizerBranch->Fill() ; - gAlice->TreeS()->Write(0,TObject::kOverwrite) ; - - if(strstr(option,"deb")) - PrintSDigits(option) ; - + const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); + if (geom->GetSampling() == 0.) { + Fatal("InitParameters", "Sampling factor not set !") ; } - - if(strstr(option,"tim")){ - gBenchmark->Stop("EMCALSDigitizer"); - cout << "AliEMCALSDigitizer:" << endl ; - cout << " took " << gBenchmark->GetCpuTime("EMCALSDigitizer") << " seconds for SDigitizing " - << gBenchmark->GetCpuTime("EMCALSDigitizer")/fNevents << " seconds per event " << endl ; - cout << endl ; + + // Get the parameters from the OCDB via the loader + AliRunLoader *rl = AliRunLoader::Instance(); + AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); + AliEMCALSimParam * simParam = 0x0; + if(emcalLoader) simParam = emcalLoader->SimulationParameters(); + + if(!simParam){ + simParam = AliEMCALSimParam::GetInstance(); + AliWarning("Simulation Parameters not available in OCDB?"); } + + // + //JLK 26-Jun-2008 THIS SHOULD HAVE BEEN EXPLAINED AGES AGO: + // + //In order to be able to store SDigit Energy info into + //AliEMCALDigit, we need to convert it temporarily to an ADC amplitude + //and later when summing SDigits to form digits, convert it back to + //energy. These fA and fB parameters accomplish this through the + //Digitize() and Calibrate() methods + // + // Initializes parameters + fA = simParam->GetA(); //0; + fB = simParam->GetB(); //1.e+6; // Changed 24 Apr 2007. Dynamic range now 2 TeV + fSampling = geom->GetSampling(); + + // threshold for deposit energy of hit + fECPrimThreshold = simParam->GetECPrimaryThreshold();//0.05;// GeV // 22-may-07 was 0// 24-nov-04 - was 1.e-6; + AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName())); + AliDebug(2,Form(" fInit %i\n", int(fInit))); + AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent)); + AliDebug(2,Form(" fLastEvent %i\n", fLastEvent)); + AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data())); + AliDebug(2,Form(" with digitization parameters A = %f\n", fA)); + AliDebug(2,Form(" B = %f\n", fB)); + AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold)); + AliDebug(2,Form(" Sampling = %f\n", fSampling)); + AliDebug(2,Form("---------------------------------------------------\n")); + +} + +//____________________________________________________________________________ +void AliEMCALSDigitizer::Digitize(Option_t *option) +{ + // Collects all hit of the same tower into digits + TString o(option); o.ToUpper(); + if (strstr(option, "print") ) { + + AliDebug(2,Form("Print: \n------------------- %s -------------\n",GetName())); + AliDebug(2,Form(" fInit %i\n", int(fInit))); + AliDebug(2,Form(" fFirstEvent %i\n", fFirstEvent)); + AliDebug(2,Form(" fLastEvent %i\n", fLastEvent)); + AliDebug(2,Form(" Writing SDigits to branch with title %s\n", fEventFolderName.Data())); + AliDebug(2,Form(" with digitization parameters A = %f\n", fA)); + AliDebug(2,Form(" B = %f\n", fB)); + AliDebug(2,Form(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold)); + AliDebug(2,Form(" Sampling = %f\n", fSampling)); + AliDebug(2,Form("---------------------------------------------------\n")); + + return ; + } + + + if(strstr(option,"tim")) + gBenchmark->Start("EMCALSDigitizer"); + + AliRunLoader *rl = AliRunLoader::Instance(); + AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); + + if (!fInit) { // to prevent overwrite existing file + AliError( Form("Give a version name different from %s", fEventFolderName.Data()) ) ; + return ; + } + + if (fLastEvent == -1) + fLastEvent = rl->GetNumberOfEvents() - 1 ; + else { + fLastEvent = TMath::Min(fLastEvent, rl->GetNumberOfEvents()-1); + } + Int_t nEvents = fLastEvent - fFirstEvent + 1; + + Float_t energy=0.; // de * fSampling - 23-nov-04 + rl->LoadKinematics(); + rl->LoadHits("EMCAL"); + Int_t ievent; + for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { + rl->GetEvent(ievent); + TTree * treeS = emcalLoader->TreeS(); + if ( !treeS ) { + emcalLoader->MakeSDigitsContainer(); + treeS = emcalLoader->TreeS(); + } + + TClonesArray * sdigits = emcalLoader->SDigits() ; + sdigits->Clear("C"); + + Int_t nSdigits = 0 ; + Int_t iHit, iTrack, iSDigit; + + AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); + + TTree *treeH = emcalLoader->TreeH(); + if (treeH) { + Int_t nTrack = treeH->GetEntries(); // TreeH has array of hits for every primary + TBranch * branchH = treeH->GetBranch("EMCAL"); + //if(fHits)fHits->Clear(); + branchH->SetAddress(&fHits); + for (iTrack = 0; iTrack < nTrack; iTrack++) { + branchH->GetEntry(iTrack); + + if(fHits){ + + Int_t nHit = fHits->GetEntriesFast(); + for(iHit = 0; iHit< nHit;iHit++){ + + AliEMCALHit * hit = dynamic_cast(fHits->At(iHit)) ; + AliEMCALDigit * curSDigit = 0 ; + AliEMCALDigit * sdigit = 0 ; + Bool_t newsdigit = kTRUE; + + // hit->GetId() - Absolute Id number EMCAL segment + if(hit){ + if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId()) + energy = hit->GetEnergy() * fSampling; // 23-nov-04 + if(energy > fECPrimThreshold ) + // Assign primary number only if deposited energy is significant + curSDigit = new AliEMCALDigit(hit->GetPrimary(), + hit->GetIparent(), hit->GetId(), + Digitize(energy), hit->GetTime(),kFALSE, + -1, 0,0,energy ) ; + else + curSDigit = new AliEMCALDigit(-1, + -1, + hit->GetId(), + Digitize(energy), hit->GetTime(),kFALSE, + -1, 0,0,energy ) ; + } else { + Warning("Digitize"," abs id %i is bad \n", hit->GetId()); + newsdigit = kFALSE; + curSDigit = 0; + } + + if(curSDigit != 0){ + for(Int_t check= 0; check < nSdigits ; check++) { + sdigit = dynamic_cast(sdigits->At(check)) ; + if(sdigit){ + if( sdigit->GetId() == curSDigit->GetId()) { // Are we in the same ECAL tower ? + *sdigit = *sdigit + *curSDigit; + newsdigit = kFALSE; + } + }// sdigit exists + else { + AliWarning("Sdigit do not exist"); + newsdigit = kFALSE; + }// sdigit does not exist + }//sdigit loop + }// currsdigit exists + + if (newsdigit) { + new((*sdigits)[nSdigits]) AliEMCALDigit(*curSDigit); + nSdigits++ ; + } + delete curSDigit ; + + }// hit exists + else AliFatal("Hit is NULL!"); + + } // loop over all hits (hit = deposited energy/entering particle) + + }//fHits is not NULL + else AliFatal("fHit is NULL!"); + + sdigits->Sort() ; + + nSdigits = sdigits->GetEntriesFast() ; + fSDigitsInRun += nSdigits ; + + for (iSDigit = 0 ; iSDigit < sdigits->GetEntriesFast() ; iSDigit++) { + AliEMCALDigit * sdigit = dynamic_cast(sdigits->At(iSDigit)) ; + if(sdigit)sdigit->SetIndexInList(iSDigit) ; + else AliFatal("sdigit is NULL!"); + } + if(fHits)fHits->Clear(); + }//track loop + }// tree exists + + // Now write SDigits + + Int_t bufferSize = 32000 ; + TBranch * sdigitsBranch = treeS->GetBranch("EMCAL"); + if (sdigitsBranch) + sdigitsBranch->SetAddress(&sdigits); + else + treeS->Branch("EMCAL",&sdigits,bufferSize); + + treeS->Fill(); + + emcalLoader->WriteSDigits("OVERWRITE"); + + //NEXT - SDigitizer + //emcalLoader->WriteSDigitizer("OVERWRITE"); // why in event cycle ? + + if(strstr(option,"deb")) + PrintSDigits(option) ; + } + + Unload(); + + if(strstr(option,"tim")){ + gBenchmark->Stop("EMCALSDigitizer"); + printf("\n Digitize: took %f seconds for SDigitizing %f seconds per event\n", + gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; + } } + //__________________________________________________________________ -void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){ - //Seting title to branch SDigits - if(!fSDigitsTitle.IsNull()) - cout << "AliEMCALSdigitizer: changing SDigits file from " <2147483647.0) { + //PH 2147483647 is the max. integer + //PH This apparently is a problem which needs investigation + AliWarning(Form("Too big or too small energy %f",aSignal)); + aSignal = TMath::Sign((Double_t)2147483647,aSignal); + } + + return (Float_t ) aSignal; } + //__________________________________________________________________ -void AliEMCALSDigitizer::Print(Option_t* option)const{ - cout << "------------------- "<< GetName() << " -------------" << endl ; - cout << " Writing SDigitis to branch with title " << fSDigitsTitle.Data() << endl ; - cout << " with digitization parameters A = " << fA << endl ; - cout << " B = " << fB << endl ; - cout << " Threshold for Primary assignment= " << fPrimThreshold << endl ; - cout << "---------------------------------------------------"<GetEntriesFast() << endl ; - cout << endl ; - - if(strstr(option,"all")){// print all digits + AliEMCALLoader *rl = dynamic_cast(AliRunLoader::Instance()->GetDetectorLoader("EMCAL")); + if(rl){ + const TClonesArray * sdigits = rl->SDigits() ; - //loop over digits - AliEMCALDigit * digit; - cout << "SDigit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl; - Int_t index ; - for (index = 0 ; index < fSDigits->GetEntries() ; index++) { - digit = (AliEMCALDigit * ) fSDigits->At(index) ; - cout << setw(8) << digit->GetId() << " " << setw(3) << digit->GetAmp() << " " - << setw(6) << digit->GetIndexInList() << " " - << setw(5) << digit->GetNprimary() <<" "; + printf("\n") ; + printf("event %i", rl->GetRunLoader()->GetEventNumber()); + printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast()); + if(strstr(option,"all")||strstr(option,"EMC")){ - Int_t iprimary; - for (iprimary=0; iprimaryGetNprimary(); iprimary++) - cout << setw(5) << digit->GetPrimary(iprimary+1) << " "; - cout << endl; - } - + //loop over digits + AliEMCALDigit * digit; + printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ; + Int_t index = 0; + Float_t isum = 0.; + const Int_t bufferSize = 8192; + char * tempo = new char[bufferSize]; + for (index = 0 ; index < sdigits->GetEntries() ; index++) { + digit = dynamic_cast( sdigits->At(index) ) ; + if(digit){ + snprintf(tempo, bufferSize,"\n%6d %8f %6.5e %4d %2d :", + digit->GetId(), digit->GetAmplitude(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ; + printf("%s",tempo); + isum += digit->GetAmplitude(); + + Int_t iprimary; + for (iprimary=0; iprimaryGetNprimary(); iprimary++) { + snprintf(tempo,bufferSize, "%d ",digit->GetPrimary(iprimary+1) ) ; + printf("%s",tempo); + } + } //sdigit exists + else AliFatal("SDigit is NULL!"); + }//loop + delete [] tempo ; + printf("\n** Sum %2.3f : %10.3f GeV/c **\n ", isum, Calibrate(isum)); + } else printf("\n"); + } + else AliFatal("EMCALLoader is NULL!"); +} + +//____________________________________________________________________________ +void AliEMCALSDigitizer::Unload() const +{ + // Unload Hits and SDigits from the folder + AliEMCALLoader *rl = dynamic_cast(AliRunLoader::Instance()->GetDetectorLoader("EMCAL")); + if(rl){ + rl->UnloadHits() ; + rl->UnloadSDigits() ; } + else AliFatal("EMCALLoader is NULL!"); +} + +//____________________________________________________________________________ +void AliEMCALSDigitizer::Browse(TBrowser* b) +{ + TNamed::Browse(b); }