X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSTrigger.cxx;h=9cde799aecb98b40a03596c7292703a2b9465822;hb=6483babc2388acdec73dbabd46fa083362a36f58;hp=300799032be6d6f73994481802137e153a29bca6;hpb=03ecfe88df6e8f6f5d6d55ca60bc94ee3971b636;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSTrigger.cxx b/PHOS/AliPHOSTrigger.cxx index 300799032be..9cde799aecb 100644 --- a/PHOS/AliPHOSTrigger.cxx +++ b/PHOS/AliPHOSTrigger.cxx @@ -13,118 +13,166 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ -/* $Log $ */ //_________________________________________________________________________ // Class for trigger analysis. // Digits are grouped in TRU's (Trigger Units). A TRU consist of 16x28 // crystals ordered fNTRUPhi x fNTRUZ. The algorithm searches all possible -// 4x4 crystal combinations per each TRU, adding the digits amplitude and -// finding the maximum. Maximums are transformed in ADC time samples. -// Each time bin is compared to the trigger threshold until it is larger -// and then, triggers are set. Thresholds need to be fixed. +// 2x2 and nxn (n multiple of 2) crystal combinations per each TRU, adding the +// digits amplitude and finding the maximum. If found, look if it is isolated. +// Maxima are transformed in ADC time samples. Each time bin is compared to the trigger +// threshold until it is larger and then, triggers are set. Thresholds need to be fixed. // Usage: // // //Inside the event loop -// AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger +// AliPHOSTrigger *tr = new AliPHOSTrigger();//Init Trigger // tr->SetL0Threshold(100); // tr->SetL1JetLowPtThreshold(1000); +// tr->SetL1JetMediumPtThreshold(10000); // tr->SetL1JetHighPtThreshold(20000); +// .... // tr->Trigger(); //Execute Trigger -// tr->Print(""); //Print result, with "deb" option all data members -// //are printed +// tr->Print(""); //Print data members after calculation. +// // //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) ////////////////////////////////////////////////////////////////////////////// // --- ROOT system --- -//#include "TMatrixD.h" +#include "TMath.h" // --- ALIROOT system --- +#include "AliConfig.h" #include "AliPHOS.h" #include "AliPHOSTrigger.h" #include "AliPHOSGeometry.h" -#include "AliPHOSGetter.h" +#include "AliPHOSDigit.h" +#include "AliPHOSLoader.h" +#include "AliPHOSPulseGenerator.h" #include "AliTriggerInput.h" -//#include "AliLog.h" + ClassImp(AliPHOSTrigger) //______________________________________________________________________ AliPHOSTrigger::AliPHOSTrigger() : AliTriggerDetector(), - f2x2MaxAmp(-1), f2x2CrystalPhi(-1), f2x2CrystalEta(-1), - f4x4MaxAmp(-1), f4x4CrystalPhi(-1), f4x4CrystalEta(-1), - fL0Threshold(50), fL1JetLowPtThreshold(200), fL1JetHighPtThreshold(500), - fNTRU(8), fNTRUZ(2), fNTRUPhi(4), fSimulation(kTRUE) + f2x2MaxAmp(-1), f2x2CrystalPhi(-1), f2x2CrystalEta(-1), f2x2SM(0), + fnxnMaxAmp(-1), fnxnCrystalPhi(-1), fnxnCrystalEta(-1), fnxnSM(0), + fADCValuesHighnxn(0), fADCValuesLownxn(0), + fADCValuesHigh2x2(0), fADCValuesLow2x2(0), fDigitsList(0), + fAmptrus(0), fAmpmods(0), fTimeRtrus(0), + fL0Threshold(50), fL1JetLowPtThreshold(200), fL1JetMediumPtThreshold(500), + fL1JetHighPtThreshold(1000), + fNTRU(8), fNTRUZ(2), fNTRUPhi(4), + fNCrystalsPhi(16), + fNCrystalsZ(28), + fPatchSize(1), fIsolPatchSize(1), + f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1), + f2x2AmpOutOfPatchThres(2), fnxnAmpOutOfPatchThres(2), //2 GeV out of patch + fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE), + fSimulation(kTRUE), fIsolateInModule(kTRUE) { //ctor - - fADCValuesHigh4x4 = 0x0; //new Int_t[fTimeBins]; - fADCValuesLow4x4 = 0x0; //new Int_t[fTimeBins]; + fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins]; + fADCValuesLownxn = 0x0; //new Int_t[fTimeBins]; fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins]; fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins]; - fDigitsList = 0x0 ; - SetName("PHOS"); CreateInputs(); - //Print("") ; + fAmptrus = new TClonesArray("TMatrixD",1000); + fAmpmods = new TClonesArray("TMatrixD",1000); + fTimeRtrus = new TClonesArray("TMatrixD",1000); } - - //____________________________________________________________________________ -AliPHOSTrigger::AliPHOSTrigger(const AliPHOSTrigger & trig) - : AliTriggerDetector(trig) +AliPHOSTrigger::AliPHOSTrigger(const AliPHOSTrigger & trig) : + AliTriggerDetector(trig), + f2x2MaxAmp(trig.f2x2MaxAmp), + f2x2CrystalPhi(trig.f2x2CrystalPhi), + f2x2CrystalEta(trig.f2x2CrystalEta), + f2x2SM(trig.f2x2SM), + fnxnMaxAmp(trig.fnxnMaxAmp), + fnxnCrystalPhi(trig.fnxnCrystalPhi), + fnxnCrystalEta(trig.fnxnCrystalEta), + fnxnSM(trig.fnxnSM), + fADCValuesHighnxn(trig.fADCValuesHighnxn), + fADCValuesLownxn(trig.fADCValuesLownxn), + fADCValuesHigh2x2(trig.fADCValuesHigh2x2), + fADCValuesLow2x2(trig.fADCValuesLow2x2), + fDigitsList(trig.fDigitsList), + fAmptrus(trig.fAmptrus), fAmpmods(trig.fAmpmods), fTimeRtrus(trig.fTimeRtrus), + fL0Threshold(trig.fL0Threshold), + fL1JetLowPtThreshold(trig.fL1JetLowPtThreshold), + fL1JetMediumPtThreshold(trig.fL1JetMediumPtThreshold), + fL1JetHighPtThreshold(trig.fL1JetHighPtThreshold), + fNTRU(trig.fNTRU), + fNTRUZ(trig.fNTRUZ), + fNTRUPhi(trig.fNTRUPhi), + fNCrystalsPhi(trig.fNCrystalsPhi), + fNCrystalsZ(trig. fNCrystalsZ), + fPatchSize(trig.fPatchSize), + fIsolPatchSize(trig.fIsolPatchSize), + f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch), + fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch), + f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres), + fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres), + fIs2x2Isol(trig.fIs2x2Isol), + fIsnxnIsol(trig.fIsnxnIsol), + fSimulation(trig.fSimulation), + fIsolateInModule(trig.fIsolateInModule) { - // cpy ctor +} - f2x2MaxAmp = trig.f2x2MaxAmp ; - f4x4MaxAmp = trig.f4x4MaxAmp ; - f2x2CrystalPhi = trig.f2x2CrystalPhi ; - f4x4CrystalPhi = trig.f4x4CrystalPhi ; - f2x2CrystalEta = trig.f2x2CrystalEta ; - f4x4CrystalEta = trig.f4x4CrystalEta ; - fADCValuesHigh4x4 = trig.fADCValuesHigh4x4 ; - fADCValuesLow4x4 = trig.fADCValuesLow4x4 ; - fADCValuesHigh2x2 = trig.fADCValuesHigh2x2 ; - fADCValuesLow2x2 = trig.fADCValuesLow2x2 ; - fDigitsList = trig.fDigitsList ; - fL0Threshold = trig.fL0Threshold ; - fL1JetLowPtThreshold = trig.fL1JetLowPtThreshold ; - fL1JetHighPtThreshold = trig.fL1JetHighPtThreshold ; - fNTRU = trig.fNTRU ; - fNTRUZ = trig.fNTRUZ ; - fNTRUPhi = trig.fNTRUPhi ; - fSimulation = trig.fSimulation ; - +//_________________________________________________________________________ +AliPHOSTrigger::~AliPHOSTrigger() +{ + // dtor + + if(fADCValuesHighnxn)delete [] fADCValuesHighnxn; + if(fADCValuesLownxn)delete [] fADCValuesLownxn; + if(fADCValuesHigh2x2)delete [] fADCValuesHigh2x2; + if(fADCValuesLow2x2)delete [] fADCValuesLow2x2; + // fDigitsList is now ours... + if(fAmptrus) { fAmptrus->Delete() ; delete fAmptrus ; } + if(fAmpmods) { fAmpmods->Delete() ; delete fAmpmods ; } + if(fTimeRtrus) { fTimeRtrus->Delete(); delete fTimeRtrus; } } //_________________________________________________________________________ +AliPHOSTrigger & AliPHOSTrigger::operator = (const AliPHOSTrigger &) +{ + Fatal("operator =", "no implemented"); + return *this; +} + void AliPHOSTrigger::CreateInputs() { // inputs // Do not create inputs again!! if( fInputs.GetEntriesFast() > 0 ) return; + + TString name = GetName(); - fInputs.AddLast( new AliTriggerInput( "PHOS_L0", "PHOS L0", 0x02 ) ); - fInputs.AddLast( new AliTriggerInput( "PHOS_JetHPt_L1","PHOS Jet High Pt L1", 0x04 ) ); - fInputs.AddLast( new AliTriggerInput( "PHOS_JetLPt_L1","PHOS Jet Low Pt L1", 0x08 ) ); + fInputs.AddLast( new AliTriggerInput( "PHOS_L0", name, 0 ) ); + fInputs.AddLast( new AliTriggerInput( "PHOS_JetHPt_L1",name, 1 ) ); + fInputs.AddLast( new AliTriggerInput( "PHOS_JetMPt_L1",name, 1 ) ); + fInputs.AddLast( new AliTriggerInput( "PHOS_JetLPt_L1",name, 1 ) ); } //____________________________________________________________________________ -void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry * geom, TClonesArray * ampmatrix, TClonesArray * timeRmatrix) const { +void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry * geom) const { //Orders digits ampitudes list and times in fNTRU TRUs (28x16 crystals) //per module. Each TRU is a TMatrixD, and they are kept in TClonesArrays. //In a module, the number of TRU in phi is fNTRUPhi, and the number of - //TRU in eta is fNTRUZ. + //TRU in eta is fNTRUZ. Also fill a matrix with all amplitudes in module for isolation studies. //Check data members @@ -133,8 +181,6 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry //Initilize and declare variables Int_t nModules = geom->GetNModules(); - Int_t nCrystalsPhi = geom->GetNPhi()/fNTRUPhi ;// 64/4=16 - Int_t nCrystalsZ = geom->GetNZ()/fNTRUZ ;// 56/2=28 Int_t relid[4] ; Float_t amp = -1; Float_t timeR = -1; @@ -142,26 +188,39 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry //List of TRU matrices initialized to 0. for(Int_t k = 0; k < fNTRU*nModules ; k++){ - TMatrixD * amptrus = new TMatrixD(nCrystalsPhi,nCrystalsZ) ; - TMatrixD * timeRtrus = new TMatrixD(nCrystalsPhi,nCrystalsZ) ; - for(Int_t i = 0; i < nCrystalsPhi; i++){ - for(Int_t j = 0; j < nCrystalsZ; j++){ - (*amptrus)(i,j) = 0.0; - (*timeRtrus)(i,j) = 0.0; + TMatrixD amptrus(fNCrystalsPhi,fNCrystalsZ) ; + TMatrixD timeRtrus(fNCrystalsPhi,fNCrystalsZ) ; + for(Int_t i = 0; i < fNCrystalsPhi; i++){ + for(Int_t j = 0; j < fNCrystalsZ; j++){ + amptrus(i,j) = 0.0; + timeRtrus(i,j) = 0.0; + } + } + new((*fAmptrus)[k]) TMatrixD(amptrus) ; + new((*fTimeRtrus)[k]) TMatrixD(timeRtrus) ; + } + + //List of Modules matrices initialized to 0. + Int_t nmodphi = geom->GetNPhi(); + Int_t nmodz = geom->GetNZ(); + + for(Int_t k = 0; k < nModules ; k++){ + TMatrixD ampmods(nmodphi,nmodz) ; + for(Int_t i = 0; i < nmodphi; i++){ + for(Int_t j = 0; j < nmodz; j++){ + ampmods(i,j) = 0.0; } } - new((*ampmatrix)[k]) TMatrixD(*amptrus) ; - new((*timeRmatrix)[k]) TMatrixD(*timeRtrus) ; + new((*fAmpmods)[k]) TMatrixD(ampmods) ; } AliPHOSDigit * dig ; //Digits loop to fill TRU matrices with amplitudes. - for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){ dig = static_cast(digits->At(idig)) ; - amp = dig->GetAmp() ; // Energy of the digit (arbitrary units) + amp = dig->GetEnergy() ; // Energy of the digit id = dig->GetId() ; // Id label of the cell timeR = dig->GetTimeR() ; // Earliest time of the digit geom->AbsToRelNumbering(id, relid) ; @@ -172,36 +231,40 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry //relid[3] = column <= 56 (fNZ) if(relid[1] == 0){//Not CPV, Only EMC digits - + //############# TRU ################### //Check to which TRU in the supermodule belongs the crystal. //Supermodules are divided in a TRU matrix of dimension //(fNTRUPhi,fNTRUZ). - //Each TRU is a crystal matrix of dimension (nCrystalsPhi,nCrystalsZ) + //Each TRU is a crystal matrix of dimension (fNCrystalsPhi,fNCrystalsZ) //First calculate the row and column in the supermodule //of the TRU to which the crystal belongs. - Int_t col = (relid[3]-1)/nCrystalsZ+1; - Int_t row = (relid[2]-1)/nCrystalsPhi+1; + Int_t col = (relid[3]-1)/fNCrystalsZ+1; + Int_t row = (relid[2]-1)/fNCrystalsPhi+1; //Calculate label number of the TRU Int_t itru = (row-1) + (col-1)*fNTRUPhi + (relid[0]-1)*fNTRU ; //Fill TRU matrix with crystal values - TMatrixD * amptrus = dynamic_cast(ampmatrix->At(itru)) ; - TMatrixD * timeRtrus = dynamic_cast(timeRmatrix->At(itru)) ; + TMatrixD * amptrus = dynamic_cast(fAmptrus ->At(itru)) ; + TMatrixD * timeRtrus = dynamic_cast(fTimeRtrus->At(itru)) ; //Calculate row and column of the crystal inside the TRU with number itru - Int_t irow = (relid[2]-1) - (row-1) * nCrystalsPhi; - Int_t icol = (relid[3]-1) - (col-1) * nCrystalsZ; + Int_t irow = (relid[2]-1) - (row-1) * fNCrystalsPhi; + Int_t icol = (relid[3]-1) - (col-1) * fNCrystalsZ; (*amptrus)(irow,icol) = amp ; (*timeRtrus)(irow,icol) = timeR ; + + //####################MODULE MATRIX ################## + TMatrixD * ampmods = dynamic_cast(fAmpmods->At(relid[0]-1)) ; + (*ampmods)(relid[2]-1,relid[3]-1) = amp ; } } } //______________________________________________________________________ -void AliPHOSTrigger::GetCrystalPhiEtaIndexInModuleFromTRUIndex(const Int_t itru,const Int_t iphitru,const Int_t ietatru,Int_t &iphiMod,Int_t &ietaMod,const AliPHOSGeometry* geom) const +void AliPHOSTrigger::GetCrystalPhiEtaIndexInModuleFromTRUIndex(const Int_t itru,const Int_t iphitru,const Int_t ietatru,Int_t &iphiMod,Int_t &ietaMod) const { // This method transforms the (eta,phi) index of a crystals in a // TRU matrix into Super Module (eta,phi) index. @@ -212,105 +275,200 @@ void AliPHOSTrigger::GetCrystalPhiEtaIndexInModuleFromTRUIndex(const Int_t itru, Int_t row = itru - (col-1)*fNTRUPhi + 1; //Calculate the (eta,phi) index in SM - Int_t nCrystalsPhi = geom->GetNPhi()/fNTRUPhi; - Int_t nCrystalsZ = geom->GetNZ()/fNTRUZ; - iphiMod = nCrystalsPhi*(row-1) + iphitru + 1 ; - ietaMod = nCrystalsZ*(col-1) + ietatru + 1 ; + iphiMod = fNCrystalsPhi*(row-1) + iphitru + 1 ; + ietaMod = fNCrystalsZ*(col-1) + ietatru + 1 ; } + //____________________________________________________________________________ -void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t imod, TMatrixD *ampmax2, TMatrixD *ampmax4, const AliPHOSGeometry *geom){ - //Sums energy of all possible 2x2 (L0) and 4x4 (L1) crystals per each TRU. +Bool_t AliPHOSTrigger::IsPatchIsolated(Int_t iPatchType, const Int_t imod, const Int_t mtru, const Float_t maxamp, const Int_t maxphi, const Int_t maxeta) { + + //Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch, + //inside isolation patch . iPatchType = 0 means calculation for 2x2 patch, + //iPatchType = 1 means calculation for nxn patch. + //In the next table there is an example of the different options of patch size and isolation patch size: + // Patch Size (fPatchSize) + // 0 1 2 + // fIsolPatchSize 2x2 (not overlap) 4x4 (overlapped) 6x6(overlapped) ... + // 1 4x4 8x8 10x10 + // 2 6x6 12x12 14x14 + // 3 8x8 16x16 18x18 + + Bool_t b = kFALSE; + Float_t amp = 0; + + //Get matrix of TRU or Module with maximum amplitude patch. + Int_t itru = mtru+imod*fNTRU ; //number of tru, min 0 max 8*5. + TMatrixD * ampmatrix = 0x0; + Int_t colborder = 0; + Int_t rowborder = 0; + + if(fIsolateInModule){ + ampmatrix = dynamic_cast(fAmpmods->At(imod)) ; + rowborder = fNCrystalsPhi*fNTRUPhi; + colborder = fNCrystalsZ*fNTRUZ; + AliDebug(2,"Isolate trigger in Module"); + } + else{ + ampmatrix = dynamic_cast(fAmptrus->At(itru)) ; + rowborder = fNCrystalsPhi; + colborder = fNCrystalsZ; + AliDebug(2,"Isolate trigger in TRU"); + } + + //Define patch cells + Int_t isolcells = fIsolPatchSize*(1+iPatchType); + Int_t ipatchcells = 2*(1+fPatchSize*iPatchType); + Int_t minrow = maxphi - isolcells; + Int_t mincol = maxeta - isolcells; + Int_t maxrow = maxphi + isolcells + ipatchcells; + Int_t maxcol = maxeta + isolcells + ipatchcells; + + AliDebug(2,Form("Number of added Isol Cells %d, Patch Size %d",isolcells, ipatchcells)); + AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol)); + + if(minrow < 0 || mincol < 0 || maxrow > rowborder || maxcol > colborder){ + AliDebug(1,Form("Out of Module/TRU range, cannot isolate patch")); + return kFALSE; + } + + //Add amplitudes in all isolation patch + for(Int_t irow = minrow ; irow < maxrow; irow ++) + for(Int_t icol = mincol ; icol < maxcol ; icol ++) + amp += (*ampmatrix)(irow,icol); + + AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp)); + + if(TMath::Nint(amp*1E5) < TMath::Nint(maxamp*1E5)){ + AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp)); + return kFALSE; + } + else + amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch. + + AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp)); + + //Fill isolation amplitude data member and say if patch is isolated. + if(iPatchType == 0){ //2x2 case + f2x2AmpOutOfPatch = amp; + if(amp < f2x2AmpOutOfPatchThres) + b=kTRUE; + } + else if(iPatchType == 1){ //nxn case + fnxnAmpOutOfPatch = amp; + if(amp < fnxnAmpOutOfPatchThres) + b=kTRUE; + } + + return b; + +} + + +//____________________________________________________________________________ +void AliPHOSTrigger::MakeSlidingCell(const Int_t imod, TMatrixD &max2, TMatrixD &maxn) +{ + //Sums energy of all possible 2x2 (L0) and nxn (L1) crystals per each TRU. //Fast signal in the experiment is given by 2x2 crystals, //for this reason we loop inside the TRU crystals by 2. //Declare and initialize varibles - Int_t nCrystalsPhi = geom->GetNPhi()/fNTRUPhi ;// 64/4=16 - Int_t nCrystalsZ = geom->GetNZ()/fNTRUZ ;// 56/2=28 Float_t amp2 = 0 ; - Float_t amp4 = 0 ; - for(Int_t i = 0; i < 3; i++){ + Float_t ampn = 0 ; + for(Int_t i = 0; i < 4; i++){ for(Int_t j = 0; j < fNTRU; j++){ - (*ampmax2)(i,j) = -1; - (*ampmax4)(i,j) = -1; + ampmax2(i,j) = -1; + ampmaxn(i,j) = -1; } } //Create matrix that will contain 2x2 amplitude sums - //used to calculate the 4x4 sums - TMatrixD * tru2x2 = new TMatrixD(nCrystalsPhi/2,nCrystalsZ/2) ; - - for(Int_t i = 0; i < nCrystalsPhi/2; i++) - for(Int_t j = 0; j < nCrystalsZ/2; j++) - (*tru2x2)(i,j) = 0.0; + //used to calculate the nxn sums + TMatrixD tru2x2(fNCrystalsPhi/2,fNCrystalsZ/2) ; + for(Int_t i = 0; i < fNCrystalsPhi/2; i++) + for(Int_t j = 0; j < fNCrystalsZ/2; j++) + tru2x2(i,j) = -1.; //Loop over all TRUS in a module - for(Int_t itru = 0 + (imod - 1) * fNTRU ; itru < imod*fNTRU ; itru++){ - TMatrixD * amptru = dynamic_cast(amptrus->At(itru)) ; - TMatrixD * timeRtru = dynamic_cast(timeRtrus->At(itru)) ; - Int_t mtru = itru-(imod-1)*fNTRU ; //Number of TRU in Module + for(Int_t itru = 0 + imod * fNTRU ; itru < (imod+1)*fNTRU ; itru++){ + TMatrixD * amptru = dynamic_cast(fAmptrus ->At(itru)) ; + TMatrixD * timeRtru = dynamic_cast(fTimeRtrus->At(itru)) ; + Int_t mtru = itru-imod*fNTRU ; //Number of TRU in Module //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP) - for(Int_t irow = 0 ; irow < nCrystalsPhi; irow += 2){ - for(Int_t icol = 0 ; icol < nCrystalsZ ; icol += 2){ + for(Int_t irow = 0 ; irow < fNCrystalsPhi; irow += 2){ + for(Int_t icol = 0 ; icol < fNCrystalsZ ; icol += 2){ amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+ (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1); - //Fill new matrix with added 2x2 crystals for use in 4x4 sums - (*tru2x2)(irow/2,icol/2) = amp2 ; + //Fill new matrix with added 2x2 crystals for use in nxn sums + tru2x2(irow/2,icol/2) = amp2 ; //Select 2x2 maximum sums to select L0 - if(amp2 > (*ampmax2)(0,mtru)){ - (*ampmax2)(0,mtru) = amp2 ; - (*ampmax2)(1,mtru) = irow; - (*ampmax2)(2,mtru) = icol; + if(amp2 > ampmax2(0,mtru)){ + ampmax2(0,mtru) = amp2 ; + ampmax2(1,mtru) = irow; + ampmax2(2,mtru) = icol; } } } //Find most recent time in the selected 2x2 cell - (*ampmax2)(3,mtru) = 1 ; - Int_t row2 = static_cast ((*ampmax2)(1,mtru)); - Int_t col2 = static_cast ((*ampmax2)(2,mtru)); + ampmax2(3,mtru) = 1 ; + Int_t row2 = static_cast (ampmax2(1,mtru)); + Int_t col2 = static_cast (ampmax2(2,mtru)); for(Int_t i = 0; i<2; i++){ for(Int_t j = 0; j<2; j++){ if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){ - if((*timeRtru)(row2+i,col2+j) < (*ampmax2)(3,mtru) ) - (*ampmax2)(3,mtru) = (*timeRtru)(row2+i,col2+j); + if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) ) + ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); } } } - //Sliding 4x4, add 4x4 amplitudes (OVERLAP) - for(Int_t irow = 0 ; irow < nCrystalsPhi/2; irow++){ - for(Int_t icol = 0 ; icol < nCrystalsZ/2 ; icol++){ - if( (irow+1) < nCrystalsPhi/2 && (icol+1) < nCrystalsZ/2){//Avoid exit the TRU - amp4 = (*tru2x2)(irow,icol)+(*tru2x2)(irow+1,icol)+ - (*tru2x2)(irow,icol+1)+(*tru2x2)(irow+1,icol+1); - //Select 4x4 maximum sums to select L1 - if(amp4 > (*ampmax4)(0,mtru)){ - (*ampmax4)(0,mtru) = amp4 ; - (*ampmax4)(1,mtru) = irow*2; - (*ampmax4)(2,mtru) = icol*2; + //Sliding nxn, add nxn amplitudes (OVERLAP) + if(fPatchSize > 0){ + for(Int_t irow = 0 ; irow < fNCrystalsPhi/2; irow++){ + for(Int_t icol = 0 ; icol < fNCrystalsZ/2 ; icol++){ + ampn = 0; + if( (irow+fPatchSize) < fNCrystalsPhi/2 && (icol+fPatchSize) < fNCrystalsZ/2){//Avoid exit the TRU + for(Int_t i = 0 ; i <= fPatchSize ; i++) + for(Int_t j = 0 ; j <= fPatchSize ; j++) + ampn += tru2x2(irow+i,icol+j); + //Select nxn maximum sums to select L1 + if(ampn > ampmaxn(0,mtru)){ + ampmaxn(0,mtru) = ampn ; + ampmaxn(1,mtru) = irow*2; + ampmaxn(2,mtru) = icol*2; + } } } } - } - - //Find most recent time in selected 4x4 cell - (*ampmax4)(3,mtru) = 1 ; - Int_t row4 = static_cast ((*ampmax4)(1,mtru)); - Int_t col4 = static_cast ((*ampmax4)(2,mtru)); - for(Int_t i = 0; i<4; i++){ - for(Int_t j = 0; j<4; j++){ - if((*amptru)(row4+i,col4+j) > 0 && (*timeRtru)(row4+i,col4+j)> 0){ - if((*timeRtru)(row4+i,col4+j) < (*ampmax4)(3,mtru) ) - (*ampmax4)(3,mtru) = (*timeRtru)(row4+i,col4+j); + + //Find most recent time in selected nxn cell + ampmaxn(3,mtru) = 1 ; + Int_t rown = static_cast (ampmaxn(1,mtru)); + Int_t coln = static_cast (ampmaxn(2,mtru)); + for(Int_t i = 0; i<4*fPatchSize; i++){ + for(Int_t j = 0; j<4*fPatchSize; j++){ + if( (rown+i) < fNCrystalsPhi && (coln+j) < fNCrystalsZ/2){//Avoid exit the TRU + if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){ + if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) ) + ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); + } + } } } } + else { + ampmaxn(0,mtru) = ampmax2(0,mtru); + ampmaxn(1,mtru) = ampmax2(1,mtru); + ampmaxn(2,mtru) = ampmax2(2,mtru); + ampmaxn(3,mtru) = ampmax2(3,mtru); + } } } + //____________________________________________________________________________ void AliPHOSTrigger::Print(const Option_t * opt) const { @@ -325,9 +483,20 @@ void AliPHOSTrigger::Print(const Option_t * opt) const printf( " -2x2 crystals sum (not overlapped): %10.2f, in Super Module %d\n", f2x2MaxAmp,f2x2SM) ; printf( " -2x2 from row %d to row %d and from column %d to column %d\n", f2x2CrystalPhi, f2x2CrystalPhi+2, f2x2CrystalEta, f2x2CrystalEta+2) ; - printf( " -4x4 crystals sum (overlapped) : %10.2f, in Super Module %d\n", - f4x4MaxAmp,f4x4SM) ; - printf( " -4x4 from row %d to row %d and from column %d to column %d\n", f4x4CrystalPhi, f4x4CrystalPhi+4, f4x4CrystalEta, f4x4CrystalEta+4) ; + printf( " -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n", + 2*fIsolPatchSize+2, 2*fIsolPatchSize+2, f2x2AmpOutOfPatch, f2x2AmpOutOfPatchThres,static_cast (fIs2x2Isol)) ; + if(fPatchSize > 0){ + printf( " Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)); + printf( " -nxn crystals sum (overlapped) : %10.2f, in Super Module %d\n", + fnxnMaxAmp,fnxnSM) ; + printf( " -nxn from row %d to row %d and from column %d to column %d\n", fnxnCrystalPhi, fnxnCrystalPhi+4*fPatchSize, fnxnCrystalEta, fnxnCrystalEta+4*fPatchSize) ; + printf( " -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n", + 4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) , fnxnAmpOutOfPatch, fnxnAmpOutOfPatchThres,static_cast (fIsnxnIsol) ) ; + } + + printf( " Isolate in Module? %d\n", + fIsolateInModule) ; + printf( " Threshold for LO %10.1f\n", fL0Threshold) ; @@ -341,157 +510,211 @@ void AliPHOSTrigger::Print(const Option_t * opt) const if(in->GetValue()) printf( " *** PHOS Jet Low Pt for L1 is set ***\n") ; + printf( " Jet Medium Pt Threshold for L1 %10.2f\n", fL1JetMediumPtThreshold) ; + in = (AliTriggerInput*)fInputs.FindObject( "PHOS_JetMPt_L1" ); + if(in->GetValue()) + printf( " *** PHOS Jet Medium Pt for L1 is set ***\n") ; + printf( " Jet High Pt Threshold for L1 %10.2f\n", fL1JetHighPtThreshold) ; in = (AliTriggerInput*) fInputs.FindObject( "PHOS_JetHPt_L1" ); if(in->GetValue()) printf( " *** PHOS Jet High Pt for L1 is set ***\n") ; - + } //____________________________________________________________________________ -void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD * ampmax2, const TMatrixD * ampmax4, const AliPHOSGeometry *geom) +void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD & ampmax2, const TMatrixD & ampmaxn) { - //Checks the 2x2 and 4x4 maximum amplitude per each TRU and compares - //with the different L0 and L1 triggers thresholds + //Checks the 2x2 and nxn maximum amplitude per each TRU and compares + //with the different L0 and L1 triggers thresholds. It finds if maximum amplitudes are isolated. //Initialize variables Float_t max2[] = {-1,-1,-1,-1} ; - Float_t max4[] = {-1,-1,-1,-1} ; - Int_t itru2 = -1 ; - Int_t itru4 = -1 ; + Float_t maxn[] = {-1,-1,-1,-1} ; + Int_t mtru2 = -1 ; + Int_t mtrun = -1 ; //Find maximum summed amplitude of all the TRU //in a Module for(Int_t i = 0 ; i < fNTRU ; i++){ - if(max2[0] < (*ampmax2)(0,i) ){ - max2[0] = (*ampmax2)(0,i) ; // 2x2 summed max amplitude - max2[1] = (*ampmax2)(1,i) ; // corresponding phi position in TRU - max2[2] = (*ampmax2)(2,i) ; // corresponding eta position in TRU - max2[3] = (*ampmax2)(3,i) ; // corresponding most recent time - itru2 = i ; // TRU number + if(max2[0] < ampmax2(0,i) ){ + max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude + max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU + max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU + max2[3] = ampmax2(3,i) ; // corresponding most recent time + mtru2 = i ; // TRU number in module } - if(max4[0] < (*ampmax4)(0,i) ){ - max4[0] = (*ampmax4)(0,i) ; // 4x4 summed max amplitude - max4[1] = (*ampmax4)(1,i) ; // corresponding phi position in TRU - max4[2] = (*ampmax4)(2,i) ; // corresponding eta position in TRU - max4[3] = (*ampmax4)(3,i) ; // corresponding most recent time - itru4 = i ; // TRU number + if(maxn[0] < ampmaxn(0,i) ){ + maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude + maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU + maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU + maxn[3] = ampmaxn(3,i) ; // corresponding most recent time + mtrun = i ; // TRU number in module } } //Set max amplitude if larger than in other Modules Float_t maxtimeR2 = -1 ; - Float_t maxtimeR4 = -1 ; - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - AliPHOS * phos = gime->PHOS(); - Int_t nTimeBins = phos->GetRawFormatTimeBins() ; + Float_t maxtimeRn = -1 ; + // Create a shaper pulse object + AliPHOSPulseGenerator pulse ; + Int_t nTimeBins = pulse.GetRawFormatTimeBins() ; //Set max 2x2 amplitude and select L0 trigger if(max2[0] > f2x2MaxAmp ){ f2x2MaxAmp = max2[0] ; f2x2SM = iMod ; maxtimeR2 = max2[3] ; - GetCrystalPhiEtaIndexInModuleFromTRUIndex(itru2,static_cast(max2[1]),static_cast(max2[2]),f2x2CrystalPhi,f2x2CrystalEta,geom) ; + GetCrystalPhiEtaIndexInModuleFromTRUIndex(mtru2, + static_cast(max2[1]), + static_cast(max2[2]), + f2x2CrystalPhi,f2x2CrystalEta) ; + //Isolated patch? + if(fIsolateInModule) + fIs2x2Isol = IsPatchIsolated(0, iMod, mtru2, f2x2MaxAmp, f2x2CrystalPhi,f2x2CrystalEta) ; + else + fIs2x2Isol = IsPatchIsolated(0, iMod, mtru2, f2x2MaxAmp, static_cast(max2[1]), static_cast(max2[2])) ; + //Transform digit amplitude in Raw Samples - fADCValuesLow2x2 = new Int_t[nTimeBins]; - fADCValuesHigh2x2 = new Int_t[nTimeBins]; + if (fADCValuesLow2x2 == 0) { + fADCValuesLow2x2 = new Int_t[nTimeBins]; + } + if(!fADCValuesHigh2x2) fADCValuesHigh2x2 = new Int_t[nTimeBins]; + - phos->RawSampledResponse(maxtimeR2, f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ; + pulse.SetAmplitude(f2x2MaxAmp); + pulse.SetTZero(maxtimeR2); + pulse.MakeSamples(); + pulse.GetSamples(fADCValuesHigh2x2, fADCValuesLow2x2) ; //Set Trigger Inputs, compare ADC time bins until threshold is attained //Set L0 for(Int_t i = 0 ; i < nTimeBins ; i++){ - if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){ + if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold) { SetInput("PHOS_L0") ; break; } } -// for(Int_t i = 0 ; i < 256 ; i++) -// if(fADCValuesLow2x2[i]!=0||fADCValuesHigh2x2[i]!=0) -// cout<< "2x2 Time Bin "< f4x4MaxAmp ){ - f4x4MaxAmp = max4[0] ; - f4x4SM = iMod ; - maxtimeR4 = max4[3] ; - GetCrystalPhiEtaIndexInModuleFromTRUIndex(itru4,static_cast(max4[1]),static_cast(max4[2]),f4x4CrystalPhi,f4x4CrystalEta,geom) ; + //Set max nxn amplitude and select L1 triggers + if(maxn[0] > fnxnMaxAmp && fPatchSize > 0){ + fnxnMaxAmp = maxn[0] ; + fnxnSM = iMod ; + maxtimeRn = maxn[3] ; + GetCrystalPhiEtaIndexInModuleFromTRUIndex(mtrun, + static_cast(maxn[1]), + static_cast(maxn[2]), + fnxnCrystalPhi,fnxnCrystalEta) ; + //Isolated patch? + if(fIsolateInModule) + fIsnxnIsol = IsPatchIsolated(1, iMod, mtrun, fnxnMaxAmp, fnxnCrystalPhi, fnxnCrystalEta) ; + else + fIsnxnIsol = IsPatchIsolated(1, iMod, mtrun, fnxnMaxAmp, static_cast(maxn[1]), static_cast(maxn[2])) ; + //Transform digit amplitude in Raw Samples - fADCValuesHigh4x4 = new Int_t[nTimeBins]; - fADCValuesLow4x4 = new Int_t[nTimeBins]; - phos->RawSampledResponse(maxtimeR4, f4x4MaxAmp, fADCValuesHigh4x4, fADCValuesLow4x4) ; + if (fADCValuesHighnxn == 0) { + fADCValuesHighnxn = new Int_t[nTimeBins]; + fADCValuesLownxn = new Int_t[nTimeBins]; + } + + pulse.SetAmplitude(fnxnMaxAmp); + pulse.SetTZero(maxtimeRn); + pulse.MakeSamples(); + pulse.GetSamples(fADCValuesHighnxn, fADCValuesLownxn) ; //Set Trigger Inputs, compare ADC time bins until threshold is attained //SetL1 Low for(Int_t i = 0 ; i < nTimeBins ; i++){ - if(fADCValuesHigh4x4[i] >= fL1JetLowPtThreshold || fADCValuesLow4x4[i] >= fL1JetLowPtThreshold){ + if(fADCValuesHighnxn[i] >= fL1JetLowPtThreshold || fADCValuesLownxn[i] >= fL1JetLowPtThreshold){ SetInput("PHOS_JetLPt_L1") ; break; } } + //SetL1 Medium + for(Int_t i = 0 ; i < nTimeBins ; i++){ + if(fADCValuesHighnxn[i] >= fL1JetMediumPtThreshold || fADCValuesLownxn[i] >= fL1JetMediumPtThreshold){ + SetInput("PHOS_JetMPt_L1") ; + break; + } + } //SetL1 High for(Int_t i = 0 ; i < nTimeBins ; i++){ - if(fADCValuesHigh4x4[i] >= fL1JetHighPtThreshold || fADCValuesLow4x4[i] >= fL1JetHighPtThreshold){ + if(fADCValuesHighnxn[i] >= fL1JetHighPtThreshold || fADCValuesLownxn[i] >= fL1JetHighPtThreshold){ SetInput("PHOS_JetHPt_L1") ; break; } } -// for(Int_t i = 0 ; i < 256 ; i++) -// if(fADCValuesLow4x4[i]!=0||fADCValuesHigh4x4[i]!=0) -// cout<< "4x4 Time Bin "<GetRunLoader(); - //Getter - AliPHOSGetter * gime = AliPHOSGetter::Instance( rl->GetFileName() ) ; - //AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - //Get Geometry - const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; + fDigitsList = digits; + DoIt() ; +} + +//____________________________________________________________________________ +void AliPHOSTrigger::DoIt() +{ + // does the trigger job + + AliRunLoader* rl = AliRunLoader::GetRunLoader() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + + // Get PHOS Geometry object + AliPHOSGeometry *geom; + if (!(geom = AliPHOSGeometry::GetInstance())) + geom = AliPHOSGeometry::GetInstance("IHEP",""); //Define parameters Int_t nModules = geom->GetNModules(); + fNCrystalsPhi = geom->GetNPhi()/fNTRUPhi ;// 64/4=16 + fNCrystalsZ = geom->GetNZ()/fNTRUZ ;// 56/2=28 //Intialize data members each time the trigger is called in event loop f2x2MaxAmp = -1; f2x2CrystalPhi = -1; f2x2CrystalEta = -1; - f4x4MaxAmp = -1; f4x4CrystalPhi = -1; f4x4CrystalEta = -1; + fnxnMaxAmp = -1; fnxnCrystalPhi = -1; fnxnCrystalEta = -1; //Take the digits list if simulation if(fSimulation) - fDigitsList = gime->Digits() ; + fDigitsList = phosLoader->Digits() ; if(!fDigitsList) AliFatal("Digits not found !") ; //Fill TRU Matrix - TClonesArray * amptrus = new TClonesArray("TMatrixD",1000); - TClonesArray * timeRtrus = new TClonesArray("TMatrixD",1000); - FillTRU(fDigitsList,geom,amptrus, timeRtrus) ; +// TClonesArray * amptrus = new TClonesArray("TMatrixD",1000); +// TClonesArray * ampmods = new TClonesArray("TMatrixD",1000); +// TClonesArray * timeRtrus = new TClonesArray("TMatrixD",1000); + FillTRU(fDigitsList,geom) ; //Do Crystal Sliding and select Trigger //Initialize varible that will contain maximum amplitudes and //its corresponding cell position in eta and phi, and time. - TMatrixD * ampmax2 = new TMatrixD(4,fNTRU) ; - TMatrixD * ampmax4 = new TMatrixD(4,fNTRU) ; + TMatrixD ampmax2(4,fNTRU) ; + TMatrixD ampmaxn(4,fNTRU) ; + + for(Int_t imod = 0 ; imod < nModules ; imod++) { - for(Int_t imod = 1 ; imod <= nModules ; imod++) { - //Do 2x2 and 4x4 sums, select maximums. - MakeSlidingCell(amptrus, timeRtrus, imod, ampmax2, ampmax4, geom); + //Do 2x2 and nxn sums, select maximums. + MakeSlidingCell(imod, ampmax2, ampmaxn); //Set the trigger - SetTriggers(imod,ampmax2,ampmax4, geom) ; + SetTriggers(imod,ampmax2,ampmaxn) ; } + + fAmptrus->Delete(); +// delete amptrus; amptrus=0; + fAmpmods->Delete(); +// delete ampmods; ampmods=0; + fTimeRtrus->Delete(); +// delete timeRtrus; timeRtrus=0; + //Print(); + }