From 61e0abb5c9bbd6e578ed6769602e4d9bb80ba13e Mon Sep 17 00:00:00 2001 From: hristov Date: Thu, 22 Nov 2001 17:41:00 +0000 Subject: [PATCH] New version of the EMCAL code: bug fixes, new digitisation classes (B.Nilsen) --- EMCAL/AliEMCAL.cxx | 17 +- EMCAL/AliEMCALDigit.cxx | 248 +++++++++++++ EMCAL/AliEMCALDigit.h | 71 ++++ EMCAL/AliEMCALDigitizer.cxx | 674 +++++++++++++++++++++++++++++++++++ EMCAL/AliEMCALDigitizer.h | 83 +++++ EMCAL/AliEMCALHit.cxx | 34 +- EMCAL/AliEMCALHit.h | 24 +- EMCAL/AliEMCALLinkDef.h | 12 - EMCAL/AliEMCALSDigitizer.cxx | 398 +++++++++++++++++++++ EMCAL/AliEMCALSDigitizer.h | 64 ++++ EMCAL/AliEMCALv0.cxx | 17 +- EMCAL/AliEMCALv1.cxx | 36 +- EMCAL/AliEMCALv1.h | 4 +- EMCAL/EMCALLinkDef.h | 3 + EMCAL/Makefile | 3 +- EMCAL/libEMCAL.pkg | 10 +- 16 files changed, 1643 insertions(+), 55 deletions(-) create mode 100644 EMCAL/AliEMCALDigit.cxx create mode 100644 EMCAL/AliEMCALDigit.h create mode 100644 EMCAL/AliEMCALDigitizer.cxx create mode 100644 EMCAL/AliEMCALDigitizer.h delete mode 100644 EMCAL/AliEMCALLinkDef.h create mode 100644 EMCAL/AliEMCALSDigitizer.cxx create mode 100644 EMCAL/AliEMCALSDigitizer.h diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index 1063b605d94..cda534cd451 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -111,7 +111,7 @@ void AliEMCAL::CreateMaterials(){ // The scintillator of the CPV made of Polystyrene scintillator -> idtmed[1601] AliMedium(2, "CPV scint. $", 2, 1, - isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, 0, 0) ; + isxfld, sxmgmx, 10.0, 0.001, 0.1, 0.001, 0.001, 0, 0) ; // Various Aluminium parts made of Al -> idtmed[1602] AliMedium(3, "Al parts $", 3, 0, @@ -123,17 +123,28 @@ void AliEMCAL::CreateMaterials(){ // --- Set decent energy thresholds for gamma and electron tracking // Tracking threshold for photons and electrons in Lead - gMC->Gstpar(idtmed[1600], "CUTGAM",0.5E-4) ; - gMC->Gstpar(idtmed[1600], "CUTELE",1.0E-4) ; + // gMC->Gstpar(idtmed[1600], "CUTGAM",0.5E-4) ; + // gMC->Gstpar(idtmed[1600], "CUTELE",1.0E-4) ; // --- Generate explicitly delta rays in Lead --- gMC->Gstpar(idtmed[1600], "LOSS",3.) ; gMC->Gstpar(idtmed[1600], "DRAY",1.) ; + gMC->Gstpar(idtmed[1600],"CUTGAM",0.00008) ; + gMC->Gstpar(idtmed[1600],"CUTELE",0.001) ; + gMC->Gstpar(idtmed[1600],"BCUTE",0.0001) ; + // --- and in aluminium parts --- gMC->Gstpar(idtmed[1602], "LOSS",3.) ; gMC->Gstpar(idtmed[1602], "DRAY",1.) ; + +// --- and finally in the scintillator --- + gMC->Gstpar(idtmed[1601],"CUTGAM",0.00008) ; + gMC->Gstpar(idtmed[1601],"CUTELE",0.001) ; + gMC->Gstpar(idtmed[1601],"BCUTE",0.0001) ; + + } //____________________________________________________________________________ void AliEMCAL::SetTreeAddress() diff --git a/EMCAL/AliEMCALDigit.cxx b/EMCAL/AliEMCALDigit.cxx new file mode 100644 index 00000000000..49ecc78f3cb --- /dev/null +++ b/EMCAL/AliEMCALDigit.cxx @@ -0,0 +1,248 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//_________________________________________________________________________ +// EMCAL digit: Id +// energy +// 3 identifiers for the primary particle(s) at the origine of the digit +// The digits are made in FinishEvent() by summing all the hits in a single EMCAL crystal or PPSD gas cell +// It would be nice to replace the 3 identifiers by an array, but, because digits are kept in a TClonesQArray, +// it is not possible to stream such an array... (beyond my understqnding!) +// +//*-- Author: Laurent Aphecetche & Yves Schutz (SUBATECH) + + +// --- ROOT system --- + +// --- Standard library --- + +#include + +// --- AliRoot header files --- + +#include "AliEMCALDigit.h" + + +ClassImp(AliEMCALDigit) + +//____________________________________________________________________________ + AliEMCALDigit::AliEMCALDigit() +{ + // default ctor + + fIndexInList = -1 ; + fNprimary = 0 ; + fNMaxPrimary = 5 ; + fNiparent = 0 ; + fNMaxiparent = fNMaxPrimary*10; +} + +//____________________________________________________________________________ +AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Int_t DigEnergy, Int_t index) +{ + // ctor with all data + + fNMaxPrimary = 5 ; + fNMaxiparent = fNMaxPrimary*10; + fAmp = DigEnergy ; + fId = id ; + fIndexInList = index ; + if( primary != -1){ + fNprimary = 1 ; + fPrimary[0] = primary ; + fNiparent = 1 ; + fIparent[0] = iparent ; + + } + else{ //If the contribution of this primary smaller than fDigitThreshold (AliEMCALv1) + fNprimary = 0 ; + fPrimary[0] = -1 ; + fNiparent = 0 ; + fIparent[0] = -1 ; + + } + Int_t i ; + for ( i = 1; i < fNMaxPrimary ; i++) + fPrimary[i] = -1 ; + + for ( Int_t j =1; j< fNMaxiparent ; j++) + fIparent[j] = -1 ; +} + +//____________________________________________________________________________ +AliEMCALDigit::AliEMCALDigit(const AliEMCALDigit & digit) +{ + // copy ctor + + + fNMaxPrimary = digit.fNMaxPrimary ; + fNMaxiparent = digit.fNMaxiparent ; + Int_t i ; + for ( i = 0; i < fNMaxPrimary ; i++) + fPrimary[i] = digit.fPrimary[i] ; + Int_t j ; + for (j = 0; j< fNMaxiparent ; j++) + fIparent[j] = digit.fIparent[j] ; + fAmp = digit.fAmp ; + fId = digit.fId; + fIndexInList = digit.fIndexInList ; + fNprimary = digit.fNprimary ; + fNiparent = digit.fNiparent ; +} + +//____________________________________________________________________________ +AliEMCALDigit::~AliEMCALDigit() +{ + // Delete array of primiries if any + +} + +//____________________________________________________________________________ +Int_t AliEMCALDigit::Compare(const TObject * obj) const +{ + // Compares two digits with respect to its Id + // to sort according increasing Id + + Int_t rv ; + + AliEMCALDigit * digit = (AliEMCALDigit *)obj ; + + Int_t iddiff = fId - digit->GetId() ; + + if ( iddiff > 0 ) + rv = 1 ; + else if ( iddiff < 0 ) + rv = -1 ; + else + rv = 0 ; + + return rv ; + +} + +//____________________________________________________________________________ +Int_t AliEMCALDigit::GetPrimary(Int_t index) const +{ + // retrieves the primary particle number given its index in the list + Int_t rv = -1 ; + if ( index <= fNprimary ){ + rv = fPrimary[index-1] ; + } + + return rv ; + +} + +//____________________________________________________________________________ +Int_t AliEMCALDigit::GetIparent(Int_t index) const +{ + // retrieves the primary particle number given its index in the list + Int_t rv = -1 ; + if ( index <= fNiparent ){ + rv = fIparent[index-1] ; + } + + return rv ; + +} + + +//____________________________________________________________________________ +void AliEMCALDigit::ShiftPrimary(Int_t shift){ + //shifts primary nimber to BIG offset, to separate primary in different TreeK + Int_t index ; + for(index = 0; index fNMaxPrimary) { + cout << "AliEMCALDigit >> Increase NMaxPrimary "<< endl ; + return *this ; + } + } + } + + for (index = 0 ; index < digit.fNiparent ; index++){ + Bool_t dejavu = kTRUE ; + Int_t old ; + for ( old = 0 ; (old < max2) && dejavu; old++) { //already have this primary? + if(fIparent[old] == (digit.fIparent)[index]) + dejavu = kFALSE; + } + if(dejavu){ + fIparent[fNiparent] = (digit.fIparent)[index] ; + fNiparent++ ; + if(fNiparent>fNMaxiparent) { + cout << "AliEMCALDigit >> Increase NMaxiparent "<< endl ; + return *this ; + } + } + } + + return *this ; +} + +//____________________________________________________________________________ +ostream& operator << ( ostream& out , const AliEMCALDigit & digit) +{ + // Prints the data of the digit + + out << "ID " << digit.fId << " Energy = " << digit.fAmp << endl ; + Int_t i,j ; + for(i=0;iExecuteTask() +// Warning in : object already instantiated +// //Digitizes SDigitis in all events found in file galice.root +// +// root[2] AliEMCALDigitizer * d1 = new AliEMCALDigitizer("galice1.root") ; +// // Will read sdigits from galice1.root +// root[3] d1->MixWith("galice2.root") +// Warning in : object already instantiated +// // Reads another portion of sdigits from galice2.root +// root[3] d1->MixWith("galice3.root") +// // Reads another portion of sdigits from galice3.root +// root[4] d->ExecuteTask("deb timing") +// // Reads SDigits from files galice1.root, galice2.root .... +// // mixes them and stores produced Digits in file galice1.root +// // deb - prints number of produced digits +// // deb all - prints list of produced digits +// // timing - prints time used for digitization +// + +// --- ROOT system --- +#include "TFile.h" +#include "TTree.h" +#include "TSystem.h" +#include "TROOT.h" +#include "TFolder.h" +#include "TObjString.h" +#include "TBenchmark.h" +// --- Standard library --- +#include + +// --- AliRoot header files --- + +#include "AliRun.h" +#include "AliEMCALDigit.h" +#include "AliEMCALHit.h" +#include "AliEMCALv1.h" +#include "AliEMCALDigitizer.h" +#include "AliEMCALSDigitizer.h" +#include "AliEMCALGeometry.h" + +ClassImp(AliEMCALDigitizer) + + +//____________________________________________________________________________ + AliEMCALDigitizer::AliEMCALDigitizer():TTask("AliEMCALDigitizer","") +{ + // ctor + + fSDigitizer = 0 ; + fNinputs = 1 ; + fPinNoise = 0.01 ; + fEMCDigitThreshold = 0.01 ; + fInitialized = kFALSE ; + + fHeaderFiles = 0; + fSDigitsTitles = 0; + fSDigits = 0 ; + fDigits = 0; + +} +//____________________________________________________________________________ +void AliEMCALDigitizer::Init(){ +// Makes all memory allocations + + if(!fInitialized){ + + fHeaderFiles = new TClonesArray("TObjString",1) ; + new((*fHeaderFiles)[0]) TObjString("galice.root") ; + + //Test, if this file already open + + TFile *file = (TFile*) gROOT->GetFile(((TObjString *) fHeaderFiles->At(0))->GetString() ) ; + + if(file == 0){ + file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString(),"update") ; + gAlice = (AliRun *) file->Get("gAlice") ; + } + else + file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString()) ; + + file->cd() ; + + fSDigitsTitles = new TClonesArray("TObjString",1); + new((*fSDigitsTitles)[0]) TObjString("") ; + + fSDigits = new TClonesArray("TClonesArray",1) ; + new((*fSDigits)[0]) TClonesArray("AliEMCALDigit",1000) ; + + fSDigitizer = 0 ; + + fDigitsTitle = "" ; + + fDigits = new TClonesArray("AliEMCALDigit",200000) ; + + fIevent = new TArrayI(1) ; + fIevent->AddAt(-1,0 ) ; + fIeventMax = new TArrayI(1) ; + + fIeventMax->AddAt((Int_t) gAlice->TreeE()->GetEntries(), 0 ); + + // add Task to //root/Tasks folder + TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; + roottasks->Add(this) ; + + fInitialized = kTRUE ; + } + +} + +//____________________________________________________________________________ +AliEMCALDigitizer::AliEMCALDigitizer(const char *headerFile,const char *sDigitsTitle): + TTask("AliEMCALDigitizer","") +{ + // ctor + fHeaderFiles = new TClonesArray("TObjString",1) ; + new((*fHeaderFiles)[0]) TObjString(headerFile) ; + + // Header file, where result will be stored + TFile * file = (TFile*) gROOT->GetFile(((TObjString *) fHeaderFiles->At(0))->GetString() ) ; + if(file==0){ + if(((TObjString *) fHeaderFiles->At(0))->GetString().Contains("rfio")) + file = TFile::Open(((TObjString *) fHeaderFiles->At(0))->GetString(),"update") ; + else + file = new TFile(((TObjString *) fHeaderFiles->At(0))->GetString(),"update") ; + gAlice = (AliRun *) file->Get("gAlice") ; //If not read yet + } + + file->cd() ; + + fSDigitsTitles = new TClonesArray("TObjString",1); // Title name of the SDigits branch + new((*fSDigitsTitles)[0]) TObjString(sDigitsTitle) ; + + fSDigits = new TClonesArray("TClonesArray",1) ; // here list of SDigits wil be stored + new((*fSDigits)[0]) TClonesArray("AliEMCALDigit",1000) ; + + fDigits = new TClonesArray("AliEMCALDigit",200000) ; + + fDigitsTitle = "" ; + + + fSDigitizer = 0 ; + + fIevent = new TArrayI(1) ; + fIevent->AddAt(-1,0 ) ; + fIeventMax = new TArrayI(1) ; + + // Get number of events to process + fIeventMax->AddAt((Int_t) gAlice->TreeE()->GetEntries(), 0 ); + + fNinputs = 1 ; + + fPinNoise = 0.01 ; + fEMCDigitThreshold = 0.01 ; + + // add Task to //root/Tasks folder + TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; + roottasks->Add(this) ; + fInitialized = kTRUE ; + +} + +//____________________________________________________________________________ + AliEMCALDigitizer::~AliEMCALDigitizer() +{ + // dtor + + if(fHeaderFiles) delete fHeaderFiles ; + if(fSDigitsTitles) delete fSDigitsTitles ; + if(fSDigits) delete fSDigits ; + if(fDigits) delete fDigits ; +} +//____________________________________________________________________________ +void AliEMCALDigitizer::Reset() { + //sets current event number to the first simulated event + + if(!fInitialized) + Init() ; + + Int_t inputs ; + for(inputs = 0; inputs < fNinputs ;inputs++) + fIevent->AddAt(-1, inputs ) ; + +} +//____________________________________________________________________________ +Bool_t AliEMCALDigitizer::Combinator() { + + //Makes all desirable combinations Signal+Background, + // returns kFALSE when all combinations are made + // May be useful to introduce some options like "One-to-One", "All-to-One" and "All-to-All" ? + + //realizing "One-to-One" option... + + if(!fInitialized) + Init() ; + + Int_t inputs ; + Bool_t endNotReached = kTRUE ; + + for(inputs = 0; (inputs < fNinputs) && endNotReached ;inputs++){ + if(fIevent->At(inputs)+1 < fIeventMax->At(inputs)) + fIevent->AddAt(fIevent->At(inputs)+1, inputs ) ; + else + if(inputs == 0) + endNotReached = kFALSE ; + else //for inputs other than base one start from the beginning + fIevent->AddAt(0, inputs ) ; + + } + return endNotReached ; + +} + +//____________________________________________________________________________ +void AliEMCALDigitizer::Digitize(Option_t *option) { + + // Makes the digitization of the collected summable digits + // for this it first creates the array of all EMCAL modules + // filled with noise (different for EMC, CPV and PPSD) and + // after that adds contributions from SDigits. This design + // helps to avoid scanning over the list of digits to add + // contribution of any new SDigit. + + if(!fInitialized) + Init() ; + + fDigits->Clear() ; + + AliEMCAL * EMCAL = (AliEMCAL *) gAlice->GetDetector("EMCAL") ; + AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance( EMCAL->GetGeometry()->GetName(), EMCAL->GetGeometry()->GetTitle() ); + + //Making digits with noise, first EMC + Int_t nEMC = geom->GetNPhi()*geom->GetNZ(); + Int_t absID ; + TString name = geom->GetName() ; + + for(absID = 1; absID <= nEMC; absID++){ + Float_t noise = gRandom->Gaus(0., fPinNoise) ; + new((*fDigits)[absID-1]) AliEMCALDigit( -1,-1, absID,fSDigitizer->Digitize(noise) ) ; + } + + + // Now look throught (unsorted) list of SDigits and add corresponding digits + AliEMCALDigit *curSDigit ; + AliEMCALDigit *digit ; + + Int_t inputs; + for(inputs = 0; inputs< fNinputs ; inputs++){ //loop over (possible) merge sources + + TClonesArray * sdigits= (TClonesArray *)fSDigits->At(inputs) ; + Int_t isdigit ; + + Int_t nSDigits = sdigits->GetEntries() ; + for(isdigit=0;isdigit< nSDigits; isdigit++){ + curSDigit = (AliEMCALDigit *)sdigits->At(isdigit) ; + if(inputs) //Shift primaries for non-background sdigits + curSDigit->ShiftPrimary(inputs) ; + digit = (AliEMCALDigit *)fDigits->At(curSDigit->GetId() - 1); + *digit = *digit + *curSDigit ; + } + } + + + //remove digits below thresholds + for(absID = 0; absID < nEMC ; absID++) + if(fSDigitizer->Calibrate(((AliEMCALDigit*)fDigits->At(absID))->GetAmp()) < fEMCDigitThreshold) + fDigits->RemoveAt(absID) ; + + fDigits->Compress() ; + + Int_t ndigits = fDigits->GetEntriesFast() ; + + fDigits->Expand(ndigits) ; + + + //Set indexes in list of digits + Int_t i ; + for (i = 0 ; i < ndigits ; i++) { + AliEMCALDigit * digit = (AliEMCALDigit *) fDigits->At(i) ; + digit->SetIndexInList(i) ; + } +} +//____________________________________________________________________________ +void AliEMCALDigitizer::WriteDigits(){ + + // Made TreeD in the output file. Check if branch already exists: if yes, exits + // without writing: ROOT TTree does not suppert overwriting/updating of the + // already existing branches. Creates branch with Digits, named "EMCAL", title "...", + // and branch "AliEMCALDigitizer", with the same title to keep all the parameters + // and names of files, from which digits are made. + + gAlice->GetEvent(fIevent->At(0)) ; // Suitable only for One-To-One mixing + gAlice->SetEvent(fIevent->At(0)) ; // for all-to-all will produce a lot of branches in TreeD + + if(gAlice->TreeD()==0) + gAlice->MakeTree("D") ; + + //Check, if this branch already exits? + TBranch * digitsBranch = 0; + TBranch * digitizerBranch = 0; + + TObjArray * branches = gAlice->TreeD()->GetListOfBranches() ; + Int_t ibranch; + Bool_t emcalNotFound = kTRUE ; + Bool_t digitizerNotFound = kTRUE ; + + for(ibranch = 0;ibranch GetEntries();ibranch++){ + + if(emcalNotFound){ + digitsBranch=(TBranch *) branches->At(ibranch) ; + if( (strcmp("EMCAL",digitsBranch->GetName())==0 ) && + (fDigitsTitle.CompareTo(digitsBranch->GetTitle()) == 0) ) + emcalNotFound = kFALSE ; + } + if(digitizerNotFound){ + digitizerBranch = (TBranch *) branches->At(ibranch) ; + if( (strcmp(digitizerBranch->GetName(),"AliEMCALDigitizer") == 0) && + (fDigitsTitle.CompareTo(digitizerBranch->GetTitle()) == 0)) + digitizerNotFound = kFALSE ; + } + } + + + if(!(digitizerNotFound && emcalNotFound)){ + cout << "AliEMCALDigitizer error: " << endl ; + cout << " can not update/overwrite existing branches "<< endl ; + cout << " do not write " << endl ; + return ; + } + + // create new branches + + //First generate file name + char * file =0; + if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name + file = new char[strlen(gAlice->GetBaseFile())+20] ; + sprintf(file,"%s/EMCAL.Digits.root",gAlice->GetBaseFile()) ; + } + + TDirectory *cwd = gDirectory; + + //First create list of sdigits + Int_t bufferSize = 32000 ; + digitsBranch = gAlice->TreeD()->Branch("EMCAL",&fDigits,bufferSize); + digitsBranch->SetTitle(fDigitsTitle.Data()); + if (file) { + digitsBranch->SetFile(file); + TIter next( digitsBranch->GetListOfBranches()); + TBranch * sbr ; + while ((sbr=(TBranch*)next())) { + sbr->SetFile(file); + } + cwd->cd(); + } + + //second - create Digitizer + Int_t splitlevel = 0 ; + AliEMCALDigitizer * d = this ; + digitizerBranch = gAlice->TreeD()->Branch("AliEMCALDigitizer","AliEMCALDigitizer", + &d,bufferSize,splitlevel); + digitizerBranch->SetTitle(fDigitsTitle.Data()); + if (file) { + digitizerBranch->SetFile(file); + TIter next( digitizerBranch->GetListOfBranches()); + TBranch * sbr; + while ((sbr=(TBranch*)next())) { + sbr->SetFile(file); + } + cwd->cd(); + } + + digitsBranch->Fill() ; + digitizerBranch->Fill() ; + + gAlice->TreeD()->Write(0,kOverwrite) ; + + //remove fSDigitizer before new event. + if(fSDigitizer){ + delete fSDigitizer ; + fSDigitizer = 0 ; + } + + +} + +//____________________________________________________________________________ +void AliEMCALDigitizer::Exec(Option_t *option) { + // Managing method + + if(!fInitialized) Init() ; + + if(strstr(option,"tim")) + gBenchmark->Start("EMCALDigitizer"); + + //reset events numbers to start from the beginnig + Reset() ; + + while(Combinator()){ + + if(!ReadSDigits()) //read sdigits event(s) evaluated by Combinator() from file(s) + return ; + + Digitize(option) ; //Add prepared SDigits to digits and add the noise + WriteDigits() ; + + if(strstr(option,"deb")) + PrintDigits(option); + + } + + if(strstr(option,"tim")){ + gBenchmark->Stop("EMCALDigitizer"); + cout << "AliEMCALDigitizer:" << endl ; + cout << " took " << gBenchmark->GetCpuTime("EMCALDigitizer") << " seconds for SDigitizing " + << gBenchmark->GetCpuTime("EMCALDigitizer")/(fIeventMax->At(0)) << " seconds per event " << endl ; + cout << endl ; + } + +} + +//__________________________________________________________________ +Bool_t AliEMCALDigitizer::ReadSDigits(){ +// Reads summable digits from the opened files for the particular set of events given by fIevent + + if(!fInitialized) Init() ; + + + Int_t inputs ; + for(inputs = fNinputs-1; inputs >= 0; inputs --){ + + Int_t event = fIevent->At(inputs) ; + + TFile * file = (TFile*) gROOT->GetFile(((TObjString *) fHeaderFiles->At(inputs))->GetString() ) ; + file->cd() ; + + // Get SDigits Tree header from file + char treeName[20]; + sprintf(treeName,"TreeS%d",event); + TTree * treeS = (TTree*)file->Get(treeName); + + if(treeS==0){ + cout << "Error at AliEMCALDigitizer: no "<GetName() << endl ; + cout << "Do nothing " << endl ; + return kFALSE ; + } + + TBranch * sdigitsBranch = 0; + TBranch * sdigitizerBranch = 0; + + TObjArray * branches = 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 ) && + ((TObjString*) fSDigitsTitles->At(inputs))->GetString().CompareTo(sdigitsBranch->GetTitle())== 0 ) + emcalNotFound = kFALSE ; + + } + + if(sdigitizerNotFound){ + sdigitizerBranch = (TBranch *) branches->At(ibranch) ; + if(( strcmp(sdigitizerBranch->GetName(),"AliEMCALSDigitizer") == 0) && + ((TObjString*) fSDigitsTitles->At(inputs))->GetString().CompareTo(sdigitizerBranch->GetTitle())== 0 ) + sdigitizerNotFound = kFALSE ; + + } + } + + if(sdigitizerNotFound || emcalNotFound){ + cout << "Can't find Branch with sdigits or SDigitizer in the file " ; + if( ((TObjString*)fSDigitsTitles->At(inputs))->GetString().IsNull() ) + cout << file->GetName() << endl ; + else + cout << ((TObjString*)fSDigitsTitles->At(inputs))->GetString().Data() << endl ; + cout << "Do nothing" <At(inputs) ; + sdigitsBranch->SetAddress(&sdigits) ; + + AliEMCALSDigitizer *sDigitizer = new AliEMCALSDigitizer(); + sdigitizerBranch->SetAddress(&sDigitizer) ; + + sdigitsBranch->GetEntry(0) ; + sdigitizerBranch->GetEntry(0) ; + + if(fSDigitizer == 0) + fSDigitizer = sDigitizer ; + else + if(!((*fSDigitizer)==(*sDigitizer)) ){ + cout << "AliEMCALDigitizer ERROR:" << endl ; + cout << " you are using sdigits made with different SDigitizers" << endl ; + cout << "fSD " << fSDigitizer << " SD" << sDigitizer << endl ; + fSDigitizer->Print("") ; + sDigitizer->Print("") ; + cout << "Do Nothing " << endl ; + return kFALSE ; + } + + } + fPedestal = fSDigitizer->GetPedestalParameter() ; + fSlope = fSDigitizer->GetCalibrationParameter() ; + + return kTRUE ; + +} +//__________________________________________________________________ +void AliEMCALDigitizer::MixWith(char* HeaderFile, char* sDigitsTitle){ + // Alows produce digits by superimposing background and signal event. + // It is assumed, that headers file with SIGNAL events is opened in + // constructor, and now we set the BACKGROUND event, with which we + // will mix. Thus we avoid writing (changing) huge and expencive + // backgound files: all output will be writen into SIGNAL, i.e. + // opened in constructor file. + // + // One can open as many files to mix with as one wants. + + + if(!fInitialized) + Init() ; + + + if(HeaderFile == 0){ + cout << "Specify at least header file to merge"<< endl ; + return ; + } + + Int_t inputs ; + for(inputs = 0; inputs < fNinputs ; inputs++){ + if(strcmp(((TObjString *)fHeaderFiles->At(inputs))->GetString(),HeaderFile) == 0 ){ + if(sDigitsTitle == 0){ + if(((TObjString*)fSDigitsTitles->At(inputs))->GetString().CompareTo("") == 0){ + cout << "Entry already exists, do not add" << endl ; + return ; + } + } + else + if(((TObjString*)fSDigitsTitles->At(inputs))->GetString().CompareTo(sDigitsTitle)){ + cout << "Entry already exists, do not add" << endl ; + return; + } + } + } + + fHeaderFiles->Expand(fNinputs+1) ; + new((*fHeaderFiles)[fNinputs]) TObjString(HeaderFile) ; + + + TFile * file = new TFile(((TObjString *) fHeaderFiles->At(fNinputs))->GetString()) ; + + file->cd() ; + + fSDigitsTitles->Expand(fNinputs+1) ; + new((*fSDigitsTitles)[fNinputs]) TObjString(sDigitsTitle) ; + + fSDigits->Expand(fNinputs+1) ; + new((*fSDigits)[fNinputs]) TClonesArray("AliEMCALDigit",1000) ; + + fIevent->Set(fNinputs+1) ; + fIevent->AddAt(-1, fNinputs) ; + + fIeventMax->Set(fNinputs+1) ; + + TTree * te = (TTree *) file->Get("TE") ; + fIeventMax->AddAt((Int_t) te->GetEntries(), fNinputs ); + + fNinputs++ ; + +} +//__________________________________________________________________ +void AliEMCALDigitizer::Print(Option_t* option)const { + + if(fInitialized){ + + cout << "------------------- "<< GetName() << " -------------" << endl ; + cout << "Digitizing sDigits from file(s): " <At(input))->GetString() << + " Branch title:" << ((TObjString *) fSDigitsTitles->At(input))->GetString() << endl ; + } + cout << endl ; + cout << "Writing digits to " << ((TObjString *) fHeaderFiles->At(0))->GetString() << endl ; + + cout << endl ; + cout << "With following parameters: " << endl ; + cout << " Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ; + cout << " Threshold in EMC (fEMCDigitThreshold) = " << fEMCDigitThreshold << endl ; ; + cout << "---------------------------------------------------" << endl ; + } + else + cout << "AliEMCALDigitizer not initialized " << endl ; + +} +//__________________________________________________________________ +void AliEMCALDigitizer::PrintDigits(Option_t * option){ + + cout << "AliEMCALDigitiser:"<< endl ; + cout << " Number of entries in Digits list " << fDigits->GetEntriesFast() << endl ; + cout << endl ; + if(strstr(option,"all")){ + + //loop over digits + AliEMCALDigit * digit; + cout << "Digit Id " << " Amplitude " << " Index " << " Nprim " << " Primaries list " << endl; + Int_t index ; + for (index = 0 ; index < fDigits->GetEntries() ; index++) { + digit = (AliEMCALDigit * ) fDigits->At(index) ; + cout << setw(8) << digit->GetId() << " " << setw(3) << digit->GetAmp() << " " + << setw(6) << digit->GetIndexInList() << " " + << setw(5) << digit->GetNprimary() <<" "; + + Int_t iprimary; + for (iprimary=0; iprimaryGetNprimary(); iprimary++) + cout << setw(5) << digit->GetPrimary(iprimary+1) << " "; + cout << endl; + } + + } +} +//__________________________________________________________________ +void AliEMCALDigitizer::SetSDigitsBranch(const char* title){ + // we set title (comment) of the SDigits branch in the first! header file + if(!fInitialized) Init() ; + + ((TObjString*) fSDigitsTitles->At(0) )->SetString((char*)title) ; + +} +//__________________________________________________________________ +void AliEMCALDigitizer::SetDigitsBranch(const char* title){ + //Sets the title (comment) of the branch to which Digits branch + if(!fInitialized) Init() ; + + fDigitsTitle = title ; + +} diff --git a/EMCAL/AliEMCALDigitizer.h b/EMCAL/AliEMCALDigitizer.h new file mode 100644 index 00000000000..cb6978c1f7f --- /dev/null +++ b/EMCAL/AliEMCALDigitizer.h @@ -0,0 +1,83 @@ +#ifndef ALIEMCALDigitizer_H +#define ALIEMCALDigitizer_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//_________________________________________________________________________ +// Task Class for making SDigits in EMCAL +// +//*-- Author: Dmitri Peressounko(SUBATECH & KI) + + +// --- ROOT system --- +#include "TTask.h" +#include "TObjString.h" +class TArrayI ; +// --- Standard library --- + +// --- AliRoot header files --- +class AliEMCALSDigitizer ; + + +class AliEMCALDigitizer: public TTask { + +public: + AliEMCALDigitizer() ; // ctor + AliEMCALDigitizer(const char *headerFile,const char * sDigitsBranchTitle = 0) ; + virtual ~AliEMCALDigitizer() ; + + void Digitize(Option_t *option); // Make Digits from SDigits stored in fSDigits + void Exec(Option_t *option); // Supervising method + + Float_t GetEMCThreshold() const { return fEMCDigitThreshold;} + Float_t GetPedestal() const { return fPedestal; } + Float_t GetPinNoise() const { return fPinNoise;} + Float_t GetSlope() const { return fSlope; } + char * GetDigitsBranch ()const { return (char*)fDigitsTitle.Data() ;} + TClonesArray * GetHeadersFiles(){ return fHeaderFiles ;} + TArrayI* GetCurrentEvents() { return fIevent ;} + + void MixWith(char* HeaderFile, char* SDigitsTitle =0) ; // Add another one file to mix + virtual void Print(Option_t* option)const ; + void Reset() ; //restarts starts event processing from 0 event(s) + + void SetEMCThreshold(Float_t EMCThreshold) {fEMCDigitThreshold = EMCThreshold;} + void SetPinNoise(Float_t PinNoise ) {fPinNoise = PinNoise;} + + void SetDigitsBranch (const char* file) ; + void SetSDigitsBranch(const char* file) ; + +private: + Bool_t Combinator() ; // makes all desirable combination sig+Bg + void Init(); + Bool_t ReadSDigits() ; // Read sdigits for particular events + void WriteDigits() ; // Writes Digits for particular event + void PrintDigits(Option_t * option) ; + +private: + TClonesArray * fSDigitsTitles ; // Titles of sdigits branches + TClonesArray * fHeaderFiles ; // Names of files with headers to merge + TString fDigitsTitle ; // Title of the Digits Branch + TClonesArray * fSDigits ; // ! Lists of SDigits + TClonesArray * fDigits ; // ! Final list of digits + AliEMCALSDigitizer * fSDigitizer ; // ! SDigitizer to extarct some sdigitizing parameters + Int_t fNinputs ; // Number of input files + Bool_t fInitialized ; // + TArrayI * fIevent ; // events to read at the next ReadSDigits() call + TArrayI * fIeventMax ; // Maximal number of events in each input file + + Float_t fPedestal ; // Calibration parameters + Float_t fSlope ; // read from SDigitizer + + Float_t fPinNoise ; // Electronics noise in EMC + Float_t fEMCDigitThreshold ; // Threshold for storing digits in EMC + + + ClassDef(AliEMCALDigitizer,1) // description + +}; + + +#endif // AliEMCALDigitizer_H diff --git a/EMCAL/AliEMCALHit.cxx b/EMCAL/AliEMCALHit.cxx index 5616adf23f8..f0f522c2fef 100644 --- a/EMCAL/AliEMCALHit.cxx +++ b/EMCAL/AliEMCALHit.cxx @@ -51,7 +51,12 @@ AliEMCALHit::AliEMCALHit(){ fX = 0.0; fY = 0.0; fZ = 0.0; - fP = 0.0; + fPx = 0.0; + fPy = 0.0; + fPz = 0.0; + fPe = 0.0; + fIparent = 0; + fIenergy = 0.0; } //______________________________________________________________________ AliEMCALHit::AliEMCALHit(const AliEMCALHit & hit){ @@ -64,22 +69,31 @@ AliEMCALHit::AliEMCALHit(const AliEMCALHit & hit){ fX = hit.fX; fY = hit.fY; fZ = hit.fZ; - fP = hit.fP; + fPx = hit.fPx; + fPy = hit.fPy; + fPz = hit.fPz; + fPe = hit.fPe; + fIparent = hit.fIparent; + fIenergy = hit.fIenergy; } //______________________________________________________________________ -AliEMCALHit::AliEMCALHit(Int_t shunt, Int_t primary, Int_t track,Int_t id, - Float_t *hits,TLorentzVector *p):AliHit(shunt, track){ +AliEMCALHit::AliEMCALHit(Int_t shunt, Int_t primary, Int_t track,Int_t iparent, Float_t ienergy, Int_t id, + Float_t *hits,Float_t *p):AliHit(shunt, track){ // // Create a CPV hit object // - fX = hits[0]; fY = hits[1]; fZ = hits[2]; fId = id; fELOS = hits[3]; fPrimary = primary; - fP = *p; + fPx = p[0]; + fPy = p[1]; + fPz = p[2]; + fPe = p[3]; + fIparent = iparent; + fIenergy = ienergy; } //______________________________________________________________________ Bool_t AliEMCALHit::operator==(AliEMCALHit const &rValue) const{ @@ -87,7 +101,7 @@ Bool_t AliEMCALHit::operator==(AliEMCALHit const &rValue) const{ // from the same primary Bool_t rv = kFALSE; - if ( (fId == rValue.GetId()) && ( fPrimary == rValue.GetPrimary()) ) + if ( (fId == rValue.GetId()) && ( fPrimary == rValue.GetIparent()) ) rv = kTRUE; return rv; @@ -111,8 +125,10 @@ ostream& operator << (ostream& out,AliEMCALHit& hit){ out << "GeV , Track no.=" << hit.GetPrimary(); out << ", (xyz)=(" << hit.X()<< ","<< hit.Y()<< ","< +//class ostream; class AliEMCALHit : public AliHit { friend ostream& operator << (ostream&,AliEMCALHit&); public: AliEMCALHit(); // default ctor AliEMCALHit(const AliEMCALHit & hit); - AliEMCALHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t id, - Float_t *hits,TLorentzVector *p); + AliEMCALHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, Int_t id, + Float_t *hits,Float_t *p); virtual ~AliEMCALHit(void) {}// dtor //returns the energy loss for this hit Float_t GetEnergy(void) const{return fELOS;} @@ -33,7 +33,13 @@ class AliEMCALHit : public AliHit { // returns the primary particle id at the origine of this hit Int_t GetPrimary(void) const{return fPrimary;} // returns the energy/momentum LorentzVector of the enetering particle. - TLorentzVector& GetP(void) {return fP;} + Int_t GetIparent(void) const{return fIparent;} + Float_t GetIenergy(void) const{return fIenergy;} + + Float_t GetPx(void) const{return fPx;} + Float_t GetPy(void) const{return fPy;} + Float_t GetPz(void) const{return fPz;} + Float_t GetPe(void) const{return fPe;} Bool_t operator == (AliEMCALHit const &rValue) const; AliEMCALHit operator + (const AliEMCALHit& rValue); @@ -41,9 +47,13 @@ class AliEMCALHit : public AliHit { Int_t fId; // Absolute Id number EMCAL segment Float_t fELOS; // Energy deposited Int_t fPrimary; // Primary particles at the origine of the hit - TLorentzVector fP; // Primary partical enetrence momentum/energy - - ClassDef(AliEMCALHit,1) // Hit for EMCAL + Float_t fPx; // Primary partical enetrence momentum/energy + Float_t fPy; // Primary partical enetrence momentum/energy + Float_t fPz; // Primary partical enetrence momentum/energy + Float_t fPe; // Primary partical enetrence momentum/energy + Int_t fIparent; // Parent particle that enterred emcal + Float_t fIenergy; // Initial energy of parent particle that enterred the emcal + ClassDef(AliEMCALHit,2) // Hit for EMCAL }; #endif // ALIEMCALHIT_H diff --git a/EMCAL/AliEMCALLinkDef.h b/EMCAL/AliEMCALLinkDef.h deleted file mode 100644 index c43a0e76375..00000000000 --- a/EMCAL/AliEMCALLinkDef.h +++ /dev/null @@ -1,12 +0,0 @@ -#ifdef __CINT__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ class AliEMCAL+; -#pragma link C++ class AliEMCALGeometry+; -#pragma link C++ class AliEMCALv0+; - -#endif - diff --git a/EMCAL/AliEMCALSDigitizer.cxx b/EMCAL/AliEMCALSDigitizer.cxx new file mode 100644 index 00000000000..0961a29c87f --- /dev/null +++ b/EMCAL/AliEMCALSDigitizer.cxx @@ -0,0 +1,398 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//_________________________________________________________________________ +// This is a TTask that makes SDigits out of Hits +// A Summable Digits is the sum of all hits originating +// from one primary in one active cell +// A threshold for assignment of the primary to SDigit is applied +// SDigits are written to TreeS, branch "EMCAL" +// AliEMCALSDigitizer with all current parameters is written +// to TreeS branch "AliEMCALSDigitizer". +// Both branches have the same title. If necessary one can produce +// another set of SDigits with different parameters. Two versions +// can be distunguished using titles of the branches. +// User case: +// root [0] AliEMCALSDigitizer * s = new AliEMCALSDigitizer("galice.root") +// Warning in : object already instantiated +// root [1] s->ExecuteTask() +// // 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") +// // available parameters: +// deb - print # of produced SDigitis +// deb all - print # and list of produced SDigits +// tim - print benchmarking information +// +//*-- Author : Dmitri Peressounko (SUBATECH & KI) +////////////////////////////////////////////////////////////////////////////// + + +// --- ROOT system --- +#include "TFile.h" +#include "TTask.h" +#include "TTree.h" +#include "TSystem.h" +#include "TROOT.h" +#include "TFolder.h" +#include "TBenchmark.h" +// --- Standard library --- +#include + +// --- AliRoot header files --- +#include "AliRun.h" +#include "AliEMCALDigit.h" +#include "AliEMCALHit.h" +#include "AliEMCALSDigitizer.h" +#include "AliEMCALGeometry.h" +#include "AliEMCALv1.h" + +ClassImp(AliEMCALSDigitizer) + + +//____________________________________________________________________________ + AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("AliEMCALSDigitizer","") +{ + // ctor + fA = 0; + fB = 10000000. ; + fPrimThreshold = 0.001 ; + fNevents = 0 ; + fSDigits = 0 ; + fHits = 0 ; + fIsInitialized = kFALSE ; + +} + +//____________________________________________________________________________ +AliEMCALSDigitizer::AliEMCALSDigitizer(const char* headerFile, const char *sDigitsTitle):TTask("AliEMCALSDigitizer","") +{ + // ctor + fA = 0; + fB = 10000000.; + fPrimThreshold = 0.001 ; + fNevents = 0 ; + fSDigitsTitle = sDigitsTitle ; + fHeadersFile = headerFile ; + fSDigits = new TClonesArray("AliEMCALDigit",1000); + fHits = new TClonesArray("AliEMCALHit",1000); + + 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() +{ + // dtor + if(fSDigits) + delete fSDigits ; + if(fHits) + delete fHits ; +} +//____________________________________________________________________________ +void AliEMCALSDigitizer::Init(){ + //Initialization can not be done in the default constructor + + if(!fIsInitialized){ + + if(fHeadersFile.IsNull()) + fHeadersFile="galice.root" ; + + 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 ; + } +} +//____________________________________________________________________________ +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); + AliEMCAL * EMCAL = (AliEMCAL *) 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 = AliEMCALDigit(1,1,1,1); + AliEMCALDigit *sdigit ; + + // Assign primary number only if contribution is significant + + if( hit->GetEnergy() > fPrimThreshold) + curSDigit = AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; + else + curSDigit = AliEMCALDigit( -1 , -1 ,(((hit->GetId()/geom->GetNPhi())%geom->GetNZ()+1) * (hit->GetId()%(geom->GetNPhi()+1))), Digitize( hit->GetEnergy() ) ) ; + + for(Int_t check= 0; check < nSdigits; check++) { + sdigit = (AliEMCALDigit *)fSDigits->At(check); + if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ())) && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi()))) + { + *sdigit = *sdigit + curSDigit ; + } + else + { new((*fSDigits)[nSdigits]) AliEMCALDigit(curSDigit); + nSdigits++ ; } + } + + + if( hit->GetEnergy() > fPrimThreshold) + curSDigit = 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 = 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())) + { + for(Int_t check= 0; check < nSdigits; check++) { + sdigit = (AliEMCALDigit *)fSDigits->At(check); + if( (((hit->GetId()/geom->GetNPhi())%geom->GetNZ()) == ((sdigit->GetId()/geom->GetNPhi())%geom->GetNZ())) && ((hit->GetId()%geom->GetNPhi()) == (sdigit->GetId()%geom->GetNPhi()))) + { + *sdigit = *sdigit + curSDigit ; + } + else + { new((*fSDigits)[nSdigits]) AliEMCALDigit(curSDigit); + nSdigits++ ; } + } + } + } + } // 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) ; + } + + 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(!(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; + } + + sdigitsBranch->Fill() ; + sdigitizerBranch->Fill() ; + gAlice->TreeS()->Write(0,TObject::kOverwrite) ; + + if(strstr(option,"deb")) + PrintSDigits(option) ; + + } + + 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 ; + } + + +} +//__________________________________________________________________ +void AliEMCALSDigitizer::SetSDigitsBranch(const char * title ){ + //Seting title to branch SDigits + if(!fSDigitsTitle.IsNull()) + cout << "AliEMCALSdigitizer: changing SDigits file from " <GetEntriesFast() << endl ; + cout << endl ; + + if(strstr(option,"all")){// print all digits + + //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() <<" "; + + Int_t iprimary; + for (iprimary=0; iprimaryGetNprimary(); iprimary++) + cout << setw(5) << digit->GetPrimary(iprimary+1) << " "; + cout << endl; + } + + } +} diff --git a/EMCAL/AliEMCALSDigitizer.h b/EMCAL/AliEMCALSDigitizer.h new file mode 100644 index 00000000000..543b9e45c81 --- /dev/null +++ b/EMCAL/AliEMCALSDigitizer.h @@ -0,0 +1,64 @@ +#ifndef ALIEMCALSDigitizer_H +#define ALIEMCALSDigitizer_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//_________________________________________________________________________ +// Task Class for making SDigits in EMCAL +// +//*-- Author: Dmitri Peressounko(SUBATECH & KI) + + +// --- ROOT system --- +#include "TTask.h" +#include "TString.h" +// --- Standard library --- + +// --- AliRoot header files --- + +class AliEMCALSDigitizer: public TTask { + +public: + AliEMCALSDigitizer() ; // ctor + AliEMCALSDigitizer(const char* HeaderFile,const char *SdigitsTitle = 0) ; + virtual ~AliEMCALSDigitizer() ; // dtor + + Float_t Calibrate(Int_t amp)const {return (amp - fA)/fB ; } + Int_t Digitize(Float_t Energy)const { return (Int_t ) ( fA + Energy*fB); } + + virtual void Exec(Option_t *option); + + Float_t GetPedestalParameter()const {return fA;} + Float_t GetCalibrationParameter()const{return fB;} + char * GetSDigitsBranch()const{return (char*) fSDigitsTitle.Data();} + + virtual void Print(Option_t* option) const ; + + void SetPedestalParameter(Float_t A){fA = A ;} + void SetSlopeParameter(Float_t B){fB = B ;} + void SetSDigitsBranch(const char * title ) ; + + Bool_t operator == (const AliEMCALSDigitizer & sd) const ; + +private: + void Init() ; + void PrintSDigits(Option_t * option) ; + +private: + Float_t fA ; //Pedestal parameter + Float_t fB ; //Slope Digitizition parameters + Int_t fNevents ; // Number of events to digitize + Float_t fPrimThreshold ; // To store primari if Elos > threshold + TString fSDigitsTitle ; // title of SDigits branch + TString fHeadersFile ; //input file + Bool_t fIsInitialized ; + TClonesArray * fSDigits ; //! list of SDigits + TClonesArray * fHits ; //! + + ClassDef(AliEMCALSDigitizer,1) // description + +}; + +#endif // AliEMCALSDigitizer_H diff --git a/EMCAL/AliEMCALv0.cxx b/EMCAL/AliEMCALv0.cxx index a8f75449d7c..112577e883f 100644 --- a/EMCAL/AliEMCALv0.cxx +++ b/EMCAL/AliEMCALv0.cxx @@ -56,6 +56,7 @@ AliEMCALv0::AliEMCALv0(const char *name, const char *title): if (strcmp(GetTitle(),"") != 0 ) fGeom = AliEMCALGeometry::GetInstance(GetTitle(), "") ; + } //______________________________________________________________________ void AliEMCALv0::BuildGeometry(){ @@ -114,12 +115,12 @@ void AliEMCALv0::CreateGeometry(){ envelopC[3] = 2; envelopB[3] = 2; - envelopB[4] = (fGeom->GetEnvelop(0) + fGeom->GetGap2Active() + 1.59) / - (tan(2*atan(exp(-0.7)))) ; + envelopB[4] = (fGeom->GetEnvelop(0) + fGeom->GetGap2Active()) / + (tan(2*atan(exp(0.7)))) ; envelopB[5] = fGeom->GetEnvelop(0) + fGeom->GetGap2Active(); //rmin envelopD[6] = envelopB[6] = envelopB[5] + 3.18; //rmax - envelopB[7] = (fGeom->GetEnvelop(0) + fGeom->GetGap2Active()+ 1.59) / - (tan(2*atan(exp(0.7)))) ; + envelopB[7] = (fGeom->GetEnvelop(0) + fGeom->GetGap2Active()) / + (tan(2*atan(exp(-0.7)))) ; envelopB[8] = envelopB[5] ; envelopB[9] = envelopB[6] ; @@ -129,14 +130,14 @@ void AliEMCALv0::CreateGeometry(){ gMC->Gsvolu("XPST", "PGON", idtmed[1601], 0, 0) ; gMC->Gsvolu("XPBX", "PGON", idtmed[1600], 0, 0) ; // filled with Lead gMC->Gsdvn("XPHI", "XPST", fGeom->GetNPhi(), 2) ; // Naming Phi divisions - + Int_t idrotm = 1; AliMatrix(idrotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ; // Position ENV1 container in ALIC - gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ; + gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "MANY") ; // Position ARM1 into ENV1 - gMC->Gspos("XALU", 1, "XEN1", 0.0, 0.0, 0.0 , idrotm, "MANY") ; + gMC->Gspos("XALU", 1, "XEN1", 0.0, 0.0, 0.0 , idrotm, "ONLY") ; for (int i = 0; i < (fGeom->GetNLayers()); i++ ){ envelopC[5] = envelopD[6] ; //rmin @@ -165,7 +166,7 @@ void AliEMCALv0::CreateGeometry(){ 0.0, 0.0, 0.0 , idrotm, "MANY", envelopD, 10) ; } // end for j } // end if i - } // for i + } // for i } //______________________________________________________________________ void AliEMCALv0::Init(void){ diff --git a/EMCAL/AliEMCALv1.cxx b/EMCAL/AliEMCALv1.cxx index ba5cddc51b9..c401e430042 100644 --- a/EMCAL/AliEMCALv1.cxx +++ b/EMCAL/AliEMCALv1.cxx @@ -80,8 +80,8 @@ AliEMCALv1::~AliEMCALv1(){ } } //______________________________________________________________________ -void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, - Int_t id, Float_t * hits,TLorentzVector *p){ +void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, + Int_t id, Float_t * hits,Float_t * p){ // Add a hit to the hit list. // A PHOS hit is the sum of all hits in a single crystal // or in a single PPSD gas cell @@ -91,11 +91,10 @@ void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, AliEMCALHit *curHit; Bool_t deja = kFALSE; - newHit = new AliEMCALHit(shunt, primary, tracknumber, id, hits, p); - + newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p); for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) { curHit = (AliEMCALHit*) (*fHits)[hitCounter]; - // We add hits with the same primary, while GEANT treats + // We add hits with the same tracknumber, while GEANT treats // primaries succesively if(curHit->GetPrimary() != primary) break; if( *curHit == *newHit ) { @@ -119,26 +118,41 @@ void AliEMCALv1::StepManager(void){ Int_t absid; // position wrt MRS and energy deposited Float_t xyze[4]={0.,0.,0.,0.}; + Float_t pmom[4]={0.,0.,0.,0.}; TLorentzVector pos; // Lorentz vector of the track current position. TLorentzVector mom; // Lorentz vector of the track current momentum. Int_t tracknumber = gAlice->CurrentTrack(); - static Int_t primary; + Int_t primary; + static Int_t iparent = 0; + static Float_t ienergy = 0; + Int_t copy; + + + if(gMC->IsTrackEntering() && (strcmp(gMC->CurrentVolName(),"XALU") == 0)) + { + iparent = tracknumber; + gMC->TrackMomentum(mom); + ienergy = mom[3]; +} + if(gMC->CurrentVolID(copy) == gMC->VolId("XPHI") ) { // We are inside a PBWO crysta gMC->CurrentVolOffID(1, id[0]); // get the POLY copy number; gMC->CurrentVolID(id[1]); // get the phi number inside the layer primary = gAlice->GetPrimary(tracknumber); - if(gMC->IsTrackEntering()&&id[0]==1) primary = tracknumber; -// tracknumber = primary; gMC->TrackPosition(pos); gMC->TrackMomentum(mom); xyze[0] = pos[0]; xyze[1] = pos[1]; xyze[2] = pos[2]; xyze[3] = gMC->Edep(); - + pmom[0] = mom[0]; + pmom[1] = mom[1]; + pmom[2] = mom[2]; + pmom[3] = mom[3]; + if(xyze[3] > 0.){// Track is inside the crystal and deposits some energy absid = (id[0]-1)*(fGeom->GetNPhi()) + id[1]; - AddHit(fIshunt, primary,tracknumber, absid, xyze, &mom); + AddHit(fIshunt, primary,tracknumber, iparent, ienergy, absid, xyze, pmom); } // there is deposited energy - +} } diff --git a/EMCAL/AliEMCALv1.h b/EMCAL/AliEMCALv1.h index 1e94de26cc3..186c0b3fe5e 100644 --- a/EMCAL/AliEMCALv1.h +++ b/EMCAL/AliEMCALv1.h @@ -28,8 +28,8 @@ class AliEMCALv1 : public AliEMCALv0 { // requested by the Coding Convention AliEMCALv1(const AliEMCALv0 & emcal) {abort();} virtual ~AliEMCALv1(void) ; - virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, - Int_t id, Float_t *hits, TLorentzVector *p); + virtual void AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t iparent, Float_t ienergy, + Int_t id, Float_t *hits, Float_t *p); // Gives the version number virtual Int_t IsVersion(void) const {return 1;} virtual void StepManager(void) ; diff --git a/EMCAL/EMCALLinkDef.h b/EMCAL/EMCALLinkDef.h index ad897608ebd..2655c88e4a1 100644 --- a/EMCAL/EMCALLinkDef.h +++ b/EMCAL/EMCALLinkDef.h @@ -9,4 +9,7 @@ #pragma link C++ class AliEMCALv0+; #pragma link C++ class AliEMCALv1+; #pragma link C++ class AliEMCALHit+; +#pragma link C++ class AliEMCALDigit+; +#pragma link C++ class AliEMCALSDigitizer+; +#pragma link C++ class AliEMCALDigitizer+; #endif diff --git a/EMCAL/Makefile b/EMCAL/Makefile index 3913ac2baa2..c423be621f7 100644 --- a/EMCAL/Makefile +++ b/EMCAL/Makefile @@ -10,7 +10,8 @@ PACKAGE = EMCAL # C++ sources SRCS = AliEMCAL.cxx AliEMCALGeometry.cxx AliEMCALv0.cxx \ - AliEMCALv1.cxx AliEMCALHit.cxx + AliEMCALv1.cxx AliEMCALHit.cxx AliEMCALDigit.cxx \ + AliEMCALSDigitizer.cxx AliEMCALDigitizer.cxx # C++ Headers diff --git a/EMCAL/libEMCAL.pkg b/EMCAL/libEMCAL.pkg index 8647874cc64..fd7826cc3cb 100644 --- a/EMCAL/libEMCAL.pkg +++ b/EMCAL/libEMCAL.pkg @@ -3,13 +3,19 @@ AliEMCAL.cxx \ AliEMCALGeometry.cxx \ AliEMCALHit.cxx \ AliEMCALv0.cxx \ -AliEMCALv1.cxx +AliEMCALv1.cxx \ +AliEMCALDigit.cxx \ +AliEMCALSDigitizer.cxx \ +AliEMCALDigitizer.cxx HDRS= \ AliEMCAL.h \ AliEMCALGeometry.h \ AliEMCALHit.h \ AliEMCALv0.h \ -AliEMCALv1.h +AliEMCALv1.h \ +AliEMCALDigit.h \ +AliEMCALSDigitizer.h \ +AliEMCALDigitizer.h DHDR:=EMCALLinkDef.h -- 2.39.3