X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSReconstructor.cxx;h=1c5074a1209d9b0c0cb116fe8ef16d87f1c211d7;hb=2a4540e6ca39757f922322686fd718b882ef89cb;hp=6c6c764efbed75225ccb83cef5d9dedfbc922596;hpb=6e0b8dc1c50d591614bbebb383009dd4b4437382;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSReconstructor.cxx b/PHOS/AliPHOSReconstructor.cxx index 6c6c764efbe..1c5074a1209 100644 --- a/PHOS/AliPHOSReconstructor.cxx +++ b/PHOS/AliPHOSReconstructor.cxx @@ -16,8 +16,8 @@ /* $Id$ */ //_________________________________________________________________________ -//*-- -//*-- Yves Schutz (SUBATECH) +//-- +//-- Yves Schutz (SUBATECH) // Reconstruction class. Redesigned from the old AliReconstructionner class and // derived from STEER/AliReconstructor. // @@ -27,8 +27,10 @@ // --- AliRoot header files --- #include "AliLog.h" +#include "AliAltroMapping.h" #include "AliESDEvent.h" #include "AliESDCaloCluster.h" +#include "AliESDCaloCells.h" #include "AliPHOSReconstructor.h" #include "AliPHOSClusterizerv1.h" #include "AliPHOSTrackSegmentMakerv1.h" @@ -45,6 +47,8 @@ #include "AliPHOSEmcRecPoint.h" #include "AliPHOSRecParticle.h" #include "AliPHOSRawDecoder.h" +#include "AliPHOSRawDecoderv1.h" +#include "AliPHOSRawDecoderv2.h" #include "AliPHOSRawDigiProducer.h" #include "AliPHOSPulseGenerator.h" @@ -123,7 +127,7 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, else pid->TrackSegments2RecParticles("") ; - + // This function creates AliESDtracks from AliPHOSRecParticles // and // writes them to the ESD @@ -133,12 +137,34 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, esd->SetNumberOfPHOSClusters(nOfRecParticles) ; esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ; - + AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption())); + + // Read digits array + TBranch *branch = digitsTree->GetBranch("PHOS"); + if (!branch) { + AliError("can't get the branch with the PHOS digits !"); + return; + } + TClonesArray *digitsArray = new TClonesArray("AliPHOSDigit",100); + branch->SetAddress(&digitsArray); + 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 *emcRecPoints = new TObjArray(100) ; + emcbranch->SetAddress(&emcRecPoints); + emcbranch->GetEntry(0); //#########Calculate trigger and set trigger info########### - + AliPHOSTrigger tr ; // tr.SetPatchSize(1);//create 4x4 patches tr.Trigger(); @@ -155,8 +181,10 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 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)); + 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; @@ -187,42 +215,46 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 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; + //######################################## + //############# Fill CaloCells ########### + //######################################## + + Int_t nDigits = digitsArray->GetEntries(); + Int_t idignew = 0 ; + AliDebug(1,Form("%d digits",nDigits)); + + const Int_t knEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); + AliESDCaloCells &phsCells = *(esd->GetPHOSCells()); + phsCells.CreateContainer(nDigits); + phsCells.SetType(AliESDCaloCells::kPHOSCell); + + // Add to CaloCells only EMC digits with non-zero energy + for (Int_t idig = 0 ; idig < nDigits ; idig++) { + const AliPHOSDigit * dig = (const AliPHOSDigit*)digitsArray->At(idig); + if(dig->GetId() <= knEMC && dig->GetEnergy() > 0 ){ + //printf("i %d; id %d; amp %f; time %e\n", + //idignew,dig->GetId(),dig->GetEnergy(), dig->GetTime()); + phsCells.SetCell(idignew,dig->GetId(), dig->GetEnergy(), dig->GetTime()); + idignew++; + } } + phsCells.SetNumberOfCells(idignew); + phsCells.Sort(); - 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 + //######################################## + //############## Fill CaloClusters ####### + //######################################## for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) { - AliPHOSRecParticle * rp = dynamic_cast(recParticles->At(recpart)); + AliPHOSRecParticle *rp = dynamic_cast(recParticles->At(recpart)); if (Debug()) rp->Print(); // Get track segment and EMC rec.point associated with this rec.particle - AliPHOSTrackSegment *ts = static_cast(tsm->GetTrackSegments()->At(rp->GetPHOSTSIndex())); + AliPHOSTrackSegment *ts = static_cast(tsm->GetTrackSegments() + ->At(rp->GetPHOSTSIndex())); - AliPHOSEmcRecPoint *emcRP = static_cast(fEmcRecPoints->At(ts->GetEmcIndex())); + AliPHOSEmcRecPoint *emcRP = static_cast(emcRecPoints->At(ts->GetEmcIndex())); AliESDCaloCluster *ec = new AliESDCaloCluster() ; Float_t xyz[3]; @@ -230,34 +262,30 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, xyz[ixyz] = rp->GetPos()[ixyz]; 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(); - 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(digit->GetEnergy()*gev500,kBigShort)); // Energy in units of GeV/500 - timeList[iDigit] = - (Short_t)(TMath::Min(digit->GetTime()*nsec100,kBigShort)); // time in units of 0.01 ns - digiList[iDigit] = (Short_t)(digit->GetId()); + + // Create cell lists + + Int_t cellMult = emcRP->GetDigitsMultiplicity(); + Int_t *digitsList = emcRP->GetDigitsList(); + Float_t *rpElist = emcRP->GetEnergiesList() ; + UShort_t *absIdList = new UShort_t[cellMult]; + Double_t *fracList = new Double_t[cellMult]; + + for (Int_t iCell=0; iCell(digitsArray->At(digitsList[iCell])); + absIdList[iCell] = (UShort_t)(digit->GetId()); + if (digit->GetEnergy() > 0) + fracList[iCell] = rpElist[iCell]/digit->GetEnergy(); + else + fracList[iCell] = 0; } - + //Primaries Int_t primMult = 0; - Int_t *primInts = emcRP->GetPrimaries(primMult); - Short_t *primList = new Short_t[primMult]; - for (Int_t ipr=0; iprGetPrimaries(primMult); + // fills the ESDCaloCluster - - ec->SetPHOS(kTRUE); + 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 @@ -265,23 +293,20 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, ec->SetM02(emcRP->GetM2x()) ; //second moment M2x ec->SetM20(emcRP->GetM2z()) ; //second moment M2z ec->SetNExMax(emcRP->GetNExMax()); //number of local maxima - ec->SetEmcCpvDistance(-1); //not yet implemented + 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); + //Cells contributing to clusters + ec->SetNCells(cellMult); + ec->SetCellsAbsId(absIdList); + ec->SetCellsAmplitudeFraction(fracList); //Distance to the nearest bad crystal ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal()); //Array of MC indeces - TArrayS arrayPrim(primMult,primList); + TArrayI arrayPrim(primMult,primList); ec->AddLabels(arrayPrim); //Array of tracks uncomment when available in future @@ -290,17 +315,14 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, //ec->AddTracksMatched(arrayTrackMatched); // add the track to the esd object + esd->AddCaloCluster(ec); delete ec; - delete [] primList; - delete [] amplList; - delete [] timeList; - delete [] digiList; } - fDigitsArr ->Delete(); - delete fDigitsArr; - fEmcRecPoints->Delete(); - delete fEmcRecPoints; + digitsArray ->Delete(); + delete digitsArray; + emcRecPoints->Delete(); + delete emcRecPoints; delete tsm; delete pid; } @@ -312,6 +334,7 @@ AliTracker* AliPHOSReconstructor::CreateTracker() const return new AliPHOSTracker(); } +//____________________________________________________________________________ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const { // Converts raw data to @@ -320,30 +343,40 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits rawReader->Reset() ; - AliPHOSRawDecoder dc(rawReader); - TString option = GetOption(); - if (option.Contains("OldRCUFormat")) - dc.SetOldRCUFormat(kTRUE); - else - dc.SetOldRCUFormat(kFALSE); + 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 + if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v2")==0) + dc=new AliPHOSRawDecoderv2(rawReader,mapping); + else + dc=new AliPHOSRawDecoder(rawReader,mapping); + + dc->SetOldRCUFormat(fgkRecoParamEmc->IsOldRCUFormat()); - dc.SubtractPedestals(fgkRecoParamEmc->SubtractPedestals()); + 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); + AliPHOSRawDigiProducer pr(fgkRecoParamEmc,fgkRecoParamCpv); + pr.MakeDigits(digits,dc); + + delete dc ; - //ADC counts -> GeV - for(Int_t i=0; iGetEntries(); i++) { - AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i); - digit->SetEnergy(digit->GetEnergy()/AliPHOSPulseGenerator::GeV2ADC()); - } - //!!!!for debug!!! +/* Int_t modMax=-111; Int_t colMax=-111; Int_t rowMax=-111; @@ -364,7 +397,7 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits 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;