// --- ROOT system ---
#include "TFile.h"
+#include "TRandom3.h"
//---- ANALYSIS system ----
-#include "AliCaloTrackReader.h"
#include "AliMCEvent.h"
#include "AliAODMCHeader.h"
#include "AliGenPythiaEventHeader.h"
-#include "AliVEvent.h"
+#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliVTrack.h"
#include "AliVParticle.h"
#include "AliMixedEvent.h"
#include "AliESDtrack.h"
-#include "AliEMCALRecoUtils.h"
#include "AliESDtrackCuts.h"
+#include "AliTriggerAnalysis.h"
+#include "AliESDVZERO.h"
+
+//---- PartCorr/EMCAL ---
+#include "AliEMCALRecoUtils.h"
+#include "AliCaloTrackReader.h"
ClassImp(AliCaloTrackReader)
//____________________________________________________________________________
AliCaloTrackReader::AliCaloTrackReader() :
- TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0),
+ TObject(), fEventNumber(-1), //fCurrentFileName(""),
+ fDataType(0), fDebug(0),
fFiducialCut(0x0), fCheckFidCut(kFALSE), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0), fAODBranchList(new TList ),
- fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
+ fCTSTracks(new TObjArray()), fEMCALClusters(new TObjArray()), fPHOSClusters(new TObjArray()),
fEMCALCells(0x0), fPHOSCells(0x0),
fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
- fFillEMCALCells(0),fFillPHOSCells(0),
+ fFillEMCALCells(0),fFillPHOSCells(0),
+ fRemoveSuspiciousClusters(kFALSE), fSmearClusterEnergy(kFALSE), fRandom(),
// fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
// fSecondInputFileName(""),fSecondInputFirstEvent(0),
-// fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
-// fAODPHOSNormalInputEntries(0),
+// fCTSTracksNormalInputEntries(0), fEMCALClustersNormalInputEntries(0),
+// fPHOSClustersNormalInputEntries(0),
fTrackStatus(0), fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.8),
fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0),
fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL),
fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),fCaloFilterPatch(kFALSE),
- fEMCALClustersListName(""),fZvtxCut(0.)
+ fEMCALClustersListName(""),fZvtxCut(0.),
+ fDoEventSelection(kFALSE), fDoV0ANDEventSelection(kFALSE),
+ fTriggerAnalysis (new AliTriggerAnalysis), fCentralityClass("V0M"),fCentralityOpt(10)
+
{
//Ctor
delete fAODBranchList ;
}
- if(fAODCTS){
- if(fDataType!=kMC)fAODCTS->Clear() ;
- else fAODCTS->Delete() ;
- delete fAODCTS ;
+ if(fCTSTracks){
+ if(fDataType!=kMC)fCTSTracks->Clear() ;
+ else fCTSTracks->Delete() ;
+ delete fCTSTracks ;
}
- if(fAODEMCAL){
- if(fDataType!=kMC)fAODEMCAL->Clear("C") ;
- else fAODEMCAL->Delete() ;
- delete fAODEMCAL ;
+ if(fEMCALClusters){
+ if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
+ else fEMCALClusters->Delete() ;
+ delete fEMCALClusters ;
}
- if(fAODPHOS){
- if(fDataType!=kMC)fAODPHOS->Clear("C") ;
- else fAODPHOS->Delete() ;
- delete fAODPHOS ;
+ if(fPHOSClusters){
+ if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
+ else fPHOSClusters->Delete() ;
+ delete fPHOSClusters ;
}
// if(fEMCALCells){
delete [] fVertex ;
}
- if(fESDtrackCuts) delete fESDtrackCuts;
+ if(fESDtrackCuts) delete fESDtrackCuts;
+ if(fTriggerAnalysis) delete fTriggerAnalysis;
// Pointers not owned, done by the analysis frame
// if(fInputEvent) delete fInputEvent ;
//Compare jet pT and pt Hard
//if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
- printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
- nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
- return kFALSE;
+ printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
+ nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
+ return kFALSE;
}
}
if(jet) delete jet;
fZvtxCut = 10.;
+ //Centrality
+ fCentralityBin[0]=fCentralityBin[1]=-1;
+
+ //Cluster smearing
+ fSmearClusterEnergy = kFALSE;
+ fSmearClusterParam[0] = 0.07; // * sqrt E term
+ fSmearClusterParam[1] = 0.02; // * E term
+ fSmearClusterParam[2] = 0.00; // constant
+
}
//________________________________________________________________
printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
+ printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
+
printf(" \n") ;
+
}
//___________________________________________________
-Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
+Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * /*currentFileName*/) {
//Fill the event counter and input lists that are needed, called by the analysis maker.
fEventNumber = iEntry;
- fCurrentFileName = TString(currentFileName);
+ //fCurrentFileName = TString(currentFileName);
if(!fInputEvent) {
if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
return kFALSE;
// }
//Fill Vertex array
-
FillVertexArray();
//Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
+ //------------------------------------------------------
+ //Event rejection depending on vertex, pileup, v0and
+ //------------------------------------------------------
+ if(fDoEventSelection){
+ if(!fCaloFilterPatch){
+ //Do not analyze events with pileup
+ Bool_t bPileup = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
+ //Bool_t bPileup = event->IsPileupFromSPD();
+ if(bPileup) return kFALSE;
+
+ if(fDoV0ANDEventSelection){
+ Bool_t bV0AND = kTRUE;
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
+ if(esd)
+ bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
+ //else bV0AND = //FIXME FOR AODs
+ if(!bV0AND) return kFALSE;
+ }
+
+ if(!CheckForPrimaryVertex()) return kFALSE;
+ }//CaloFilter patch
+ else{
+ if(fInputEvent->GetNumberOfCaloClusters() > 0) {
+ AliVCluster * calo = fInputEvent->GetCaloCluster(0);
+ if(calo->GetNLabels() == 4){
+ Int_t * selection = calo->GetLabels();
+ Bool_t bPileup = selection[0];
+ if(bPileup) return kFALSE;
+
+ Bool_t bGoodV = selection[1];
+ if(!bGoodV) return kFALSE;
+
+ if(fDoV0ANDEventSelection){
+ Bool_t bV0AND = selection[2];
+ if(!bV0AND) return kFALSE;
+ }
+
+ fTrackMult = selection[3];
+ if(fTrackMult == 0) return kFALSE;
+ } else {
+ //First filtered AODs, track multiplicity stored there.
+ fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
+ if(fTrackMult == 0) return kFALSE;
+ }
+ }//at least one cluster
+ else {
+ //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
+ //Remove events with vertex (0,0,0), bad vertex reconstruction
+ if(TMath::Abs(fVertex[0][0]) < 1.e-6 && TMath::Abs(fVertex[0][1]) < 1.e-6 && TMath::Abs(fVertex[0][2]) < 1.e-6) return kFALSE;
+
+ //First filtered AODs, track multiplicity stored there.
+ fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
+ if(fTrackMult == 0) return kFALSE;
+ }// no cluster
+ }// CaloFileter patch
+ }// Event selection
+ //------------------------------------------------------
+
+ //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
+ //If we need a centrality bin, we select only those events in the corresponding bin.
+ if(GetCentrality() && fCentralityBin[0]>=0 && fCentralityBin[1]>=0 && fCentralityOpt==100){
+ Int_t cen = GetEventCentrality();
+ if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
+ }
+
//Fill the arrays with cluster/tracks/cells data
if(fFillEMCALCells)
FillInputEMCALCells();
if(fFillCTS){
FillInputCTS();
//Accept events with at least one track
- if(fTrackMult == 0) return kFALSE;
- }
-
- //In case of data produced with calo filter, some information stored in non usual places
- if(IsCaloFilterPatchOn()){
- fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
- //printf("Track multiplicity %d \n",fTrackMult);
- if(fTrackMult == 0) return kFALSE;
+ if(fTrackMult == 0 && fDoEventSelection) return kFALSE;
}
if(fFillEMCAL)
void AliCaloTrackReader::ResetLists() {
// Reset lists, called by the analysis maker
- if(fAODCTS) fAODCTS -> Clear();
- if(fAODEMCAL) fAODEMCAL -> Clear("C");
- if(fAODPHOS) fAODPHOS -> Clear("C");
+ if(fCTSTracks) fCTSTracks -> Clear();
+ if(fEMCALClusters) fEMCALClusters -> Clear("C");
+ if(fPHOSClusters) fPHOSClusters -> Clear("C");
// if(fEMCALCells) fEMCALCells -> Clear("");
// if(fPHOSCells) fPHOSCells -> Clear("");
}
}
+//__________________________________________________
+Int_t AliCaloTrackReader::GetEventCentrality() const {
+ //Return current event centrality
+
+ if(GetCentrality()){
+ if(fCentralityOpt==100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass);
+ else if(fCentralityOpt==10) return GetCentrality()->GetCentralityClass10(fCentralityClass);
+ else if(fCentralityOpt==5) return GetCentrality()->GetCentralityClass5(fCentralityClass);
+ else {
+ printf("AliAnaPartCorrBaseClass::Unknown centrality option %d, use 5, 10 or 100\n",fCentralityOpt);
+ return 0;
+ }
+ }
+ else return 0;
+
+}
+
//____________________________________________________________________________
void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const {
//Return vertex position to be used for single event analysis
track->SetID(itrack);
}
- fAODCTS->Add(track);
+ fCTSTracks->Add(track);
}//Pt and Fiducial cut passed.
}// track loop
- //fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
+ //fCTSTracksNormalInputEntries = fCTSTracks->GetEntriesFast();
if(fDebug > 1)
- printf("AliCaloTrackReader::FillInputCTS() - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fAODCTS->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fAODCTSNormalInputEntries);
+ printf("AliCaloTrackReader::FillInputCTS() - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fCTSTracksNormalInputEntries);
// //If second input event available, add the clusters.
// if(fSecondInputAODTree && fSecondInputAODEvent){
// if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
// momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
//
- // fAODCTS->Add(track);
+ // fCTSTracks->Add(track);
//
// }//Pt and Fiducial cut passed.
// }// track loop
//
- // if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS() - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
+ // if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS() - aod normal entries %d, after second input %d\n", fCTSTracksNormalInputEntries, fCTSTracks->GetEntriesFast());
// } //second input loop
//
}
if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex))
return;
+ //Remove suspicious clusters
+ if(fRemoveSuspiciousClusters){
+ Int_t ncells = clus->GetNCells();
+ Float_t energy = clus->E();
+ Float_t minNCells = 1+energy/3;//-x*x*0.0033
+ if(ncells < minNCells) {
+ //if(energy > 2)printf("AliCaloTrackReader::FillInputEMCALAlgorithm() - Remove cluster: e %2.2f, Ncells %d, min Ncells %2.1f\n",energy,ncells,minNCells);
+ return;
+ }
+// else {
+// if(energy > 2)printf("AliCaloTrackReader::FillInputEMCALAlgorithm() - Keep cluster: e %2.2f, Ncells %d, min Ncells %2.1f\n",energy,ncells,minNCells);
+// }
+ }
+
+
TLorentzVector momentum ;
clus->GetMomentum(momentum, fVertex[vindex]);
//printf("Linearity Corrected Energy %f\n",clus->E());
}
+ //In case of MC analysis, to match resolution/calibration in real data
+ if(fSmearClusterEnergy){
+ Float_t energy = clus->E();
+ Float_t rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0]*TMath::Sqrt(energy)+
+ fSmearClusterParam[1]*energy+fSmearClusterParam[2]);
+ clus->SetE(rdmEnergy);
+ if(fDebug > 2) printf("\t Energy %f, smeared %f\n", energy, clus->E());
+ }
+
if (fMixedEvent)
clus->SetID(iclus) ;
- fAODEMCAL->Add(clus);
+ fEMCALClusters->Add(clus);
}
}
}//EMCAL cluster
}// cluster exists
}// cluster loop
+
+ //Recalculate track matching
+ if(fDataType==kESD)GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
+
}//Get the clusters from the input event
else {
TClonesArray * clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
//printf("E %f\n",clus->E());
if (clus) FillInputEMCALAlgorithm(clus, iclus);
else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
+
}// cluster loop
+
+ //Recalculate track matching, not necessary, already done in the reclusterization task
+ //GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent,clusterList);
+
}
-
- //Recalculate track matching
- GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
-
- //fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
- if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n", fAODEMCAL->GetEntriesFast());//fAODEMCALNormalInputEntries);
+
+ //fEMCALClustersNormalInputEntries = fEMCALClusters->GetEntriesFast();
+ if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n", fEMCALClusters->GetEntriesFast());//fEMCALClustersNormalInputEntries);
//If second input event available, add the clusters.
// if(fSecondInputAODTree && fSecondInputAODEvent){
//
// if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
// momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
- // fAODEMCAL->Add(clus);
+ // fEMCALClusters->Add(clus);
// }//Pt and Fiducial cut passed.
// }//EMCAL cluster
// }// cluster exists
// }// cluster loop
//
- // if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
+ // if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fEMCALClustersNormalInputEntries, fEMCALClusters->GetEntriesFast());
//
// } //second input loop
}
clus->SetID(iclus) ;
}
- fAODPHOS->Add(clus);
+ fPHOSClusters->Add(clus);
}//Pt and Fiducial cut passed.
}//PHOS cluster
}//cluster exists
}//esd cluster loop
- //fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
- if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS() - aod entries %d\n", fAODPHOS->GetEntriesFast());//fAODPHOSNormalInputEntries);
+ //fPHOSClustersNormalInputEntries = fPHOSClusters->GetEntriesFast() ;
+ if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS() - aod entries %d\n", fPHOSClusters->GetEntriesFast());//fPHOSClustersNormalInputEntries);
//If second input event available, add the clusters.
// if(fSecondInputAODTree && fSecondInputAODEvent){
//
// if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
// momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
- // fAODPHOS->Add(clus);
+ // fPHOSClusters->Add(clus);
// }//Pt and Fiducial cut passed.
// }//PHOS cluster
// }// cluster exists
// }// cluster loop
- // if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS() - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
+ // if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS() - aod normal entries %d, after second input %d\n", fPHOSClustersNormalInputEntries, fPHOSClusters->GetEntriesFast());
// } //second input loop
}
}
+//____________________________________________________________________________
+void AliCaloTrackReader::FillInputVZERO(){
+ //Fill VZERO information in data member, add all the channels information.
+ AliVVZERO* v0 = fInputEvent->GetVZEROData();
+ //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
+
+ if (v0)
+ {
+ AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
+ for (Int_t i = 0; i < 32; i++)
+ {
+ if(esdV0){//Only available in ESDs
+ fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
+ fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
+ }
+ fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
+ fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
+ }
+ if(fDebug > 0)
+ printf("V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
+ }
+ else
+ {
+ if(fDebug > 0)
+ printf("Cannot retrieve V0 ESD! Run w/ null V0 charges\n ");
+ }
+}
+
//____________________________________________________________________________
Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const {
}
+//____________________________________________________________________________
+Bool_t AliCaloTrackReader::CheckForPrimaryVertex(){
+ //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
+ //Only for ESDs ...
+
+ AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
+ if(!event) return kFALSE;
+
+ if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
+ return kTRUE;
+ }
+
+ if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
+ // SPD vertex
+ if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
+ //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
+ return kTRUE;
+
+ }
+ if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
+ // cout<<"bad vertex type::"<< event->GetPrimaryVertex()->GetName() << endl;
+ return kFALSE;
+ }
+ }
+
+ return kFALSE;
+
+}
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff