fOutputContainer(new TList ), fAnalysisContainer(new TList ),
fMakeHisto(kFALSE), fMakeAOD(kFALSE), fAnaDebug(0),
fReader(0), fCaloUtils(0),
-fCuts(new TList), fhNEvents(0x0)
+fCuts(new TList), fhNEvents(0x0), fhTrackMult(0x0)
{
//Default Ctor
if(fAnaDebug > 1 ) printf("*** Analysis Maker Constructor *** \n");
fAnaDebug(maker.fAnaDebug),
fReader(),//new AliCaloTrackReader(*maker.fReader)),
fCaloUtils(),//(new AliCalorimeterUtils(*maker.fCaloUtils)),
-fCuts(new TList()), fhNEvents(maker.fhNEvents)
+fCuts(new TList()),fhNEvents(maker.fhNEvents) ,fhTrackMult(maker.fhTrackMult)
{
// cpy ctor
}// Analysis with histograms as output on
}//Loop on analysis defined
}//Analysis list available
+
+ //General event histograms
+
fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
fOutputContainer->Add(fhNEvents);
-
+ fhTrackMult = new TH1I("hTrackMult", "Number of tracks per events" , 2000 , 0 , 2000 ) ;
+ fOutputContainer->Add(fhTrackMult);
+
return fOutputContainer;
}
}
fhNEvents->Fill(0); //Event analyzed
-
+ fhTrackMult->Fill(fReader->GetTrackMultiplicity());
+
fReader->ResetLists();
//printf(">>>>>>>>>> AFTER >>>>>>>>>>>\n");
TList * fCuts ; //! List with analysis cuts
- TH1I * fhNEvents; //! Number of events counter histogram
-
- ClassDef(AliAnaPartCorrMaker,7)
+ TH1I * fhNEvents; //! Number of events counter histogram
+ TH1I * fhTrackMult; //! Number of tracks per event histogram
+
+ ClassDef(AliAnaPartCorrMaker,8)
} ;
fEventNumber = iEntry;
fCurrentFileName = TString(currentFileName);
-
+ fTrackMult = 0;
//In case of analysis of events with jets, skip those with jet pt > 5 pt hard
if(fComparePtHardAndJetPt && GetStack()) {
- if(!ComparePtHardAndJetPt()) return kFALSE ;
+ if(!ComparePtHardAndJetPt()) return kFALSE ;
}
//Fill Vertex array
particle->Momentum(momentum);
//---------- Charged particles ----------------------
if(charge != 0){
- if(fFillCTS && (momentum.Pt() > fCTSPtMin)){
- //Particles in CTS acceptance
-
- if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
-
- if(TMath::Abs(pdg) == 11 && GetStack()->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
-
- if(fDebug > 3 && momentum.Pt() > 0.2)
- printf("AliCaloTrackMCReader::FillInputEvent() - CTS : 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());
-
- x[0] = particle->Vx(); x[1] = particle->Vy(); x[2] = particle->Vz();
- p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
- //Create object and write it to file
- AliAODTrack *aodTrack = new AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0,
- NULL,
- 0x0,//primary,
- kFALSE, // No fit performed
- kFALSE, // No fit performed
- AliAODTrack::kPrimary,
- 0);
- SetTrackChargeAndPID(pdg, aodTrack);
- fAODCTS->Add(aodTrack);//reference the selected object to the list
- }
- //Keep some charged particles in calorimeters lists
- if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg)) FillCalorimeters(iParticle, particle, momentum);
-
+
+ if(TMath::Abs(momentum.Eta())< fTrackMultEtaCut) fTrackMult++;
+
+ if(fFillCTS && (momentum.Pt() > fCTSPtMin)){
+ //Particles in CTS acceptance
+
+ if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
+
+ if(TMath::Abs(pdg) == 11 && GetStack()->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
+
+ if(fDebug > 3 && momentum.Pt() > 0.2)
+ printf("AliCaloTrackMCReader::FillInputEvent() - CTS : 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());
+
+ x[0] = particle->Vx(); x[1] = particle->Vy(); x[2] = particle->Vz();
+ p[0] = particle->Px(); p[1] = particle->Py(); p[2] = particle->Pz();
+ //Create object and write it to file
+ AliAODTrack *aodTrack = new AliAODTrack(0, iParticle, p, kTRUE, x, kFALSE,NULL, 0, 0,
+ NULL,
+ 0x0,//primary,
+ kFALSE, // No fit performed
+ kFALSE, // No fit performed
+ AliAODTrack::kPrimary,
+ 0);
+ SetTrackChargeAndPID(pdg, aodTrack);
+ fAODCTS->Add(aodTrack);//reference the selected object to the list
+ }
+ //Keep some charged particles in calorimeters lists
+ if((fFillPHOS || fFillEMCAL) && KeepChargedParticles(pdg)) FillCalorimeters(iParticle, particle, momentum);
+
}//Charged
//-------------Neutral particles ----------------------
else if(charge == 0 && (fFillPHOS || fFillEMCAL)){
- //Skip neutrinos or other neutral particles
- //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
- if(SkipNeutralParticles(pdg)) continue ;
- //Fill particle/calocluster arrays
- if(!fDecayPi0) {
- FillCalorimeters(iParticle, particle, momentum);
- }
- else {
- //Sometimes pi0 are stable for the generator, if needed decay it by hand
- if(pdg == 111 ){
- if(momentum.Pt() > fPHOSPtMin || momentum.Pt() > fEMCALPtMin){
- TLorentzVector lvGamma1, lvGamma2 ;
- //Double_t angle = 0;
-
- //Decay
- MakePi0Decay(momentum,lvGamma1,lvGamma2);//,angle);
-
- //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
- TParticle * pPhoton1 = new TParticle(22,1,iParticle,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
- lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);
- TParticle * pPhoton2 = new TParticle(22,1,iParticle,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
- lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
- //Fill particle/calocluster arrays
- FillCalorimeters(iParticle,pPhoton1,lvGamma1);
- FillCalorimeters(iParticle,pPhoton2,lvGamma2);
- }//pt cut
- }//pi0
- else FillCalorimeters(iParticle,particle, momentum); //Add the rest
- }
+ //Skip neutrinos or other neutral particles
+ //if(SkipNeutralParticles(pdg) || particle->IsPrimary()) continue ; // protection added (MG)
+ if(SkipNeutralParticles(pdg)) continue ;
+ //Fill particle/calocluster arrays
+ if(!fDecayPi0) {
+ FillCalorimeters(iParticle, particle, momentum);
+ }
+ else {
+ //Sometimes pi0 are stable for the generator, if needed decay it by hand
+ if(pdg == 111 ){
+ if(momentum.Pt() > fPHOSPtMin || momentum.Pt() > fEMCALPtMin){
+ TLorentzVector lvGamma1, lvGamma2 ;
+ //Double_t angle = 0;
+
+ //Decay
+ MakePi0Decay(momentum,lvGamma1,lvGamma2);//,angle);
+
+ //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
+ TParticle * pPhoton1 = new TParticle(22,1,iParticle,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
+ lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);
+ TParticle * pPhoton2 = new TParticle(22,1,iParticle,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
+ lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
+ //Fill particle/calocluster arrays
+ FillCalorimeters(iParticle,pPhoton1,lvGamma1);
+ FillCalorimeters(iParticle,pPhoton2,lvGamma2);
+ }//pt cut
+ }//pi0
+ else FillCalorimeters(iParticle,particle, momentum); //Add the rest
+ }
}//neutral particles
} //particle with correct status
}//particle loop
fIndex2ndPhoton = -1; //In case of overlapping studies, reset for each event
return kTRUE;
-
+
}
//________________________________________________________________
// fSecondInputFileName(""),fSecondInputFirstEvent(0),
// fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0),
// fAODPHOSNormalInputEntries(0),
- fTrackStatus(0), fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.9),
+ fTrackStatus(0), fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.8),
fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0),
//_________________________________________________________________________
Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
- // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
- // Only for PYTHIA.
- if(!fReadStack) return kTRUE; //Information not filtered to AOD
-
- if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
- TParticle * jet = 0;
- AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
- Int_t nTriggerJets = pygeh->NTriggerJets();
- Float_t ptHard = pygeh->GetPtHard();
-
- //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
+ // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
+ // Only for PYTHIA.
+ if(!fReadStack) return kTRUE; //Information not filtered to AOD
+
+ if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
+ TParticle * jet = 0;
+ AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
+ Int_t nTriggerJets = pygeh->NTriggerJets();
+ Float_t ptHard = pygeh->GetPtHard();
+
+ //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
Float_t tmpjet[]={0,0,0,0};
- for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
- pygeh->TriggerJet(ijet, tmpjet);
- jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
- //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;
- }
- }
+ for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
+ pygeh->TriggerJet(ijet, tmpjet);
+ jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
+ //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;
+ }
+ }
if(jet) delete jet;
- }
-
- return kTRUE ;
-
+ }
+
+ return kTRUE ;
+
}
//____________________________________________________________________________
//____________________________________________________________________________
AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
- //Return MC header in AOD. Do it for the corresponding input event.
+ //Return MC header in AOD. Do it for the corresponding input event.
AliAODMCHeader *mch = NULL;
- if(fDataType == kAOD){
- //Normal input AOD
- if(input == 0) {
+ if(fDataType == kAOD){
+ //Normal input AOD
+ if(input == 0) {
mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
}
-// //Second input AOD
-// else if(input == 1){
-// mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
-// }
- else {
- printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
- }
- }
- else {
- printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
- }
+ // //Second input AOD
+ // else if(input == 1){
+ // mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
+ // }
+ else {
+ printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
+ }
+ }
+ else {
+ printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
+ }
return mch;
}
//_______________________________________________________________
void AliCaloTrackReader::Init()
{
- //Init reader. Method to be called in AliAnaPartCorrMaker
-
- //Get the file with second input events if the filename is given
- //Get the tree and connect the AODEvent. Only with AODs
-
- if(fReadStack && fReadAODMCParticles){
- printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
- fReadStack = kFALSE;
- fReadAODMCParticles = kFALSE;
- }
-
+ //Init reader. Method to be called in AliAnaPartCorrMaker
+
+ //Get the file with second input events if the filename is given
+ //Get the tree and connect the AODEvent. Only with AODs
+
+ if(fReadStack && fReadAODMCParticles){
+ printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
+ fReadStack = kFALSE;
+ fReadAODMCParticles = kFALSE;
+ }
+
// if(fSecondInputFileName!=""){
// if(fDataType == kAOD){
// TFile * input2 = new TFile(fSecondInputFileName,"read");
if(fFillPHOSCells)
FillInputPHOSCells();
- if(fFillCTS)
+ if(fFillCTS){
FillInputCTS();
+ //Accept events with at least one track
+ if(fTrackMult == 0) return kFALSE;
+ }
if(fFillEMCAL)
FillInputEMCAL();
if(fFillPHOS)