From 64df000df0708eeec3eb7ab1ae7e29ee7fedd0ee Mon Sep 17 00:00:00 2001 From: hristov Date: Mon, 26 Feb 2007 16:21:03 +0000 Subject: [PATCH] Adding information from the PHOS trigger (Gustavo, Yves) --- PHOS/AliPHOSReconstructor.cxx | 73 +++++++++++- PHOS/AliPHOSTrigger.cxx | 209 +++++++++++++++++++++++++--------- PHOS/AliPHOSTrigger.h | 63 +++++++--- PWG4/AliTriggerPHOS.cxx | 14 +-- STEER/AliESD.cxx | 28 ++++- STEER/AliESD.h | 17 ++- STEER/AliESDCaloCluster.cxx | 20 +++- STEER/AliESDCaloCluster.h | 9 +- 8 files changed, 348 insertions(+), 85 deletions(-) diff --git a/PHOS/AliPHOSReconstructor.cxx b/PHOS/AliPHOSReconstructor.cxx index cdd3cae2cab..dcc2b5c996b 100644 --- a/PHOS/AliPHOSReconstructor.cxx +++ b/PHOS/AliPHOSReconstructor.cxx @@ -34,8 +34,9 @@ #include "AliPHOSGetter.h" #include "AliPHOSTracker.h" #include "AliRawReader.h" +#include "AliPHOSTrigger.h" +#include "AliPHOSGeometry.h" - ClassImp(AliPHOSReconstructor) Bool_t AliPHOSReconstructor::fgDebug = kFALSE ; @@ -121,6 +122,62 @@ void AliPHOSReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const AliDebug(2,Form("%d digits and %d rec. particles in event %d, option %s",gime->Digits()->GetEntries(),nOfRecParticles,eventNumber,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() ; + + AliPHOSGeometry * geom = gime->PHOSGeometry(); + + 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,0,iCrystalPhi2x2,iCrystalEta2x2}; + Int_t iAbsId2x2 =-1; + Int_t iRelIdnxn []= {iSMnxn,0,iCrystalPhinxn,iCrystalEtanxn}; + Int_t iAbsIdnxn =-1; + TVector3 pos2x2(-1,-1,-1); + TVector3 posnxn(-1,-1,-1); + geom->RelToAbsNumbering(iRelId2x2, iAbsId2x2); + geom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn); + geom->RelPosInAlice(iAbsId2x2, pos2x2); + geom->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); + + //###################################### + + //Fill CaloClusters for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) { AliPHOSRecParticle * rp = dynamic_cast(recParticles->At(recpart)); if (Debug()) @@ -130,13 +187,14 @@ void AliPHOSReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const AliPHOSEmcRecPoint *emcRP = gime->EmcRecPoint(ts->GetEmcIndex()); AliESDCaloCluster *ec = new AliESDCaloCluster() ; - // fills the ESDCaloCluster + // fills the ESDCaloCluster Float_t xyz[3]; for (Int_t ixyz=0; ixyz<3; ixyz++) xyz[ixyz] = rp->GetPos()[ixyz]; AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2])); + Int_t digitMult = emcRP->GetDigitsMultiplicity(); Int_t *digitsList = emcRP->GetDigitsList(); UShort_t *amplList = new UShort_t[digitMult]; @@ -156,7 +214,6 @@ void AliPHOSReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const ec->SetClusterEnergy(rp->Energy()); //total particle energy ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion ec->SetPid (rp->GetPID()) ; //array of particle identification - ec->SetPrimaryIndex (rp->GetPrimaryIndex());//index of primary particle (for simulations) ec->SetM02(emcRP->GetM2x()) ; //second moment M2x ec->SetM20(emcRP->GetM2z()) ; //second moment M2z ec->SetNExMax(emcRP->GetNExMax()); //number of local maxima @@ -168,6 +225,16 @@ void AliPHOSReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const ec->SetClusterChi2(-1); //not yet implemented ec->SetM11(-1) ; //not yet implemented + //Primaries + ec->SetPrimaryIndex(rp->GetPrimaryIndex()); + Int_t primMult = 0; + Int_t *primInts = emcRP->GetPrimaries(primMult); + ec->SetNumberOfPrimaries(primMult); //primary multiplicity + UShort_t *primList = new UShort_t[primMult]; + for (Int_t ipr=0; iprSetListOfPrimaries(primList); //primary List for a cluster + // add the track to the esd object esd->AddCaloCluster(ec); delete ec; diff --git a/PHOS/AliPHOSTrigger.cxx b/PHOS/AliPHOSTrigger.cxx index 2496928229c..f9422e3ad50 100644 --- a/PHOS/AliPHOSTrigger.cxx +++ b/PHOS/AliPHOSTrigger.cxx @@ -20,9 +20,9 @@ // Digits are grouped in TRU's (Trigger Units). A TRU consist of 16x28 // crystals ordered fNTRUPhi x fNTRUZ. The algorithm searches all possible // 2x2 and nxn (n multiple of 4) crystal combinations per each TRU, adding the -// digits amplitude and finding the maximum. 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. +// digits amplitude and finding the maximum. Iti is found is maximum 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 @@ -30,9 +30,10 @@ // tr->SetL0Threshold(100); // tr->SetL1JetLowPtThreshold(1000); // 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) ////////////////////////////////////////////////////////////////////////////// @@ -60,7 +61,14 @@ AliPHOSTrigger::AliPHOSTrigger() fADCValuesHighnxn(0), fADCValuesLownxn(0), fADCValuesHigh2x2(0), fADCValuesLow2x2(0), fDigitsList(0), fL0Threshold(50), fL1JetLowPtThreshold(200), fL1JetHighPtThreshold(500), - fNTRU(8), fNTRUZ(2), fNTRUPhi(4), fPatchSize(1), fSimulation(kTRUE) + 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) { //ctor fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins]; @@ -96,7 +104,17 @@ AliPHOSTrigger::AliPHOSTrigger(const AliPHOSTrigger & trig) : fNTRU(trig.fNTRU), fNTRUZ(trig.fNTRUZ), fNTRUPhi(trig.fNTRUPhi), - fPatchSize(trig.fPatchSize), fSimulation(trig.fSimulation) + 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) { // cpy ctor } @@ -136,8 +154,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; @@ -145,10 +161,10 @@ 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++){ + TMatrixD * amptrus = new TMatrixD(fNCrystalsPhi,fNCrystalsZ) ; + TMatrixD * timeRtrus = new TMatrixD(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; } @@ -179,12 +195,12 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry //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 ; @@ -194,8 +210,8 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry TMatrixD * timeRtrus = dynamic_cast(timeRmatrix->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 ; @@ -204,7 +220,7 @@ void AliPHOSTrigger::FillTRU(const TClonesArray * digits, const AliPHOSGeometry } //______________________________________________________________________ -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. @@ -215,22 +231,89 @@ 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 ; } + +//____________________________________________________________________________ +Bool_t AliPHOSTrigger::IsPatchIsolated(Int_t iPatchType, const TClonesArray * amptrus, const Int_t mtru, const Int_t imod, const Float_t *maxarray) { + + //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 with maximum amplitude patch. + Int_t itru = mtru+imod*fNTRU ; //number of tru, min 0 max 8*5. + TMatrixD * amptru = dynamic_cast(amptrus->At(itru)) ; + + //Define patch cells + Int_t isolcells = fIsolPatchSize*(1+iPatchType); + Int_t ipatchcells = 2*(1+fPatchSize*iPatchType); + Int_t minrow = static_cast(maxarray[1]) - isolcells; + Int_t mincol = static_cast(maxarray[2]) - isolcells; + Int_t maxrow = static_cast(maxarray[1]) + isolcells + ipatchcells; + Int_t maxcol = static_cast(maxarray[2]) + 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 > fNCrystalsPhi || maxcol > fNCrystalsZ){ + AliDebug(1,Form("Out of 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 += (*amptru)(irow,icol); + + AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxarray[0], amp)); + + if(amp < maxarray[0]){ + AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxarray[0], amp)); + return kFALSE; + } + else + amp-=maxarray[0]; //Calculate energy in isolation patch that do not comes from maximum patch. + + AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxarray[0], 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 TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t imod, TMatrixD *ampmax2, TMatrixD *ampmaxn, const AliPHOSGeometry *geom){ +void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t imod, TMatrixD *ampmax2, TMatrixD *ampmaxn){ //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 ampn = 0 ; for(Int_t i = 0; i < 4; i++){ @@ -242,23 +325,22 @@ void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClones //Create matrix that will contain 2x2 amplitude sums //used to calculate the nxn 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; + TMatrixD * tru2x2 = new TMatrixD(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++){ + for(Int_t itru = 0 + imod * fNTRU ; itru < (imod+1)*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 + 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 nxn sums (*tru2x2)(irow/2,icol/2) = amp2 ; //Select 2x2 maximum sums to select L0 @@ -285,10 +367,10 @@ void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClones //Sliding nxn, add nxn amplitudes (OVERLAP) if(fPatchSize > 0){ - for(Int_t irow = 0 ; irow < nCrystalsPhi/2; irow++){ - for(Int_t icol = 0 ; icol < nCrystalsZ/2 ; icol++){ + for(Int_t irow = 0 ; irow < fNCrystalsPhi/2; irow++){ + for(Int_t icol = 0 ; icol < fNCrystalsZ/2 ; icol++){ ampn = 0; - if( (irow+fPatchSize) < nCrystalsPhi/2 && (icol+fPatchSize) < nCrystalsZ/2){//Avoid exit the TRU + 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); @@ -308,7 +390,7 @@ void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClones 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) < nCrystalsPhi && (coln+j) < nCrystalsZ/2){//Avoid exit the TRU + 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); @@ -325,6 +407,8 @@ void AliPHOSTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClones } } } + + //____________________________________________________________________________ void AliPHOSTrigger::Print(const Option_t * opt) const { @@ -339,12 +423,15 @@ 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( " -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",4*fPatchSize, 4*fPatchSize); + 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, fnxnCrystalEta, fnxnCrystalEta+4) ; + 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( " Threshold for LO %10.1f\n", fL0Threshold) ; @@ -367,16 +454,16 @@ void AliPHOSTrigger::Print(const Option_t * opt) const } //____________________________________________________________________________ -void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD * ampmax2, const TMatrixD * ampmaxn, const AliPHOSGeometry *geom) +void AliPHOSTrigger::SetTriggers(const TClonesArray * amptrus, const Int_t iMod, const TMatrixD * ampmax2, const TMatrixD * ampmaxn) { //Checks the 2x2 and nxn maximum amplitude per each TRU and compares - //with the different L0 and L1 triggers thresholds + //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 maxn[] = {-1,-1,-1,-1} ; - Int_t itru2 = -1 ; - Int_t itrun = -1 ; + Int_t mtru2 = -1 ; + Int_t mtrun = -1 ; //Find maximum summed amplitude of all the TRU @@ -387,14 +474,14 @@ void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD * ampmax2, con 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 + mtru2 = i ; // TRU number in module } 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 - itrun = i ; // TRU number + mtrun = i ; // TRU number in module } } @@ -410,11 +497,14 @@ void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD * ampmax2, con f2x2MaxAmp = max2[0] ; f2x2SM = iMod ; maxtimeR2 = max2[3] ; - GetCrystalPhiEtaIndexInModuleFromTRUIndex(itru2, + GetCrystalPhiEtaIndexInModuleFromTRUIndex(mtru2, static_cast(max2[1]), static_cast(max2[2]), - f2x2CrystalPhi,f2x2CrystalEta,geom) ; + f2x2CrystalPhi,f2x2CrystalEta) ; + //Isolated patch? + fIs2x2Isol = IsPatchIsolated(0, amptrus, mtru2, iMod, max2) ; + //Transform digit amplitude in Raw Samples fADCValuesLow2x2 = new Int_t[nTimeBins]; fADCValuesHigh2x2 = new Int_t[nTimeBins]; @@ -435,15 +525,18 @@ void AliPHOSTrigger::SetTriggers(const Int_t iMod, const TMatrixD * ampmax2, con } //Set max nxn amplitude and select L1 triggers - if(maxn[0] > fnxnMaxAmp ){ + if(maxn[0] > fnxnMaxAmp && fPatchSize > 0){ fnxnMaxAmp = maxn[0] ; fnxnSM = iMod ; maxtimeRn = maxn[3] ; - GetCrystalPhiEtaIndexInModuleFromTRUIndex(itrun, + GetCrystalPhiEtaIndexInModuleFromTRUIndex(mtrun, static_cast(maxn[1]), static_cast(maxn[2]), - fnxnCrystalPhi,fnxnCrystalEta,geom) ; + fnxnCrystalPhi,fnxnCrystalEta) ; + //Isolated patch? + fIsnxnIsol = IsPatchIsolated(1, amptrus, mtrun, iMod, maxn) ; + //Transform digit amplitude in Raw Samples fADCValuesHighnxn = new Int_t[nTimeBins]; fADCValuesLownxn = new Int_t[nTimeBins]; @@ -486,6 +579,8 @@ void AliPHOSTrigger::Trigger() //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; @@ -509,10 +604,14 @@ void AliPHOSTrigger::Trigger() TMatrixD * ampmax2 = new TMatrixD(4,fNTRU) ; TMatrixD * ampmaxn = new TMatrixD(4,fNTRU) ; - for(Int_t imod = 1 ; imod <= nModules ; imod++) { + for(Int_t imod = 0 ; imod < nModules ; imod++) { + //Do 2x2 and nxn sums, select maximums. - MakeSlidingCell(amptrus, timeRtrus, imod, ampmax2, ampmaxn, geom); + MakeSlidingCell(amptrus, timeRtrus, imod, ampmax2, ampmaxn); //Set the trigger - SetTriggers(imod,ampmax2,ampmaxn, geom) ; + SetTriggers(amptrus,imod,ampmax2,ampmaxn) ; } + + //Print(); + } diff --git a/PHOS/AliPHOSTrigger.h b/PHOS/AliPHOSTrigger.h index 5b7fe03e36a..6d5d1878a9b 100644 --- a/PHOS/AliPHOSTrigger.h +++ b/PHOS/AliPHOSTrigger.h @@ -8,24 +8,26 @@ //____________________________________________________________ // Class for trigger analysis. + +// +//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, SUBATECH, CERN) // Digits are grouped in TRU's (Trigger Units). A TRU consist of 16x28 // crystals ordered fNTRUPhi x fNTRUZ. The algorithm searches all possible -// 2x2 and nxn (n multiple of 4) crystal combinations per each TRU, adding the -// digits amplitude and finding the maximum. 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. +// 2x2 and nxn (n multiple of 4) crystal combinations per each TRU, adding the +// digits amplitude and finding the maximum. Iti is found is maximum 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 -// AliPHOSTrigger *tr = new AliPHOSTrigger();//Init Trigger +// AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger // tr->SetL0Threshold(100); // tr->SetL1JetLowPtThreshold(1000); // tr->SetL1JetHighPtThreshold(20000); +// .... // tr->Trigger(); //Execute Trigger -// tr->Print(""); //Print results -// -//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, SUBATECH, CERN) - +// tr->Print(""); //Print data members after calculation. +// // --- ROOT system --- class TClonesArray ; @@ -64,7 +66,7 @@ class AliPHOSTrigger : public AliTriggerDetector { Int_t * GetADCValuesLowGainMaxnxnSum() {return fADCValuesLownxn; } Int_t * GetADCValuesHighGainMaxnxnSum() {return fADCValuesHighnxn; } - void GetCrystalPhiEtaIndexInModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru,Int_t &ietaMod,Int_t &iphiMod, const AliPHOSGeometry *geom) const ; + void GetCrystalPhiEtaIndexInModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru,Int_t &ietaMod,Int_t &iphiMod) const ; Float_t GetL0Threshold() const {return fL0Threshold ; } Float_t GetL1JetLowPtThreshold() const {return fL1JetLowPtThreshold ; } @@ -74,7 +76,16 @@ class AliPHOSTrigger : public AliTriggerDetector { Int_t GetNTRUZ() const {return fNTRUZ ; } Int_t GetNTRUPhi() const {return fNTRUPhi ; } - Float_t GetPatchSize() const {return fPatchSize ; } + Int_t GetPatchSize() const {return fPatchSize ; } + Int_t GetIsolPatchSize() const {return fIsolPatchSize ; } + + Float_t Get2x2AmpOutOfPatch() const {return f2x2AmpOutOfPatch; } + Float_t GetnxnAmpOutOfPatch() const {return fnxnAmpOutOfPatch; } + Float_t Get2x2AmpOutOfPatchThres() const {return f2x2AmpOutOfPatchThres; } + Float_t GetnxnAmpOutOfPatchThres() const {return fnxnAmpOutOfPatchThres; } + Bool_t Is2x2Isol() const {return fIs2x2Isol; } + Bool_t IsnxnIsol() const {return fIsnxnIsol; } + Bool_t IsSimulation() const {return fSimulation ; } //Setters @@ -95,6 +106,9 @@ class AliPHOSTrigger : public AliTriggerDetector { {fL1JetHighPtThreshold = amp ; } void SetPatchSize(Int_t ps) {fPatchSize = ps ; } + void SetIsolPatchSize(Int_t ps) {fIsolPatchSize = ps ; } + void Set2x2AmpOutOfPatchThres(Float_t th) { f2x2AmpOutOfPatchThres = th; } + void SetnxnAmpOutOfPatchThres(Float_t th) { fnxnAmpOutOfPatchThres = th; } void SetSimulation(Bool_t sim ) {fSimulation = sim ; } private: @@ -103,9 +117,11 @@ class AliPHOSTrigger : public AliTriggerDetector { void FillTRU(const TClonesArray * digits, const AliPHOSGeometry * geom, TClonesArray * amptru, TClonesArray * timeRtru) const ; - void MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, Int_t mod, TMatrixD *ampmax2, TMatrixD *ampmaxn, const AliPHOSGeometry *geom) ; + Bool_t IsPatchIsolated(Int_t iPatchType, const TClonesArray * amptrus, const Int_t mtru, const Int_t imod, const Float_t *maxarray) ; + + void MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, Int_t mod, TMatrixD *ampmax2, TMatrixD *ampmaxn) ; - void SetTriggers(Int_t iMod, const TMatrixD *ampmax2,const TMatrixD *ampmaxn, const AliPHOSGeometry *geom) ; + void SetTriggers(const TClonesArray * amptrus, Int_t iMod, const TMatrixD *ampmax2,const TMatrixD *ampmaxn) ; private: @@ -129,11 +145,26 @@ class AliPHOSTrigger : public AliTriggerDetector { Float_t fL1JetLowPtThreshold ; //! L1 Low pT trigger threshold Float_t fL1JetHighPtThreshold ; //! L1 High pT trigger threshold - Int_t fNTRU ; //! Number of TRUs per module + Int_t fNTRU ; //! Number of TRUs per module Int_t fNTRUZ ; //! Number of crystal rows per Z in one TRU - Int_t fNTRUPhi ; //! Number of crystal rows per Phi in one TRU + Int_t fNTRUPhi ; //! Number of crystal rows per Phi in one TRU + Int_t fNCrystalsPhi; //!Number of rows in a TRU + Int_t fNCrystalsZ; //!Number of columns in a TRU + Int_t fPatchSize; //! Trigger patch factor, to be multiplied to 2x2 cells - // 0 means 2x2, 1 means nxn, 2 means 8x8 ... + // 0 means 2x2, 1 means 4x4, 2 means 6x6 ... + Int_t fIsolPatchSize ; //Isolation patch size, number of rows or columns to add to + //the 2x2 or nxn maximum amplitude patch. + //1 means a patch around max amplitude of 2x2 of 4x4 and around + //max ampl patch of 4x4 of 8x8 + + Float_t f2x2AmpOutOfPatch; //Amplitude in isolation cone minus maximum amplitude of the reference patch + Float_t fnxnAmpOutOfPatch; + Float_t f2x2AmpOutOfPatchThres; //Threshold to select a trigger as isolated on f2x2AmpOutOfPatch value + Float_t fnxnAmpOutOfPatchThres; + Float_t fIs2x2Isol; //Patch is isolated if f2x2AmpOutOfPatchThres threshold is passed + Float_t fIsnxnIsol ; + Bool_t fSimulation ; //! Flag to do the trigger during simulation or reconstruction ClassDef(AliPHOSTrigger,4) } ; diff --git a/PWG4/AliTriggerPHOS.cxx b/PWG4/AliTriggerPHOS.cxx index e69bc3c8c4b..5c08f07c81b 100644 --- a/PWG4/AliTriggerPHOS.cxx +++ b/PWG4/AliTriggerPHOS.cxx @@ -125,13 +125,13 @@ void AliTriggerPHOS::Exec(Option_t *) const Float_t aNNO = static_cast(triggerAmplitudes->At(3)) ; // trigger position - const TArrayI * triggerCells = fESD->GetPHOSTriggerCells(); - const Float_t p22 = static_cast(triggerCells->At(0)) ; - const Float_t phi22 = static_cast(triggerCells->At(1)) ; - const Float_t eta22 = static_cast(triggerCells->At(2)) ; - const Float_t pNN = static_cast(triggerCells->At(3)) ; - const Float_t phiNN = static_cast(triggerCells->At(4)) ; - const Float_t etaNN = static_cast(triggerCells->At(5)) ; + const TArrayF * triggerPosition = fESD->GetPHOSTriggerPosition(); + const Float_t p22 = static_cast(triggerPosition->At(0)) ; + const Float_t phi22 = static_cast(triggerPosition->At(1)) ; + const Float_t eta22 = static_cast(triggerPosition->At(2)) ; + const Float_t pNN = static_cast(triggerPosition->At(3)) ; + const Float_t phiNN = static_cast(triggerPosition->At(4)) ; + const Float_t etaNN = static_cast(triggerPosition->At(5)) ; Int_t firstPhosCluster = fESD->GetFirstPHOSCluster() ; const Int_t numberOfPhosClusters = fESD->GetNumberOfPHOSClusters() ; diff --git a/STEER/AliESD.cxx b/STEER/AliESD.cxx index 9858babfeb6..a60b406669e 100644 --- a/STEER/AliESD.cxx +++ b/STEER/AliESD.cxx @@ -60,10 +60,15 @@ AliESD::AliESD(): fCaloClusters("AliESDCaloCluster",10000), fEMCALClusters(0), fFirstEMCALCluster(-1), + fEMCALTriggerPosition(0x0), + fEMCALTriggerAmplitudes(0x0), fPHOSClusters(0), fFirstPHOSCluster(-1), + fPHOSTriggerPosition(0x0), + fPHOSTriggerAmplitudes(0x0), fESDFMD(0x0), fESDVZERO(0x0) + { for (Int_t i=0; i<24; i++) { fT0time[i] = 0; @@ -104,8 +109,12 @@ AliESD::AliESD(const AliESD& esd): fCaloClusters(*((TClonesArray*)esd.fCaloClusters.Clone())), fEMCALClusters(esd.fEMCALClusters), fFirstEMCALCluster(esd.fFirstEMCALCluster), + fEMCALTriggerPosition(esd. fEMCALTriggerPosition), + fEMCALTriggerAmplitudes(esd.fEMCALTriggerAmplitudes), fPHOSClusters(esd.fPHOSClusters), fFirstPHOSCluster(esd.fFirstPHOSCluster), + fPHOSTriggerPosition(esd.fPHOSTriggerPosition), + fPHOSTriggerAmplitudes(esd.fPHOSTriggerAmplitudes), fESDFMD(esd.fESDFMD), fESDVZERO(esd.fESDVZERO) { @@ -157,7 +166,11 @@ AliESD & AliESD::operator=(const AliESD& source) { fFirstPHOSCluster = source.fFirstPHOSCluster; fESDFMD = source.fESDFMD; fESDVZERO = source.fESDVZERO; - + fEMCALTriggerPosition=source. fEMCALTriggerPosition; + fEMCALTriggerAmplitudes=source.fEMCALTriggerAmplitudes; + fPHOSTriggerPosition=source.fPHOSTriggerPosition; + fPHOSTriggerAmplitudes=source.fPHOSTriggerAmplitudes; + for (Int_t i=0; i<24; i++) { fT0time[i] = source.fT0time[i]; fT0amplitude[i] = source.fT0amplitude[i]; @@ -186,6 +199,15 @@ AliESD::~AliESD() fCaloClusters.Delete(); delete fESDFMD; delete fESDVZERO; +// fEMCALTriggerPosition->Delete(); +// fEMCALTriggerAmplitudes->Delete(); +// fPHOSTriggerPosition->Delete(); +// fPHOSTriggerAmplitudes->Delete(); +// delete fEMCALTriggerPosition; +// delete fEMCALTriggerAmplitudes; +// delete fPHOSTriggerPosition; +// delete fPHOSTriggerAmplitudes; + } //______________________________________________________________________________ @@ -224,6 +246,10 @@ void AliESD::Reset() fPHOSClusters=0; fFirstPHOSCluster=-1; if (fESDFMD) fESDFMD->Clear(); +// fEMCALTriggerPosition->Clear(); +// fEMCALTriggerAmplitudes->Clear(); +// fPHOSTriggerPosition->Clear(); +// fPHOSTriggerAmplitudes->Clear(); } Int_t AliESD::AddV0(const AliESDv0 *v) { diff --git a/STEER/AliESD.h b/STEER/AliESD.h index cee3b1e4a21..403da7b23c9 100644 --- a/STEER/AliESD.h +++ b/STEER/AliESD.h @@ -16,6 +16,7 @@ #include #include +#include #include "AliESDMuonTrack.h" #include "AliESDPmdTrack.h" @@ -121,6 +122,11 @@ public: return fCaloClusters.GetEntriesFast()-1; } + void AddPHOSTriggerPosition(TArrayF array) { fPHOSTriggerPosition = new TArrayF(array) ; } + void AddPHOSTriggerAmplitudes(TArrayF array) { fPHOSTriggerAmplitudes = new TArrayF(array) ; } + void AddEMCALTriggerPosition(TArrayF array) { fEMCALTriggerPosition = new TArrayF(array) ; } + void AddEMCALTriggerAmplitudes(TArrayF array){ fEMCALTriggerAmplitudes = new TArrayF(array) ; } + void SetVertex(const AliESDVertex *vertex) { new (&fSPDVertex) AliESDVertex(*vertex); } @@ -158,11 +164,16 @@ public: void SetNumberOfEMCALClusters(Int_t clus) {fEMCALClusters = clus;} Int_t GetFirstEMCALCluster() const {return fFirstEMCALCluster;} void SetFirstEMCALCluster(Int_t index) {fFirstEMCALCluster = index;} + TArrayF *GetEMCALTriggerPosition() const {return fEMCALTriggerPosition;} + TArrayF *GetEMCALTriggerAmplitudes() const {return fEMCALTriggerAmplitudes;} Int_t GetNumberOfPHOSClusters() const {return fPHOSClusters;} void SetNumberOfPHOSClusters(Int_t part) { fPHOSClusters = part ; } void SetFirstPHOSCluster(Int_t index) { fFirstPHOSCluster = index ; } Int_t GetFirstPHOSCluster() const { return fFirstPHOSCluster ; } + TArrayF *GetPHOSTriggerPosition() const {return fPHOSTriggerPosition;} + TArrayF *GetPHOSTriggerAmplitudes() const {return fPHOSTriggerAmplitudes;} + Float_t GetT0zVertex() const {return fT0zVertex;} void SetT0zVertex(Float_t z) {fT0zVertex=z;} @@ -242,10 +253,14 @@ protected: TClonesArray fCaloClusters; // Calorimeter clusters for PHOS/EMCAL Int_t fEMCALClusters; // Number of EMCAL clusters (subset of caloclusters) Int_t fFirstEMCALCluster; // First EMCAL cluster in the fCaloClusters list + TArrayF *fEMCALTriggerPosition; ///(x,y,z of 2x2 and x,y,z of nxn) not position of centroid but of patch corner + TArrayF *fEMCALTriggerAmplitudes; //(2x2 max ampl, 2x2 amp out of patch, nxn max ampl, nxn amp out of patch) Int_t fPHOSClusters; // Number of PHOS clusters (subset of caloclusters) Int_t fFirstPHOSCluster; // First PHOS cluster in the fCaloClusters list - + TArrayF *fPHOSTriggerPosition; //(x,y,z of 2x2 and x,y,z of nxn), not position of centroid but of patch corner + TArrayF *fPHOSTriggerAmplitudes; //(2x2 max ampl, 2x2 amp out of patch, nxn max ampl, nxn amp out of patch) + AliESDFMD *fESDFMD; // FMD object containing rough multiplicity AliESDVZERO *fESDVZERO; // VZERO object containing rough multiplicity diff --git a/STEER/AliESDCaloCluster.cxx b/STEER/AliESDCaloCluster.cxx index cfd2b2daa3f..3c491c75e97 100644 --- a/STEER/AliESDCaloCluster.cxx +++ b/STEER/AliESDCaloCluster.cxx @@ -44,6 +44,8 @@ AliESDCaloCluster::AliESDCaloCluster() : fM11(0), fNExMax(0), fEmcCpvDistance(9999), + fNumberOfPrimaries(-1), + fListOfPrimaries(0x0), fNumberOfDigits(0), fDigitAmplitude(0x0), fDigitTime(0x0), @@ -72,6 +74,8 @@ AliESDCaloCluster::AliESDCaloCluster(const AliESDCaloCluster& clus) : fM11(clus.fM11), fNExMax(clus.fNExMax), fEmcCpvDistance(clus.fEmcCpvDistance), + fNumberOfPrimaries(clus.fNumberOfPrimaries), + fListOfPrimaries(0x0), fNumberOfDigits(clus.fNumberOfDigits), fDigitAmplitude(0x0), fDigitTime(0x0), @@ -102,6 +106,11 @@ AliESDCaloCluster::AliESDCaloCluster(const AliESDCaloCluster& clus) : for (Int_t i=0; i