X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PHOS%2FAliPHOSReconstructor.cxx;h=3d928ee334485e0ed4e1273cd5d07bf81ec2518a;hb=532be2d489fe24d4eae739a099e7947365b38a81;hp=67e0ffdb450a976874ec8d06c8b90059d8e830b6;hpb=1d76a1f426f0e8ca4264ecf3f7ab5df7473048bd;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSReconstructor.cxx b/PHOS/AliPHOSReconstructor.cxx index 67e0ffdb450..3d928ee3344 100644 --- a/PHOS/AliPHOSReconstructor.cxx +++ b/PHOS/AliPHOSReconstructor.cxx @@ -52,8 +52,12 @@ #include "AliPHOSRawFitterv1.h" #include "AliPHOSRawFitterv2.h" #include "AliPHOSRawFitterv3.h" +#include "AliPHOSRawFitterv4.h" #include "AliPHOSRawDigiProducer.h" #include "AliPHOSPulseGenerator.h" +#include "AliPHOSTriggerRawDigit.h" +#include "AliPHOSTriggerRawDigiProducer.h" +#include "AliPHOSTriggerParameters.h" ClassImp(AliPHOSReconstructor) @@ -61,37 +65,42 @@ Bool_t AliPHOSReconstructor::fgDebug = kFALSE ; TClonesArray* AliPHOSReconstructor::fgDigitsArray = 0; // Array of PHOS digits TObjArray* AliPHOSReconstructor::fgEMCRecPoints = 0; // Array of EMC rec.points AliPHOSCalibData * AliPHOSReconstructor::fgCalibData = 0 ; - +TClonesArray* AliPHOSReconstructor::fgTriggerDigits = 0; // Array of PHOS trigger digits //____________________________________________________________________________ AliPHOSReconstructor::AliPHOSReconstructor() : - fGeom(NULL),fClusterizer(NULL),fTSM(NULL),fPID(NULL) + fGeom(NULL),fClusterizer(NULL),fTSM(NULL),fPID(NULL),fTmpDigLG(NULL) { // ctor - fGeom = AliPHOSGeometry::GetInstance("IHEP",""); - fClusterizer = new AliPHOSClusterizerv1 (fGeom); - fTSM = new AliPHOSTrackSegmentMakerv1(fGeom); - fPID = new AliPHOSPIDv1 (fGeom); - fgDigitsArray = new TClonesArray("AliPHOSDigit",100); - fgEMCRecPoints= new TObjArray(100) ; + fGeom = AliPHOSGeometry::GetInstance("IHEP",""); + fClusterizer = new AliPHOSClusterizerv1 (fGeom); + fTSM = new AliPHOSTrackSegmentMakerv1(fGeom); + fPID = new AliPHOSPIDv1 (fGeom); + fTmpDigLG = new TClonesArray("AliPHOSDigit",100); + fgDigitsArray = new TClonesArray("AliPHOSDigit",100); + fgEMCRecPoints = new TObjArray(100) ; if (!fgCalibData) fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number - + + fgTriggerDigits = new TClonesArray("AliPHOSTriggerRawDigit",100); + AliInfo(Form("PHOS bad channel map contains %d bad channel(s).\n", fgCalibData->GetNumOfEmcBadChannels())); } //____________________________________________________________________________ - AliPHOSReconstructor::~AliPHOSReconstructor() +AliPHOSReconstructor::~AliPHOSReconstructor() { // dtor delete fGeom; delete fClusterizer; delete fTSM; delete fPID; + delete fTmpDigLG; delete fgDigitsArray; delete fgEMCRecPoints; + delete fgTriggerDigits; } //____________________________________________________________________________ @@ -128,8 +137,6 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, else fTSM->Clusters2TrackSegments("") ; - fPID->SetEnergyCorrectionOn(GetRecoParam()->GetEMCEnergyCorrectionOn()); - fPID->SetInput(clustersTree, fTSM->GetTrackSegments()) ; fPID->SetESD(esd) ; if ( Debug() ) @@ -140,9 +147,6 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, TClonesArray *recParticles = fPID->GetRecParticles(); Int_t nOfRecParticles = recParticles->GetEntriesFast(); - esd->SetNumberOfPHOSClusters(nOfRecParticles) ; - esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ; - AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption())); // Read digits array @@ -165,60 +169,89 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, emcbranch->SetAddress(&fgEMCRecPoints); emcbranch->GetEntry(0); + + + // Trigger + + TBranch *tbranch = digitsTree->GetBranch("TPHOS"); + if (tbranch) { + + tbranch->SetAddress(&fgTriggerDigits); + tbranch->GetEntry(0); + + AliESDCaloTrigger* trgESD = esd->GetCaloTrigger("PHOS"); + + if (trgESD) { + trgESD->Allocate(fgTriggerDigits->GetEntriesFast()); + + for (Int_t i = 0; i < fgTriggerDigits->GetEntriesFast(); i++) { + AliPHOSTriggerRawDigit* tdig = (AliPHOSTriggerRawDigit*)fgTriggerDigits->At(i); + + Int_t mod,modX,modZ; + tdig->GetModXZ(mod,modX,modZ); + + const Int_t relId[4] = {5-mod,0,modX+1,modZ+1}; + Int_t absId; + + fGeom->RelToAbsNumbering(relId,absId); + trgESD->Add(mod,absId,tdig->GetAmp(),0.,(Int_t*)NULL,0,0,0); + } + } + } + +// //#########Calculate trigger and set trigger info########### - //#########Calculate trigger and set trigger info########### - - AliPHOSTrigger tr ; - // tr.SetPatchSize(1);//create 4x4 patches - tr.SetSimulation(kFALSE); - tr.Trigger(fgDigitsArray); +// AliPHOSTrigger tr ; +// // tr.SetPatchSize(1);//create 4x4 patches +// tr.SetSimulation(kFALSE); +// tr.Trigger(fgDigitsArray); - 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)); - - // Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger. - Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2}; - Int_t iAbsId2x2 =-1; - Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn}; - 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); +// 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)); + +// // Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger. +// Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2}; +// Int_t iAbsId2x2 =-1; +// Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn}; +// 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); //######################################## @@ -252,7 +285,7 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, //######################################## for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) { - AliPHOSRecParticle *rp = dynamic_cast(recParticles->At(recpart)); + AliPHOSRecParticle *rp = static_cast(recParticles->At(recpart)); if (Debug()) rp->Print(); // Get track segment and EMC rec.point associated with this rec.particle @@ -289,23 +322,27 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, Int_t primMult = 0; Int_t *primList = emcRP->GetPrimaries(primMult); - Float_t energy; + Float_t energy=0.; if (GetRecoParam()->EMCEcore2ESD()) energy = emcRP->GetCoreEnergy(); else energy = rp->Energy(); + //Apply nonlinearity correction + if(GetRecoParam()->GetEMCEnergyCorrectionOn()) + energy=CorrectNonlinearity(energy) ; // fills the ESDCaloCluster - ec->SetClusterType(AliESDCaloCluster::kPHOSCluster); + ec->SetType(AliVCluster::kPHOSNeutral); ec->SetPosition(xyz); //rec.point position in MARS ec->SetE(energy); //total or core particle energy - ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion - ec->SetPid(rp->GetPID()) ; //array of particle identification + ec->SetDispersion(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->SetTrackDistance(ts->GetCpvDistance("x"),ts->GetCpvDistance("z")); + ec->SetChi2(-1); //not yet implemented ec->SetTOF(emcRP->GetTime()); //Time of flight - already calibrated in EMCRecPoint //Cells contributing to clusters @@ -325,14 +362,28 @@ void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, arrayTrackMatched[0]= ts->GetTrackIndex(); ec->AddTracksMatched(arrayTrackMatched); - esd->AddCaloCluster(ec); + Int_t index = esd->AddCaloCluster(ec); + + //Set pointer to this cluster in ESD track + Int_t nt=esd->GetNumberOfTracks(); + for (Int_t itr=0; itrGetTrack(itr); + if(!esdTrack->IsPHOS()) + continue ; + if(esdTrack->GetPHOScluster()==-recpart){ //we store negative cluster number + esdTrack->SetPHOScluster(index) ; +//no garatie that only one track matched this cluster +// break ; + } + } + delete ec; delete [] fracList; delete [] absIdList; } - fgDigitsArray ->Delete(); - fgEMCRecPoints->Delete(); - recParticles ->Delete(); + fgDigitsArray ->Clear(); + fgEMCRecPoints->Clear("C"); + recParticles ->Clear(); //Store PHOS misalignment matrixes FillMisalMatrixes(esd) ; @@ -364,14 +415,16 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits mapping[i] = (AliAltroMapping*)maps->At(i); } - if (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0) + if (strcmp(GetRecoParam()->EMCFitterVersion(),"v0")==0) + fitter=new AliPHOSRawFitterv0(); + else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0) fitter=new AliPHOSRawFitterv1(); else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0) fitter=new AliPHOSRawFitterv2(); else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v3")==0) fitter=new AliPHOSRawFitterv3(); else - fitter=new AliPHOSRawFitterv0(); + fitter=new AliPHOSRawFitterv4(); fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals()); fitter->SetAmpOffset (GetRecoParam()->GetGlobalAltroOffset()); @@ -387,10 +440,22 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits rdp.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC rdp.SetCpvMinAmp(GetRecoParam()->GetCPVMinE()); rdp.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut()); - rdp.MakeDigits(digits,fitter); + rdp.MakeDigits(digits,fTmpDigLG,fitter); delete fitter ; + TClonesArray *tdigits = new TClonesArray("AliPHOSTriggerRawDigit",1); + tdigits->SetName("TDIGITS"); + digitsTree->Branch("TPHOS", &tdigits, bufsize); + + rawReader->Reset(); + AliPHOSTriggerRawDigiProducer tdp(rawReader); + + AliPHOSTriggerParameters* parameters = (AliPHOSTriggerParameters*)AliPHOSRecoParam::GetTriggerParameters(); + + tdp.SetTriggerParameters(parameters); + tdp.ProcessEvent(tdigits); + if (AliLog::GetGlobalDebugLevel() == 1) { Int_t modMax=-111; Int_t colMax=-111; @@ -413,10 +478,14 @@ 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; + + tdigits->Delete(); + delete tdigits; } //================================================================================== Float_t AliPHOSReconstructor::Calibrate(Float_t amp, Int_t absId)const{ @@ -478,7 +547,7 @@ void AliPHOSReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{ char path[255] ; TGeoHMatrix * m ; for(Int_t mod=0; mod<5; mod++){ - sprintf(path,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5 + snprintf(path,255,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5 if (gGeoManager->cd(path)){ m = gGeoManager->GetCurrentMatrix() ; esd->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ; @@ -489,5 +558,123 @@ void AliPHOSReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{ } } +//================================================================================== +Float_t AliPHOSReconstructor::CorrectNonlinearity(Float_t en){ + + //For backward compatibility, if no RecoParameters found + if(!GetRecoParam()){ + return 0.0241+1.0504*en+0.000249*en*en ; + } + + if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"NoCorrection")==0){ + return en ; + } + if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"Gustavo2005")==0){ + const Float_t *par=GetRecoParam()->GetNonlinearityParams() ; + return par[0]+par[1]*en + par[2]*en*en ; + } + if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"Henrik2010")==0){ + const Float_t *par=GetRecoParam()->GetNonlinearityParams() ; + return en*(par[0]+par[1]*TMath::Exp(-en*par[2]))*(1.+par[3]*TMath::Exp(-en*par[4]))*(1.+par[6]/(en*en+par[5])) ; + } + //For backward compatibility + if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"")==0){ + return 0.0241+1.0504*en+0.000249*en*en ; + } + return en ; +} + +void AliPHOSReconstructor::readTRUParameters(AliPHOSTriggerParameters* parameters) const +{ + //Read trigger parameters. + + TString path(gSystem->Getenv("ALICE_ROOT")); + path += "/PHOS/macros/Trigger/OCDB/"; + + for (Int_t mod = 2; mod < 5; ++mod) { // module + for (Int_t tru = 0; tru < 4; tru++) { // tru row + for (Int_t branch = 0; branch < 2; branch++) { // branch + + // Open the Appropriate pedestal file + TString fileName = path; + fileName += "pedestal_m"; + fileName = fileName += mod; + fileName+="_r"; + fileName+=tru; + fileName+="_b"; + fileName+=branch; + fileName+=".dat"; + std::ifstream instream; + instream.open(fileName.Data()); + + // Read pedestals from file + if( ! instream.is_open() ) + Printf("E-TRUPedestals: could not open %s", fileName.Data()); + else + { + Int_t ped[112]; + + char ch_s[36]; + char *ch_s_p = ch_s; + //Int_t nlines = 0; + + Int_t t_ped_0 =0; + Int_t t_ped_1 =0; + Int_t t_ped_2 =0; + + for(Int_t n=0; n<112; n++) + { + instream.getline(ch_s_p,36); + if (ch_s_p[23]<='9' && ch_s_p[23]>='0') + { + t_ped_0 = ch_s_p[23]-'0'; + } + else if (ch_s_p[23]>='A' && ch_s_p[23]<='Z') + { + t_ped_0 = ch_s_p[23]-'A'+10; + + } + + if (ch_s_p[22]<='9' && ch_s_p[22]>='0') + { + t_ped_1 = ch_s_p[22]-'0'; + } + else if (ch_s_p[22]<='Z' && ch_s_p[22]>='A') + { + t_ped_1 = ch_s_p[22]-'A'+10; + } + + if (ch_s_p[21]<='9' && ch_s_p[21]>='0') + { + t_ped_2 = ch_s_p[21]-'0'; + } + else if (ch_s_p[21]<='Z' && ch_s_p[21]>='A') + { + t_ped_2 = ch_s_p[21]-'A'+10; + } + + ped[n]=t_ped_2*256+t_ped_1*16+t_ped_0; + + + } + for (Int_t xrow = 0; xrow < 8; xrow++){ + for (Int_t zcol = 0; zcol < 14; zcol++){ + Int_t pedestal = ped[zcol*8+xrow]; + + if( pedestal < 612 && pedestal > 412 ) // resonable + parameters->SetTRUPedestal(pedestal, mod, tru, branch, xrow, zcol); + else // unresonable + continue; + } + } + } // Ends read of pedestals from branch from file. + instream.close(); + }// end branch + }// end tru + + }// end for mod +} + +