X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALJetFinderAlgoOmni.cxx;h=d48a1debf5a51922a3baecabb946de36b4ea9513;hb=c2ae193faae2f3bc31f2bbaaac5a6f1b6e968e87;hp=fe69f0f6c3ae4871cb486965162068e2ff1b9ae0;hpb=ee6b678fd901cf80c9ad3b23e94efbb0496980ab;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALJetFinderAlgoOmni.cxx b/EMCAL/AliEMCALJetFinderAlgoOmni.cxx index fe69f0f6c3a..d48a1debf5a 100644 --- a/EMCAL/AliEMCALJetFinderAlgoOmni.cxx +++ b/EMCAL/AliEMCALJetFinderAlgoOmni.cxx @@ -17,7 +17,48 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ +/* + +$Log$ +Revision 1.14 2004/04/02 17:11:33 mhorner +Marco's bug - fixes implemented + +Revision 1.13 2004/03/26 01:00:38 mhorner +Implementing Marco's optimisations + +Revision 1.12 2004/03/15 19:59:37 mhorner +Pyhtia comparison extended + +Revision 1.11 2004/03/09 17:06:38 mhorner +Made more robust + +Revision 1.10 2004/03/06 01:56:09 mhorner +Pythai comp code + +Revision 1.9 2004/02/23 20:37:32 mhorner +changed geometry call + +Revision 1.8 2004/01/29 23:28:44 mhorner +Jet Finder - hard coded geom parameters removed + +Revision 1.7 2004/01/21 22:27:47 mhorner +Cleaning up coding conventions + +Revision 1.6 2003/10/28 13:54:30 schutz +Compilation warnings fixed + +Revision 1.5 2003/09/23 13:31:41 mhorner +Changed coordinate system + +Revision 1.4 2003/09/19 13:16:20 mhorner +Added additional jet energy info + + +Revision 1.3 2003/09/04 12:49:56 mhorner +Changed hadron correction and added saving EMCAL and track contributions + +*/ + //*--Author: Sarah Blyth (LBL) //*--Based on UA1 jet algorithm from LUND JETSET called from EMC-erj @@ -28,9 +69,9 @@ #include "AliEMCALJetFinderAlgo.h" #include "AliEMCALJetFinderAlgoOmni.h" #include "AliEMCALJetFinderAlgoUA1Unit.h" +#include "AliEMCALGetter.h" #include "AliEMCALGeometry.h" #include "AliEMCAL.h" -#include "AliEMCALGetter.h" #include "AliEMCALDigit.h" #include "TParticle.h" #include "AliRun.h" @@ -45,12 +86,15 @@ ClassImp(AliEMCALJetFinderAlgoOmni) //Default constructor if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); +// AliEMCALGetter * gime = AliEMCALGetter::Instance() ; +// AliEMCALGeometry * geom = gime->EMCALGeometry(); + AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","EMCAL"); fNumIter = 0; - fNumUnits = 13824; //Number of towers in EMCAL + fNumUnits = geom->GetNTowers(); //Number of towers in EMCAL fESeed = 5.0; //Default value fConeRad = 0.3; //Default value fJetEMin = 10.0; //Default value - fEtMin = 0.28; //Default value + fEtMin = 0.0; //Default value fMinMove = 0.05; //From original UA1 JetFinder fMaxMove = 0.15; //From original UA1 JetFinder fBGMaxMove = 0.035; //From original UA1 JetFinder @@ -126,7 +170,8 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); void AliEMCALJetFinderAlgoOmni::InitUnitArray() { //Initialises unit arrays - if(fArrayInitialised) delete[] fUnit; + if(fArrayInitialised) delete [] fUnit; + if(fArrayInitialised) delete [] fUnitNoCuts; fUnit = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits]; fUnitNoCuts = new AliEMCALJetFinderAlgoUA1Unit[fNumUnits]; fArrayInitialised = 1; @@ -134,16 +179,22 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); void AliEMCALJetFinderAlgoOmni::FillUnitArray(AliEMCALJetFinderAlgoUA1FillUnitFlagType_t flag) { + // Fill the unit array if (fDebug>1) Info("FillUnitArray","Beginning FillUnitArray"); - // AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); - + // AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL"); // if (pEMCAL){ // AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), ""); // }else // { - //AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance("EMCAL_5655_21", ""); - AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry() ; - // } + + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + AliEMCALGeometry * geom; + if (gime) + geom = gime->EMCALGeometry(); + else + geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","EMCAL"); + + // } AliEMCALJetFinderAlgoUA1FillUnitFlagType_t option = flag; Int_t numTracks, numDigits; @@ -156,7 +207,8 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); numDigits = fInputPointer->GetNDigits(); TParticle *myPart; AliEMCALDigit *myDigit; - + if (!fPythiaComparison) + { //Fill units with Track info if appropriate if(option==kFillTracksOnly || option ==kFillAll) { @@ -170,17 +222,20 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); Float_t unitEnergy = fUnit[towerID-1].GetUnitEnergy(); Float_t unitEnergyNoCuts = fUnitNoCuts[towerID-1].GetUnitEnergy(); - /* - //OLD WAY: //Do Hadron Correction + + //OLD WAY: //Do Hadron Correction if(fHadCorr != 0) { Double_t fullP = myPart->P(); Double_t hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta); unitEnergy -= hCEnergy*TMath::Sin(myPart->Theta()); + unitEnergyNoCuts -= hCEnergy*TMath::Sin(myPart->Theta()); fUnit[towerID-1].SetUnitEnergy(unitEnergy); + fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts); } //end Hadron Correction loop - */ - + + + /* //Do Hadron Correction with propagate phi for the track if(fHadCorr != 0) { @@ -197,7 +252,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); phi += deltaPhi; //Get new tower id for cell that track would curve into Int_t towerID2; - if(phi<(1.0/3.0)*TMath::Pi() || phi>TMath::Pi()) + if(phi>(TMath::Pi()/180.0)*geom->GetArm1PhiMax() || phi<(TMath::Pi()/180.0)*geom->GetArm1PhiMin()) { towerID2 = -1; } @@ -218,6 +273,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); fUnitNoCuts[towerID2-1].SetUnitEnergy(unitEnergy2NoCuts); }//end if for towerID2 }//end Hadron Correction loop + */ fUnitNoCuts[towerID-1].SetUnitEnergy(unitEnergyNoCuts + pT); //Do Pt cut on tracks @@ -251,15 +307,16 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); for(Int_t i=0; i10) Info("FillUnitArray","Setting all units outside jets"); - fUnit[i].SetUnitFlag(kOutJet); //Set all units to be outside a jet initially + //Set all units to be outside a jet initially + fUnit[i].SetUnitFlag(kOutJet); fUnit[i].SetUnitID(i+1); Float_t eta; Float_t phi; geom->EtaPhiFromIndex(fUnit[i].GetUnitID(), eta, phi); fUnit[i].SetUnitEta(eta); fUnit[i].SetUnitPhi(phi*TMath::Pi()/180.0); - - fUnitNoCuts[i].SetUnitFlag(kOutJet); //Set all units to be outside a jet initially + //Set all units to be outside a jet initially + fUnitNoCuts[i].SetUnitFlag(kOutJet); fUnitNoCuts[i].SetUnitID(i+1); eta = 0.0; phi = 0.0; @@ -269,6 +326,47 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); // if(i>13000) cout<<"!!!!!!!!!!!!!!!!!For unit0, eta="<1.0e-7) + { + fJetEta = fEtaInit + (fJetEtaSum / fJetESum); + fJetPhi = fPhiInit + (fJetPhiSum / fJetESum); + } } @@ -438,7 +539,14 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); //Loop over all unit objects in the array and find if within cone radius Float_t dEta = fUnit[i].GetUnitEta() - fJetEta; Float_t dPhi = fUnit[i].GetUnitPhi() - fJetPhi; - Float_t rad = TMath::Sqrt( (dEta*dEta) + (dPhi*dPhi) ); + Float_t rad; + if ((dEta*dEta) + (dPhi*dPhi)>1.e-7) + { + rad = TMath::Sqrt( (dEta*dEta) + (dPhi*dPhi) ); + }else + { + rad=0.0; + } if(fUnit[i].GetUnitFlag()==kOutJet && rad<= fConeRad) { @@ -457,7 +565,12 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); { //Stores the resulting jet information in appropriate storage structure (TO BE DECIDED!!!!) if (fDebug>1) Info("StoreJetInfo","Storing Jet Information"); - + AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + AliEMCALGeometry * geom; + if (gime) + geom = gime->EMCALGeometry(); + else + geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","EMCAL"); //Store: //fJetESum is the final jet energy (background has been subtracted) //fJetEta is the final jet Eta @@ -465,9 +578,8 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); //fNumInCone is the final number of cells included in the jet cone //fEtaInit is the eta of the initiator cell //fPhiInit is the phi of the initiator cell - fJet.SetEnergy(fJetESum); - fJet.SetEta(fJetEta); - fJet.SetPhi(fJetPhi); + + AliEMCALJet jet(fJetESum,fJetPhi,fJetEta); cout<<"For iteration "<0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); cout<<"The jet phi is ---->" <GetNTracks(); + Int_t numberDigits = fInputPointer->GetNDigits(); + AliEMCALDigit *myD; TParticle *myP; Int_t numTracksInCone = 0; + Float_t trackEnergy = 0.0; + Float_t trackEnergyPtCut =0.0; + Float_t emcalEnergy = 0.0; + Float_t emcalEnergyBGSub = 0.0; for(Int_t counter=0; counter0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); Float_t phi = myP->Phi(); Float_t deta = fJetEta-eta; Float_t dphi = fJetPhi -phi; - Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi)); + Float_t rad ; + if ((deta*deta) + (dphi*dphi)>1.e-7) + { + rad = TMath::Sqrt( (deta*deta) + (dphi*dphi) ); + }else + { + rad=0.0; + } + if(rad<=fConeRad) { pTArray[index] = myP->Pt(); + //Calculate track contribution within jetcone + trackEnergy += myP->Pt(); + if(myP->Pt() >= fPtCut) trackEnergyPtCut += myP->Pt(); etaArray[index] = eta; phiArray[index] = phi; pdgArray[index] = myP->GetPdgCode(); index++; + if(fHadCorr != 0) + { + Double_t fullP = myP->P(); + Double_t hCEnergy = fHadCorr->GetEnergy(fullP, (Double_t)eta); + emcalEnergy -= hCEnergy*TMath::Sin(myP->Theta()); + emcalEnergyBGSub -= hCEnergy*TMath::Sin(myP->Theta()); + } //end Hadron Correction loop + }//end if }//end for - fJet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray); - fOutputObject.AddJet(&fJet); + //Loop over digits to find EMCal contribution within jetcone + for(Int_t counter3=0; counter3GetDigit(counter3); + //GET DIGIT ETA, PHI so that can check if inside R! + Float_t eta = 0.0; + Float_t phi = 0.0; + Int_t iID = myD->GetId(); + geom->EtaPhiFromIndex(iID, eta, phi); + Float_t deta = fJetEta-eta; + Float_t dphi = fJetPhi -(TMath::Pi()/180.0)*phi; + //Float_t rad = TMath::Sqrt( (deta*deta) + (dphi*dphi)); + Float_t rad ; + if ((deta*deta) + (dphi*dphi)>1.e-7) + { + rad = TMath::Sqrt( (deta*deta) + (dphi*dphi) ); + }else + { + rad=0.0; + } + + if(rad<=fConeRad) + { + Int_t amplitude = myD->GetAmp(); //Gets the integer valued amplitude of the digit + Float_t amp = (Float_t)amplitude; //Need to typecast to Float_t before doing real energy conversion + Float_t digitEnergy = amp/10000000.0; //Factor of 10 million needed to convert! + emcalEnergy += digitEnergy; + emcalEnergyBGSub += (digitEnergy - fEBGAve); + }//end if + }//end count3 for + + //Save in JET object + jet.SetTrackList(numTracksInCone,pTArray, etaArray, phiArray, pdgArray); + jet.SetEMCALEnergy(emcalEnergy); + jet.SetEMCALEnergyBGSub(emcalEnergyBGSub); + jet.SetTrackEnergy(trackEnergy); + jet.SetTrackEnergyPtCut(trackEnergyPtCut); + fOutputPointer->AddJet(&jet); delete[] pTArray; delete[] etaArray; delete[] phiArray; @@ -554,7 +727,7 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); //Step 4. Find the value of the average background energy FindBG(); - fOutputObject.Reset(kResetJets); //Reset output object to store info for new iteration + fOutputPointer->Reset(kResetJets); //Reset output object to store info for new iteration fNumJets=0; //Loop over the array of unit objects and flag those with energy below MinCellEt @@ -599,22 +772,33 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); fPhiInit = fJetPhi; fEtaB = fJetEta; fPhiB = fJetPhi; + Int_t testflag = 1; + do + { fJetESum = 0.0; fJetEtaSum = 0.0; fJetPhiSum = 0.0; //Step 6. Find Jet Eta and Phi //Loop over all units in the array to find the ones in the jet cone and determine contrib to Jet eta, phi - do - { for(Int_t count1=0; count11.e-7) + { + fRad = TMath::Sqrt( (fDEta*fDEta) + (fDPhi*fDPhi) ); + }else + { + fRad=0.000; + } if(fRad <= fConeRad) { FindJetEtaPhi(count1); @@ -624,14 +808,26 @@ if (fDebug>0) Info("AliEMCALJetFinderAlgoOmni","Beginning Default Constructor"); //Find the distance cone centre moved from previous cone centre if (fDebug>10) Info("FindJets","Checking if cone move small enough"); - fDistP = TMath::Sqrt( ((fJetEta-fEtaB)*(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)*(fJetPhi-fPhiB)) ); + if (((fJetEta-fEtaB)*(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)*(fJetPhi-fPhiB)) >1.e-7) + { + fDistP = TMath::Sqrt( ((fJetEta-fEtaB)*(fJetEta-fEtaB)) + ((fJetPhi-fPhiB)*(fJetPhi-fPhiB)) ); + }else + { + fDistP = 0.00; + } // if(fDistP <= fMinMove) break; //Find the distance cone centre is from initiator cell if (fDebug>10) Info("FindJets","Checking if cone move too large"); - fDistI = TMath::Sqrt( ((fJetEtaSum/fJetESum)*(fJetEtaSum/fJetESum)) + ((fJetPhiSum/fJetESum)* + if ( ((fJetEtaSum)*(fJetEtaSum))+((fJetPhiSum)*(fJetPhiSum)) >1.e-7) + { + fDistI = TMath::Sqrt( ((fJetEtaSum/fJetESum)*(fJetEtaSum/fJetESum)) + ((fJetPhiSum/fJetESum)* (fJetPhiSum/fJetESum))); + }else + { + fDistI = 0.00; + } if(fDistP>fMinMove && fDistI