#include "AliHLTEMCALUtils.h" #include "TClass.h" #include "TFile.h" #include "TGeoManager.h" #include "TUnixSystem.h" #include "TTree.h" #include "TClonesArray.h" #include "TBranch.h" #include "TArrayS.h" #include "AliCDBManager.h" #include "AliCDBEntry.h" #include "AliRawReaderMemory.h" #include "AliEMCALRecParam.h" #include "AliEMCALReconstructor.h" #include "AliEMCALGeometry.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALRawUtils.h" #include "AliESDCaloCells.h" #include "AliESDCaloCluster.h" #include "AliEMCALDigit.h" #include "AliESDEvent.h" #include "AliEMCALRecPoint.h" #include "AliEMCALPID.h" #include #include using namespace std; ClassImp(AliHLTEMCALUtils); Int_t AliHLTEMCALUtils::fgDebug = 0; AliEMCALGeometry* AliHLTEMCALUtils::fgGeom = NULL; AliEMCALClusterizerv1* AliHLTEMCALUtils::fgClusterizer = NULL; AliEMCALRecParam* AliHLTEMCALUtils::fgRecParam = NULL; AliEMCALRawUtils* AliHLTEMCALUtils::fgRawUtils = NULL; TFile* AliHLTEMCALUtils::fgGeometryFile = NULL; TGeoManager* AliHLTEMCALUtils::fgGeoManager = NULL; TTree* AliHLTEMCALUtils::fgClustersTreeTmp = NULL; TTree* AliHLTEMCALUtils::fgDigitsTreeTmp = NULL; TClonesArray* AliHLTEMCALUtils::fgDigitsArrTmp = NULL; TBranch* AliHLTEMCALUtils::fgDigitsBranchTmp = NULL; //____________________________________________________________________________ AliHLTEMCALUtils::AliHLTEMCALUtils() : TObject() { // // Default constructor // ; } //____________________________________________________________________________ AliHLTEMCALUtils::~AliHLTEMCALUtils() { // // Destructor // ; } //____________________________________________________________________________ void AliHLTEMCALUtils::Cleanup() { DeleteStaticMembers(); DeleteReconstructionTrees(); } //____________________________________________________________________________ void AliHLTEMCALUtils::DeleteStaticMembers() { if (fgRecParam != AliEMCALReconstructor::GetRecParam()) { delete fgRecParam; fgRecParam = NULL; } delete fgRawUtils; fgRawUtils = NULL; DeleteGeoManager(); delete fgGeom; fgGeom = NULL; delete fgClusterizer; fgClusterizer = NULL; } //____________________________________________________________________________ AliHLTEMCALUtils::AliHLTEMCALUtils(const AliHLTEMCALUtils & /*t*/) : TObject() { // // copy ctor not to be used // AliFatal("May not use."); } //____________________________________________________________________________ AliHLTEMCALUtils& AliHLTEMCALUtils::operator = (const AliHLTEMCALUtils & /*t*/) { // // assignement operator not to be used // AliFatal("May not use.") ; return *this ; } //____________________________________________________________________________ void AliHLTEMCALUtils::InitRecParam() { // // Please check the AliEMCALReconstructor for comparison // Check if the instance of AliEMCALRecParam exists, // if not, get it from OCDB if available, otherwise create a default one // fgRecParam = (AliEMCALRecParam*) AliEMCALReconstructor::GetRecParam(); if (fgRecParam) return; if (!fgRecParam && (AliCDBManager::Instance()->IsDefaultStorageSet())) { AliCDBEntry *entry = (AliCDBEntry*) AliCDBManager::Instance()->Get("EMCAL/Config/RecParam"); if (entry) fgRecParam = (AliEMCALRecParam*) entry->GetObject(); } if(!fgRecParam) { AliWarningClass("The Reconstruction parameters initialized to default."); fgRecParam = new AliEMCALRecParam; } if (!fgRecParam) { AliErrorClass("Unable to init the reco params. Something is really wrong. Memory?"); } AliEMCALReconstructor::SetRecParam(fgRecParam); } //____________________________________________________________________________ AliEMCALRecParam* AliHLTEMCALUtils::GetRecParam() { // // Init the parameters and reuse the Reconstructor // AliHLTEMCALUtils::InitRecParam(); return fgRecParam; } //____________________________________________________________________________ AliEMCALRawUtils* AliHLTEMCALUtils::GetRawUtils(AliEMCALGeometry *pGeometry) { // // Init EMCAL raw utils // Use the geometry or try to figure out the default // The OCDB must be initialized. // Otherwise the default contructor of the raw utils will complain. // if (fgRawUtils == NULL) { if (AliCDBManager::Instance()->IsDefaultStorageSet()) { if (pGeometry) { fgRawUtils = new AliEMCALRawUtils(pGeometry); } else { fgRawUtils = new AliEMCALRawUtils(AliHLTEMCALUtils::GetGeometry()); } AliInfoClass("Raw Utils initialized."); } else { AliErrorClass("OCDB not initialized. Unable to init raw utils."); } } return fgRawUtils; } //____________________________________________________________________________ AliEMCALClusterizerv1* AliHLTEMCALUtils::GetClusterizer() { // // Init EMCAL clusterizer // if (fgClusterizer == NULL) { AliEMCALGeometry *geom = GetGeometry(); fgClusterizer = new AliEMCALClusterizerv1(geom); AliInfoClass("ClusterizerV1 initialized."); } return fgClusterizer; } //____________________________________________________________________________ AliEMCALGeometry* AliHLTEMCALUtils::GetGeometry() { // // Init EMCAL geometry // if (fgGeom == NULL) { AliInfoClass(Form("Using default geometry")); fgGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); } return fgGeom; } //____________________________________________________________________________ void AliHLTEMCALUtils::DeleteGeoManager() { // // Delete the geom manager and the geom TFile // if (fgGeometryFile) { fgGeometryFile->Close(); delete fgGeometryFile; fgGeometryFile = NULL; } if (fgGeoManager) { delete fgGeoManager; fgGeoManager = NULL; } } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::LoadGeoManagerFromFile(const char *fname) { // // Open the geom file and get the geom manager // Bool_t bReturn = kFALSE; DeleteGeoManager(); fgGeometryFile = TFile::Open(fname); if (fgGeometryFile) { fgGeoManager = (TGeoManager *)fgGeometryFile->Get("Geometry"); } if (!fgGeoManager) { DeleteGeoManager(); AliErrorClass(Form("Unable to retrieve TGeoManager from %s", fname)); } bReturn = kTRUE; return bReturn; } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::InitFakeCDB(const char *cdbpath, Int_t runid) { // // Init fake OCDB - use with care. // For dummy usage only! // Make sure there is no conflict with the existing ocdb instance. // Bool_t bReturn = kFALSE; TString sPath = cdbpath; if (sPath.Length() <= 0) { sPath = "local://"; sPath += gSystem->Getenv("ALICE_ROOT"); AliInfoClass(Form("Setting cdbpath to %s", sPath.Data())); } AliCDBManager* pCDB = AliCDBManager::Instance(); if (!pCDB) { AliErrorClass("Unable to get the CDBManager instance."); } else { pCDB->SetRun(runid); // THIS HAS TO BE RETRIEVED !!! pCDB->SetDefaultStorage(sPath.Data()); AliInfoClass(Form("Fake CDB initialized with path %s", sPath.Data())); bReturn = kTRUE; } return bReturn; } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree, Option_t* sOption) { // This method is a reusage the corresponding one of AliEMCALReconstructor // Conversion from raw data to EMCAL digits. // Works on a single-event basis // Make sure raw utils are initialized (they depend on geometry and OCDB) ! // We Assume the rawReader, digitsTree and sOption are valid! // Here we assume the raw reader is setup properly, so we do not need extra reset! // rawReader->Reset() ; Bool_t bReturn = kFALSE; if (!fgRawUtils) { // Try to initialize - this is potentially dangerous // We will use default geometry - most probably // The OCDB must exist! otherwise we will fail. if (GetRawUtils() == NULL) { // This is no go. AliErrorClass("Failed. No RawUtils."); return bReturn; } } if (!fgRecParam) { // Again we try to init the RecParams to default // One should better do it properly before! // The ocdb has to be there! for the recparams init if (GetRecParam() == NULL) { // This is no go. AliErrorClass("Failed. No RecParams."); return bReturn; } } TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",200); Int_t bufsize = 32000; digitsTree->Branch("EMCAL", &digitsArr, bufsize); fgRawUtils->SetOption(sOption); fgRawUtils->SetRawFormatHighLowGainFactor(fgRecParam->GetHighLowGainFactor()); fgRawUtils->SetRawFormatOrder(fgRecParam->GetOrderParameter()); fgRawUtils->SetRawFormatTau(fgRecParam->GetTau()); fgRawUtils->SetNoiseThreshold(fgRecParam->GetNoiseThreshold()); fgRawUtils->SetNPedSamples(fgRecParam->GetNPedSamples()); fgRawUtils->Raw2Digits(rawReader,digitsArr); digitsTree->Fill(); digitsArr->Delete(); delete digitsArr; bReturn = kTRUE; return bReturn; } //____________________________________________________________________________ TTree *AliHLTEMCALUtils::GetDigitsTree() { // // Create the digits tree and the digits arrays // if (fgDigitsTreeTmp == NULL) { fgDigitsTreeTmp = new TTree("EMCALdigits", "EMCAL digits default"); } if (fgDigitsArrTmp == NULL) { fgDigitsArrTmp = new TClonesArray("AliEMCALDigit",200); } if (fgDigitsBranchTmp == NULL) { Int_t bufsize = 32000; fgDigitsBranchTmp = fgDigitsTreeTmp->Branch("EMCAL", &fgDigitsArrTmp, bufsize); } fgDigitsBranchTmp->SetAddress(&fgDigitsArrTmp); return fgDigitsTreeTmp; } //____________________________________________________________________________ TTree *AliHLTEMCALUtils::GetClustersTree() { // // Create the clusters tree if necessary // if (fgClustersTreeTmp == NULL) { fgClustersTreeTmp = new TTree("EMCALclusters", "EMCAL clusters"); } return fgClustersTreeTmp; } //____________________________________________________________________________ void AliHLTEMCALUtils::ResetReconstructionTrees() { // // Reset the content of the reconstruction trees (digits and clusters) // TTree *tmp = 0; tmp = GetDigitsTree(); if (tmp) { fgDigitsArrTmp->Delete(); tmp->Reset(); } tmp = GetClustersTree(); if (tmp) { tmp->Reset(); } } //____________________________________________________________________________ void AliHLTEMCALUtils::DeleteReconstructionTrees() { // // Delete the reconstruction trees (digits and clusters) // TTree *tmp = 0; tmp = GetDigitsTree(); if (tmp) { fgDigitsArrTmp->Delete(); delete fgClustersTreeTmp; fgClustersTreeTmp = NULL; } tmp = GetClustersTree(); if (tmp) { delete fgDigitsTreeTmp; fgDigitsTreeTmp = 0; } } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::Raw2Clusters(AliRawReader* rawReader, TTree* clustersTree, Option_t* sDigitsOption) { // // Here we go from raw to clusters (not realeasing digits on the way) // This method combines ConvertDigits and Reconstruct of AliEMCALReconstructor // // Here we assume the raw reader is setup properly, so we do not need extra reset! // rawReader->Reset() ; Bool_t bReturn = kFALSE; if (!fgRawUtils) { // Try to initialize - this is potentially dangerous // We will use default geometry - most probably // The OCDB must exist! otherwise we will fail. if (GetRawUtils() == NULL) { // This is no go. AliErrorClass("Failed. No RawUtils."); return bReturn; } } if (!fgRecParam) { // Again we try to init the RecParams to default // One should better do it properly before! // The ocdb has to be there! for the recparams init if (GetRecParam() == NULL) { // This is no go. AliErrorClass("Failed. No RecParams."); return bReturn; } } if (GetDigitsTree() == NULL) { // This is no go. AliErrorClass("Failed. No Digits Tree."); return bReturn; } if (GetClusterizer() == NULL) { // This is no go. AliErrorClass("Failed. No clusterizer."); return bReturn; } // raw->digits fgRawUtils->SetOption(sDigitsOption); fgRawUtils->SetRawFormatHighLowGainFactor(fgRecParam->GetHighLowGainFactor()); fgRawUtils->SetRawFormatOrder(fgRecParam->GetOrderParameter()); fgRawUtils->SetRawFormatTau(fgRecParam->GetTau()); fgRawUtils->SetNoiseThreshold(fgRecParam->GetNoiseThreshold()); fgRawUtils->SetNPedSamples(fgRecParam->GetNPedSamples()); fgDigitsArrTmp->Delete(); fgRawUtils->Raw2Digits(rawReader, fgDigitsArrTmp); fgDigitsTreeTmp->Fill(); // digits->clusters fgClusterizer->SetOutput(clustersTree); if (fgDigitsArrTmp->GetEntries() > 0) { fgClusterizer->SetInput(fgDigitsTreeTmp); if (fgDebug > 0) { fgClusterizer->Digits2Clusters("deb all") ; } else { fgClusterizer->Digits2Clusters(""); } fgClusterizer->Clear(); } bReturn = kTRUE; return bReturn; } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::RawBuffer2Clusters(UChar_t *buffer, ULong_t buffSize, Int_t eqID, TTree* clustersTree, Option_t* sDigitsOption) { // // Method provided for convenience // Take the raw (DDL) mem buffer and put the clusters into the tree // Bool_t bReturn = kFALSE; AliRawReaderMemory rawReader; rawReader.Reset(); rawReader.SetEquipmentID(eqID); rawReader.SetMemory(buffer, buffSize); bReturn = AliHLTEMCALUtils::Raw2Clusters(&rawReader, clustersTree, sDigitsOption); return bReturn; } //____________________________________________________________________________ TTree* AliHLTEMCALUtils::RawBuffer2Clusters(UChar_t *buffer, ULong_t buffSize, Int_t eqID, Option_t* sDigitsOption) { // // Method provided for convenience // Take the raw (DDL) mem buffer and put the clusters into the tree // Bool_t bReturn = kFALSE; AliRawReaderMemory rawReader; rawReader.Reset(); rawReader.SetMemory(buffer, buffSize); rawReader.SetEquipmentID(eqID); TTree *clustersTree = GetClustersTree(); bReturn = AliHLTEMCALUtils::Raw2Clusters(&rawReader, clustersTree, sDigitsOption); if (bReturn == kFALSE) { clustersTree->Delete(); delete clustersTree; clustersTree = 0; } return clustersTree; } //____________________________________________________________________________ UChar_t *AliHLTEMCALUtils::FileToMemory(const char *fname, UChar_t *outpbuffer, ULong_t &buffSize) { // // Method provided for convenience // Take the raw (DDL) file and copy it to a mem buffer // buffSize = 0; outpbuffer = NULL; ifstream inputFile(fname, ifstream::in); if (!inputFile.good()) { AliErrorClass(Form("Unable to open file", fname)); return outpbuffer; } inputFile.seekg (0, ios::end); buffSize = inputFile.tellg(); inputFile.seekg (0, ios::beg); void *buff = calloc(buffSize, sizeof(char)); if (!buff) { AliErrorClass(Form("Unable to allocate memory: %d bytes", buffSize)); return outpbuffer; } inputFile.read((char*)buff, buffSize); outpbuffer = (UChar_t*)buff; AliInfoClass(Form("File read: %s. Allocated %d bytes at 0x%x", fname, buffSize, outpbuffer)); return outpbuffer; } //____________________________________________________________________________ Bool_t AliHLTEMCALUtils::FillESD(TTree* digitsTree, TTree* clustersTree, AliESDEvent* esd) { // // Fill the ESD just like the AliEMCALReconstructor // // This is offline code taken directly from the AliEMCALReconstructor // Factorization of the function in the AliEMCALReconstructor should allow calling the code directly // Trigger part on raw data has been removed. //######################################## //##############Fill CaloCells############ //######################################## if (digitsTree == 0) { AliErrorClass("Digits tree is NULL!"); return kFALSE; } if (clustersTree == 0) { AliErrorClass("Clusters tree is NULL!"); return kFALSE; } if (esd == 0) { AliErrorClass("Pointer to ESD not set!"); return kFALSE; } TClonesArray *digits = new TClonesArray("AliEMCALDigit",1000); TBranch *branchdig = digitsTree->GetBranch("EMCAL"); if (!branchdig) { AliErrorClass("Can't get the branch with the PHOS digits !"); return kFALSE; } branchdig->SetAddress(&digits); digitsTree->GetEvent(0); Int_t nDigits = digits->GetEntries(), idignew = 0 ; AliESDCaloCells &emcCells = *(esd->GetEMCALCells()); emcCells.CreateContainer(nDigits); emcCells.SetType(AliESDCaloCells::kEMCALCell); for (Int_t idig = 0 ; idig < nDigits ; idig++) { const AliEMCALDigit * dig = (const AliEMCALDigit*)digits->At(idig); if(dig->GetAmp() > 0 ) { emcCells.SetCell(idignew,dig->GetId(),dig->GetAmp(), dig->GetTime()); idignew++; } } emcCells.SetNumberOfCells(idignew); emcCells.Sort(); //------------------------------------------------------------ //-----------------CLUSTERS----------------------------- //------------------------------------------------------------ TObjArray *clusters = new TObjArray(100); TBranch *branch = clustersTree->GetBranch("EMCALECARP"); branch->SetAddress(&clusters); clustersTree->GetEvent(0); Int_t nClusters = clusters->GetEntries(), nClustersNew=0; esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters //###################################################### //#######################TRACK MATCHING############### //###################################################### //Fill list of integers, each one is index of track to which the cluster belongs. // step 1 - initialize array of matched track indexes Int_t *matchedTrack = new Int_t[nClusters]; for (Int_t iclus = 0; iclus < nClusters; iclus++) matchedTrack[iclus] = -1; // neg. index --> no matched track // step 2, change the flag for all matched clusters found in tracks Int_t iemcalMatch = -1; Int_t endtpc = esd->GetNumberOfTracks(); for (Int_t itrack = 0; itrack < endtpc; itrack++) { AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track iemcalMatch = track->GetEMCALcluster(); if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack; } //######################################## //##############Fill CaloClusters############# //######################################## esd->SetNumberOfEMCALClusters(nClusters); for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) { const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)clusters->At(iClust); //if(clust->GetClusterType()== AliESDCaloCluster::kEMCALClusterv1) nRP++; else nPC++; // Get information from EMCAL reconstruction points Float_t xyz[3]; TVector3 gpos; clust->GetGlobalPosition(gpos); for (Int_t ixyz=0; ixyz<3; ixyz++) xyz[ixyz] = gpos[ixyz]; Float_t elipAxis[2]; clust->GetElipsAxis(elipAxis); //Create digits lists Int_t cellMult = clust->GetMultiplicity(); //TArrayS digiList(digitMult); Float_t *amplFloat = clust->GetEnergiesList(); Int_t *digitInts = clust->GetAbsId(); TArrayS absIdList(cellMult); //Uncomment when unfolding is done //TArrayD fracList(cellMult); Int_t newCellMult = 0; for (Int_t iCell=0; iCell 0) { absIdList[newCellMult] = (UShort_t)(digitInts[iCell]); //Uncomment when unfolding is done //fracList[newCellMult] = amplFloat[iCell]/emcCells.GetCellAmplitude(digitInts[iCell]); newCellMult++; } } absIdList.Set(newCellMult); //Uncomment when unfolding is done //fracList.Set(newCellMult); if(newCellMult > 0) { // accept cluster if it has some digit nClustersNew++; //Primaries Int_t parentMult = 0; Int_t *parentList = clust->GetParents(parentMult); // fills the ESDCaloCluster AliESDCaloCluster * ec = new AliESDCaloCluster() ; ec->SetClusterType(AliESDCaloCluster::kEMCALClusterv1); ec->SetPosition(xyz); ec->SetE(clust->GetEnergy()); ec->SetNCells(newCellMult); //Change type of list from short to ushort UShort_t *newAbsIdList = new UShort_t[newCellMult]; //Uncomment when unfolding is done //Double_t *newFracList = new Double_t[newCellMult]; for(Int_t i = 0; i < newCellMult ; i++) { newAbsIdList[i]=absIdList[i]; //Uncomment when unfolding is done //newFracList[i]=fracList[i]; } ec->SetCellsAbsId(newAbsIdList); //Uncomment when unfolding is done //ec->SetCellsAmplitudeFraction(newFracList); ec->SetClusterDisp(clust->GetDispersion()); ec->SetClusterChi2(-1); //not yet implemented ec->SetM02(elipAxis[0]*elipAxis[0]) ; ec->SetM20(elipAxis[1]*elipAxis[1]) ; ec->SetTOF(clust->GetTime()) ; //time-of-fligh ec->SetNExMax(clust->GetNExMax()); //number of local maxima TArrayI arrayTrackMatched(1);// Only one track, temporal solution. arrayTrackMatched[0]= matchedTrack[iClust]; ec->AddTracksMatched(arrayTrackMatched); TArrayI arrayParents(parentMult,parentList); ec->AddLabels(arrayParents); // add the cluster to the esd object esd->AddCaloCluster(ec); delete ec; delete [] newAbsIdList ; //delete [] newFracList ; } } // cycle on clusters delete [] matchedTrack; esd->SetNumberOfEMCALClusters(nClustersNew); //if(nClustersNew != nClusters) //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew ); //Fill ESDCaloCluster with PID weights AliEMCALPID *pid = new AliEMCALPID; //pid->SetPrintInfo(kTRUE); pid->SetReconstructor(kTRUE); pid->RunPID(esd); delete pid; delete digits; delete clusters; // printf(" ## AliEMCALReconstructor::FillESD() is ended : ncl %i -> %i ### \n ",nClusters, nClustersNew); return kTRUE; }