X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSReconstructor.cxx;h=f8ac67d97dc88fe26e9918ecfe2386cb161cfe9d;hb=22fbde1655ee36b3a7898d46d89bb8da22d34243;hp=9c8b0d1df58a121fbe2da4596453eb49fb24f405;hpb=85c60a8ef6c7d60a77a305e6b551e6084ad1a9d8;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSReconstructor.cxx b/PHOS/AliPHOSReconstructor.cxx index 9c8b0d1df58..f8ac67d97dc 100644 --- a/PHOS/AliPHOSReconstructor.cxx +++ b/PHOS/AliPHOSReconstructor.cxx @@ -26,116 +26,390 @@ // --- Standard library --- // --- AliRoot header files --- -#include "AliESD.h" +#include "AliLog.h" +#include "AliAltroMapping.h" +#include "AliESDEvent.h" +#include "AliESDCaloCluster.h" #include "AliPHOSReconstructor.h" #include "AliPHOSClusterizerv1.h" #include "AliPHOSTrackSegmentMakerv1.h" #include "AliPHOSPIDv1.h" -#include "AliPHOSGetter.h" #include "AliPHOSTracker.h" #include "AliRawReader.h" +#include "AliPHOSTrigger.h" +#include "AliPHOSGeometry.h" +#include "AliPHOSRecoParam.h" +#include "AliPHOSRecoParamEmc.h" +#include "AliPHOSRecoParamCpv.h" +#include "AliPHOSDigit.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSRecParticle.h" +#include "AliPHOSRawDecoder.h" +#include "AliPHOSRawDecoderv1.h" +#include "AliPHOSRawDigiProducer.h" +#include "AliPHOSPulseGenerator.h" - ClassImp(AliPHOSReconstructor) Bool_t AliPHOSReconstructor::fgDebug = kFALSE ; +AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamEmc =0; // EMC rec. parameters +AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamCpv =0; // CPV rec. parameters //____________________________________________________________________________ - AliPHOSReconstructor::AliPHOSReconstructor() +AliPHOSReconstructor::AliPHOSReconstructor() : + fGeom(NULL) { // ctor -} + if (!fgkRecoParamEmc) { + AliWarning("The Reconstruction parameters for EMC nonitialized - Used default one"); + fgkRecoParamEmc = AliPHOSRecoParamEmc::GetEmcDefaultParameters(); + } + + if (!fgkRecoParamCpv) { + AliWarning("The Reconstruction parameters for CPV nonitialized - Used default one"); + fgkRecoParamCpv = AliPHOSRecoParamCpv::GetCpvDefaultParameters(); + } + + fGeom = AliPHOSGeometry::GetInstance("IHEP",""); +} //____________________________________________________________________________ AliPHOSReconstructor::~AliPHOSReconstructor() { // dtor - + delete fGeom; } //____________________________________________________________________________ -void AliPHOSReconstructor::Reconstruct(AliRunLoader* runLoader) const +void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const { - // method called by AliReconstruction; + // 'single-event' local reco method called by AliReconstruction; // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track - // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by + // segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by // the global tracking. - - TString headerFile(runLoader->GetFileName()) ; - TString branchName(runLoader->GetEventFolder()->GetName()) ; - - AliPHOSClusterizerv1 clu(headerFile, branchName); - clu.SetEventRange(0, -1) ; // do all the events + + AliPHOSClusterizerv1 clu(fGeom); + clu.SetInput(digitsTree); + clu.SetOutput(clustersTree); if ( Debug() ) - clu.ExecuteTask("deb all") ; + clu.Digits2Clusters("deb all") ; else - clu.ExecuteTask("") ; - + clu.Digits2Clusters("") ; } //____________________________________________________________________________ -void AliPHOSReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawreader) const +void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, + AliESDEvent* esd) const { - // method called by AliReconstruction; - // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track - // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by - // the global tracking. - // Here we reconstruct from Raw Data - - rawreader->Reset() ; - TString headerFile(runLoader->GetFileName()) ; - TString branchName(runLoader->GetEventFolder()->GetName()) ; + // This method produces PHOS rec-particles, + // then it creates AliESDtracks out of them and + // write tracks to the ESD + + AliPHOSTrackSegmentMaker *tsm = new AliPHOSTrackSegmentMakerv1(fGeom); + AliPHOSPID *pid = new AliPHOSPIDv1 (fGeom); + + // do current event; the loop over events is done by AliReconstruction::Run() + tsm->SetESD(esd) ; + tsm->SetInput(clustersTree); + if ( Debug() ) + tsm->Clusters2TrackSegments("deb all") ; + else + tsm->Clusters2TrackSegments("") ; - AliPHOSClusterizerv1 clu(headerFile, branchName); - clu.SetEventRange(0, -1) ; // do all the events - clu.SetRawReader(rawreader); + pid->SetInput(clustersTree, tsm->GetTrackSegments()) ; + pid->SetESD(esd) ; if ( Debug() ) - clu.ExecuteTask("deb all") ; + pid->TrackSegments2RecParticles("deb all") ; else - clu.ExecuteTask("") ; + pid->TrackSegments2RecParticles("") ; -} -//____________________________________________________________________________ -void AliPHOSReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const -{ // This function creates AliESDtracks from AliPHOSRecParticles // and // writes them to the ESD - Int_t eventNumber = runLoader->GetEventNumber() ; - - AliPHOSGetter::Instance()->Event(eventNumber, "P") ; - TClonesArray *recParticles = AliPHOSGetter::Instance()->RecParticles(); + TClonesArray *recParticles = pid->GetRecParticles(); Int_t nOfRecParticles = recParticles->GetEntries(); + esd->SetNumberOfPHOSClusters(nOfRecParticles) ; - esd->SetFirstPHOSCluster(esd->GetNumberOfTracks()) ; + esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ; + + AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption())); + + //#########Calculate trigger and set trigger info########### + + AliPHOSTrigger tr ; + // tr.SetPatchSize(1);//create 4x4 patches + tr.Trigger(); + + Float_t maxAmp2x2 = tr.Get2x2MaxAmplitude(); + Float_t maxAmpnxn = tr.GetnxnMaxAmplitude(); + Float_t ampOutOfPatch2x2 = tr.Get2x2AmpOutOfPatch() ; + Float_t ampOutOfPatchnxn = tr.GetnxnAmpOutOfPatch() ; + + Int_t iSM2x2 = tr.Get2x2SuperModule(); + Int_t iSMnxn = tr.GetnxnSuperModule(); + Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi(); + Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi(); + Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta(); + Int_t iCrystalEtanxn = tr.GetnxnCrystalEta(); + + AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d", maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2)); + AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d", maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn)); + + Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};// PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger. + Int_t iAbsId2x2 =-1; + Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};// PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger. + Int_t iAbsIdnxn =-1; + TVector3 pos2x2(-1,-1,-1); + TVector3 posnxn(-1,-1,-1); + fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2); + fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn); + fGeom->RelPosInAlice(iAbsId2x2, pos2x2); + fGeom->RelPosInAlice(iAbsIdnxn, posnxn); + + TArrayF triggerPosition(6); + triggerPosition[0] = pos2x2(0) ; + triggerPosition[1] = pos2x2(1) ; + triggerPosition[2] = pos2x2(2) ; + triggerPosition[3] = posnxn(0) ; + triggerPosition[4] = posnxn(1) ; + triggerPosition[5] = posnxn(2) ; + + TArrayF triggerAmplitudes(4); + triggerAmplitudes[0] = maxAmp2x2 ; + triggerAmplitudes[1] = ampOutOfPatch2x2 ; + triggerAmplitudes[2] = maxAmpnxn ; + triggerAmplitudes[3] = ampOutOfPatchnxn ; + + //esd->SetPHOSTriggerCells(triggerPosition); + esd->AddPHOSTriggerPosition(triggerPosition); + esd->AddPHOSTriggerAmplitudes(triggerAmplitudes); + + //###################################### + + // Read digits array + TBranch *branch = digitsTree->GetBranch("PHOS"); + if (!branch) { + AliError("can't get the branch with the PHOS digits !"); + return; + } + TClonesArray *fDigitsArr = new TClonesArray("AliPHOSDigit",100); + branch->SetAddress(&fDigitsArr); + branch->GetEntry(0); + + // Get the clusters array + TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP"); + if (!emcbranch) { + AliError("can't get the branch with the PHOS EMC clusters !"); + return; + } + + TObjArray *fEmcRecPoints = new TObjArray(100) ; + emcbranch->SetAddress(&fEmcRecPoints); + emcbranch->GetEntry(0); + + //Fill CaloClusters + const Float_t kBigShort = std::numeric_limits::max() - 1; + const Float_t nsec100 = 1e9*100.; // units of 0.01 ns + const Float_t gev500 = 500.; // units of GeV/500 for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) { AliPHOSRecParticle * rp = dynamic_cast(recParticles->At(recpart)); if (Debug()) rp->Print(); - AliESDCaloCluster * ec = new AliESDCaloCluster() ; -// AliESDtrack * et = new AliESDtrack() ; - // fills the ESDCaloCluster + // Get track segment and EMC rec.point associated with this rec.particle + AliPHOSTrackSegment *ts = static_cast(tsm->GetTrackSegments()->At(rp->GetPHOSTSIndex())); + + AliPHOSEmcRecPoint *emcRP = static_cast(fEmcRecPoints->At(ts->GetEmcIndex())); + AliESDCaloCluster *ec = new AliESDCaloCluster() ; + Float_t xyz[3]; for (Int_t ixyz=0; ixyz<3; ixyz++) xyz[ixyz] = rp->GetPos()[ixyz]; - ec->SetGlobalPosition(xyz); - ec->SetClusterEnergy(rp->Energy()); - ec->SetPid (rp->GetPID()) ; - ec->SetPrimaryIndex (rp->GetPrimaryIndex()); + + AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2])); + + //Create digits lists + Int_t digitMult = emcRP->GetDigitsMultiplicity(); + Int_t *digitsList = emcRP->GetDigitsList(); + Float_t *rpElist = emcRP->GetEnergiesList() ; + Short_t *amplList = new Short_t[digitMult]; + Short_t *timeList = new Short_t[digitMult]; + Short_t *digiList = new Short_t[digitMult]; + + + // Convert Float_t* and Int_t* to Short_t* to save memory + for (Int_t iDigit=0; iDigit(fDigitsArr->At(digitsList[iDigit])); + amplList[iDigit] = + (Short_t)(TMath::Min(rpElist[iDigit]*gev500,kBigShort)); // Energy in units of GeV/500 +// We should add here not full energy of digit, but unfolded one, stored in RecPoint +// amplList[iDigit] = +// (Short_t)(TMath::Min(digit->GetEnergy()*gev500,kBigShort)); // Energy in units of GeV/500 + timeList[iDigit] = + (Short_t)(TMath::Max(-kBigShort,TMath::Min(digit->GetTime()*nsec100,kBigShort))); // time in units of 0.01 ns + digiList[iDigit] = (Short_t)(digit->GetId()); + } + + + + //Primaries + Int_t primMult = 0; + Int_t *primList = emcRP->GetPrimaries(primMult); + + // fills the ESDCaloCluster + + ec->SetClusterType(AliESDCaloCluster::kPHOSCluster); + ec->SetPosition(xyz); //rec.point position in MARS + ec->SetE(rp->Energy()); //total particle energy + ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion + ec->SetPid(rp->GetPID()) ; //array of particle identification + ec->SetM02(emcRP->GetM2x()) ; //second moment M2x + ec->SetM20(emcRP->GetM2z()) ; //second moment M2z + ec->SetNExMax(emcRP->GetNExMax()); //number of local maxima + ec->SetEmcCpvDistance(ts->GetCpvDistance("r")); //Only radius, what about separate x,z???? + ec->SetClusterChi2(-1); //not yet implemented + ec->SetM11(-1) ; //not yet implemented + + //Digits Lists + TArrayS arrayAmpList(digitMult,amplList); + TArrayS arrayTimeList(digitMult,timeList); + TArrayS arrayIndexList(digitMult,digiList); + ec->AddDigitAmplitude(arrayAmpList); + ec->AddDigitTime(arrayTimeList); + ec->AddDigitIndex(arrayIndexList); + + //Distance to the nearest bad crystal + ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal()); + + //Array of MC indeces + TArrayI arrayPrim(primMult,primList); + ec->AddLabels(arrayPrim); + + //Array of tracks uncomment when available in future + //TArrayS arrayTrackMatched(1);// Only one track, temporal solution. + //arrayTrackMatched[0]= (Short_t)(matchedTrack[iClust]); + //ec->AddTracksMatched(arrayTrackMatched); + // add the track to the esd object esd->AddCaloCluster(ec); - delete ec; + delete ec; + delete [] amplList; + delete [] timeList; + delete [] digiList; } + fDigitsArr ->Delete(); + delete fDigitsArr; + fEmcRecPoints->Delete(); + delete fEmcRecPoints; + delete tsm; + delete pid; } -AliTracker* AliPHOSReconstructor::CreateTracker(AliRunLoader* runLoader) const +//____________________________________________________________________________ +AliTracker* AliPHOSReconstructor::CreateTracker() const { -// creates the PHOS tracker - if (!runLoader) return NULL; - return new AliPHOSTracker(runLoader); + // creates the PHOS tracker + return new AliPHOSTracker(); } +//____________________________________________________________________________ +void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const +{ + // Converts raw data to + // PHOS digits + // Works on a single-event basis + + rawReader->Reset() ; + + AliPHOSRawDecoder * dc ; + + const TObjArray* maps = AliPHOSRecoParamEmc::GetMappings(); + if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); + + AliAltroMapping *mapping[4]; + for(Int_t i = 0; i < 4; i++) { + mapping[i] = (AliAltroMapping*)maps->At(i); + } + + if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0) + dc=new AliPHOSRawDecoderv1(rawReader,mapping); + else + dc=new AliPHOSRawDecoder(rawReader,mapping); + + TString option = GetOption(); + if (option.Contains("OldRCUFormat")) + dc->SetOldRCUFormat(kTRUE); + else + dc->SetOldRCUFormat(kFALSE); + + dc->SubtractPedestals(fgkRecoParamEmc->SubtractPedestals()); + + TClonesArray *digits = new TClonesArray("AliPHOSDigit",1); + digits->SetName("DIGITS"); + Int_t bufsize = 32000; + digitsTree->Branch("PHOS", &digits, bufsize); + + AliPHOSRawDigiProducer pr; + pr.MakeDigits(digits,dc); + + delete dc ; + + //ADC counts -> GeV + if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0){ //"Energy" calculated as fit + for(Int_t i=0; iGetEntries(); i++) { + AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i); + digit->SetEnergy(digit->GetEnergy()*0.005); //We assume here 5 MeV/ADC channel + digit->SetTime(digit->GetTime()*1.e-7) ; //Here we assume sample step==100 ns TO BE FIXED!!!!!!!!!!!!! + } + } + else{ //Digits energy calculated as maximal energy + for(Int_t i=0; iGetEntries(); i++) { + AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i); + digit->SetEnergy(digit->GetEnergy()/AliPHOSPulseGenerator::GeV2ADC()); + } + } + + // Clean up digits below the noise threshold + // Assuming the digit noise to be 4 MeV, we suppress digits within + // 3-sigma of the noise. + // This parameter should be passed via AliPHOSRecoParamEmc later + + const Double_t emcDigitThreshold = 0.012; + for(Int_t i=0; iGetEntries(); i++) { + AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i); + if(digit->GetEnergy() < emcDigitThreshold) + digits->RemoveAt(i) ; + } + digits->Compress() ; + + //!!!!for debug!!! + Int_t modMax=-111; + Int_t colMax=-111; + Int_t rowMax=-111; + Float_t eMax=-333; + //!!!for debug!!! + + Int_t relId[4]; + for(Int_t iDigit=0; iDigitGetEntries(); iDigit++) { + AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit); + if(digit->GetEnergy()>eMax) { + fGeom->AbsToRelNumbering(digit->GetId(),relId); + eMax=digit->GetEnergy(); + modMax=relId[0]; + rowMax=relId[2]; + colMax=relId[3]; + } + } + + AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n\n", + modMax,colMax,rowMax,eMax)); + + digitsTree->Fill(); + digits->Delete(); + delete digits; +}