]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Using detector quality flag (taken from ALICE logbook) to decide whether to rpodcue...
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetResponseMaker.cxx
diff --git a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
index 1d1dd1823791196f6a2d9c74539039166dcf6e5c..b8731048df988fe2dc7cac64baeec3d584b56f62 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
+++ b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
@@ -37,18 +37,12 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fDeltaEtaDeltaPhiAxis(0),
   fNEFAxis(0),
   fZAxis(0),
-  fDoJet2Histogram(0),
   fIsJet1Rho(kFALSE),
   fIsJet2Rho(kFALSE),
-  fHistLeadingJets1PtArea(0),
-  fHistLeadingJets1CorrPtArea(0),
-  fHistLeadingJets2PtArea(0),
-  fHistLeadingJets2CorrPtArea(0),
-  fHistLeadingJets2PtAreaAcceptance(0),
-  fHistLeadingJets2CorrPtAreaAcceptance(0),
+  fHistRejectionReason1(0),
+  fHistRejectionReason2(0),
   fHistJets1(0),
   fHistJets2(0),
-  fHistJets2Acceptance(0),
   fHistMatching(0),
   fHistJets1PhiEta(0),
   fHistJets1PtArea(0),
@@ -59,9 +53,6 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fHistJets2PhiEta(0),
   fHistJets2PtArea(0),
   fHistJets2CorrPtArea(0),
-  fHistJets2PhiEtaAcceptance(0),
-  fHistJets2PtAreaAcceptance(0),
-  fHistJets2CorrPtAreaAcceptance(0),
   fHistJets2NEFvsPt(0),
   fHistJets2CEFvsCEFPt(0),
   fHistJets2ZvsPt(0),
@@ -127,18 +118,12 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fDeltaEtaDeltaPhiAxis(0),
   fNEFAxis(0),
   fZAxis(0),
-  fDoJet2Histogram(0),
   fIsJet1Rho(kFALSE),
   fIsJet2Rho(kFALSE),
-  fHistLeadingJets1PtArea(0),
-  fHistLeadingJets1CorrPtArea(0),
-  fHistLeadingJets2PtArea(0),
-  fHistLeadingJets2CorrPtArea(0),
-  fHistLeadingJets2PtAreaAcceptance(0),
-  fHistLeadingJets2CorrPtAreaAcceptance(0),
+  fHistRejectionReason1(0),
+  fHistRejectionReason2(0),
   fHistJets1(0),
   fHistJets2(0),
-  fHistJets2Acceptance(0),
   fHistMatching(0),
   fHistJets1PhiEta(0),
   fHistJets1PtArea(0),
@@ -149,9 +134,6 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fHistJets2PhiEta(0),
   fHistJets2PtArea(0),
   fHistJets2CorrPtArea(0),
-  fHistJets2PhiEtaAcceptance(0),
-  fHistJets2PtAreaAcceptance(0),
-  fHistJets2CorrPtAreaAcceptance(0),
   fHistJets2NEFvsPt(0),
   fHistJets2CEFvsCEFPt(0),
   fHistJets2ZvsPt(0),
@@ -276,27 +258,6 @@ void AliJetResponseMaker::AllocateTH2()
     fOutput->Add(fHistJets2CorrPtArea);
   }
 
-  fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
-  fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
-  fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistJets2PhiEtaAcceptance);
-
-  fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 
-                                       fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
-  fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
-  fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistJets2PtAreaAcceptance);
-
-  if (fIsJet2Rho) {
-    fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 
-                                             fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
-    fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
-    fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistJets2CorrPtAreaAcceptance);
-  }
-
   fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
   fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
@@ -635,10 +596,10 @@ void AliJetResponseMaker::AllocateTHnSparse()
 {
   // Allocate THnSparse histograms.
 
-  TString title[20]= {""};
-  Int_t nbins[20]  = {0};
-  Double_t min[20] = {0.};
-  Double_t max[20] = {0.};
+  TString title[25]= {""};
+  Int_t nbins[25]  = {0};
+  Double_t min[25] = {0.};
+  Double_t max[25] = {0.};
   Int_t dim = 0;
 
   title[dim] = "#phi";
@@ -665,16 +626,26 @@ void AliJetResponseMaker::AllocateTHnSparse()
   max[dim] = 1.5;
   dim++;
 
-  title[dim] = "NEF";
-  nbins[dim] = fNbins/4;
-  min[dim] = 0;
-  max[dim] = 1.2;
-  dim++;
+  if (fNEFAxis) {
+    title[dim] = "NEF";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+  }
 
-  title[dim] = "Z";
-  nbins[dim] = fNbins/4;
+  if (fZAxis) {
+    title[dim] = "Z";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+  }
+
+  title[dim] = "p_{T,particle}^{leading} (GeV/c)";
+  nbins[dim] = 120;
   min[dim] = 0;
-  max[dim] = 1.2;
+  max[dim] = 120;
   dim++;
 
   Int_t dim1 = dim, dim2 = dim;
@@ -708,17 +679,10 @@ void AliJetResponseMaker::AllocateTHnSparse()
     dim2++;
   }
 
-  if (fDoJet2Histogram) {
-    fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
-    for (Int_t i = 0; i < dim2; i++) 
-      fHistJets2->GetAxis(i)->SetTitle(title[i]);
-    fOutput->Add(fHistJets2);
-  }
-
-  fHistJets2Acceptance = new THnSparseD("fHistJets2Acceptance","fHistJets2Acceptance",dim2,nbins,min,max);
+  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
   for (Int_t i = 0; i < dim2; i++) 
-    fHistJets2Acceptance->GetAxis(i)->SetTitle(title[i]);
-  fOutput->Add(fHistJets2Acceptance);
+    fHistJets2->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistJets2);
   
   // Matching
 
@@ -766,6 +730,18 @@ void AliJetResponseMaker::AllocateTHnSparse()
   max[dim] = 1.2;
   dim++;
 
+  title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
+  title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
   if (fDeltaPtAxis) {
     title[dim] = "#deltaA_{jet}";
     nbins[dim] = fNbins/2;
@@ -882,6 +858,20 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
   else fIsJet2Rho = kTRUE;
 
+  fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
+  fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
+  fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+  fHistRejectionReason1->GetZaxis()->SetTitle("counts");
+  SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis());
+  fOutput->Add(fHistRejectionReason1);
+
+  fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
+  fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
+  fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+  fHistRejectionReason2->GetZaxis()->SetTitle("counts");
+  SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis());
+  fOutput->Add(fHistRejectionReason2);
+
   if (fHistoType==0)
     AllocateTH2();
   else 
@@ -891,7 +881,7 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set)
+void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set)
 {
   if (fHistoType==1) {
     THnSparse *histo = 0;
@@ -899,8 +889,6 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt,
       histo = fHistJets1;
     else if (Set==2)
       histo = fHistJets2;
-    else if (Set==3)
-      histo = fHistJets2Acceptance;
 
     if (!histo)
       return;
@@ -920,11 +908,13 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt,
       else if (title=="NEF")
        contents[i] = NEF;
       else if (title=="Z")
-       contents[i] = Z;
+       contents[i] = LeadingPt/Pt;
       else if (title=="p_{T}^{corr}")
        contents[i] = CorrPt;
       else if (title=="p_{T}^{MC}")
        contents[i] = MCPt;
+      else if (title=="p_{T,particle}^{leading} (GeV/c)")
+       contents[i] = LeadingPt;
       else 
        AliWarning(Form("Unable to fill dimension %s!",title.Data()));
     }
@@ -935,7 +925,7 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt,
     if (Set == 1) {
       fHistJets1PtArea->Fill(A, Pt);
       fHistJets1PhiEta->Fill(Eta, Phi);
-      fHistJets1ZvsPt->Fill(Z, Pt);
+      fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
       fHistJets1NEFvsPt->Fill(NEF, Pt);
       fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
       if (fIsJet1Rho)
@@ -944,17 +934,11 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt,
     else if (Set == 2) {
       fHistJets2PtArea->Fill(A, Pt);
       fHistJets2PhiEta->Fill(Eta, Phi);
-      if (fIsJet2Rho)
-       fHistJets2CorrPtArea->Fill(A, CorrPt);
-    }
-    else if (Set == 3) {
-      fHistJets2PtAreaAcceptance->Fill(A, Pt);
-      fHistJets2PhiEtaAcceptance->Fill(Eta, Phi);
-      fHistJets2ZvsPt->Fill(Z, Pt);
+      fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
       fHistJets2NEFvsPt->Fill(NEF, Pt);
       fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
       if (fIsJet2Rho)
-       fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt);
+       fHistJets2CorrPtArea->Fill(A, CorrPt);
     }
   }
 }
@@ -962,7 +946,7 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt,
 //________________________________________________________________________
 void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, 
                                             Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, 
-                                            Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t Z1, Double_t Z2)
+                                            Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
 {
   if (fHistoType==1) {
     Double_t contents[20]={0};
@@ -1006,9 +990,13 @@ void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_
       else if (title=="NEF_{2}")
        contents[i] = NEF2;
       else if (title=="Z_{1}")
-       contents[i] = Z1;
+       contents[i] = LeadingPt1/Pt1;
       else if (title=="Z_{2}")
-       contents[i] = Z2;
+       contents[i] = LeadingPt2/Pt2;
+      else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
+       contents[i] = LeadingPt1;
+      else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
+       contents[i] = LeadingPt2;
       else 
        AliWarning(Form("Unable to fill dimension %s!",title.Data()));
     }
@@ -1222,7 +1210,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
       }
 
       Int_t MClabel = TMath::Abs(track->GetLabel());
-      if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
+      MClabel -= fMCLabelShift;
       if (MClabel != 0) continue;
 
       // this is not a MC particle; remove it completely
@@ -1248,7 +1236,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
        Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
 
        Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
-       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
+       MClabel -= fMCLabelShift;
        if (MClabel != 0) continue;
 
        // this is not a MC particle; remove it completely
@@ -1269,7 +1257,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
          
       Int_t MClabel = TMath::Abs(clus->GetLabel());
-      if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
+      MClabel -= fMCLabelShift;
       if (MClabel != 0) continue;
 
       // this is not a MC particle; remove it completely
@@ -1291,7 +1279,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
        continue;
       }
       Int_t MClabel = TMath::Abs(track->GetLabel());
-      if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;     
+      MClabel -= fMCLabelShift;          
       if (MClabel <= 0) continue;
 
       Int_t index = -1;
@@ -1332,7 +1320,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
          Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
        
          Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
-         if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
+         MClabel -= fMCLabelShift;
          if (MClabel <= 0) continue;
          
          Int_t index1 = -1;
@@ -1369,7 +1357,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
        clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
       
        Int_t MClabel = TMath::Abs(clus->GetLabel());
-       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;      
+       MClabel -= fMCLabelShift;           
        if (MClabel <= 0) continue;
 
        Int_t index = -1;
@@ -1644,15 +1632,19 @@ Bool_t AliJetResponseMaker::FillHistograms()
   
   jets2->ResetCurrentID();
   while ((jet2 = jets2->GetNextJet())) {
+
     AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
 
-    if (fDoJet2Histogram)
-      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), 
-                  jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 2);
+    if (jet2->Pt() < jets2->GetJetPtCut()) continue;
+
+    Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
+    Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
 
     if (jets2->AcceptJet(jet2))
-      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), 
-                  jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 3);
+      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2, 
+                  corrpt2, jet2->MCPt(), 2);
+    else
+      fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt());
 
     jet1 = jet2->MatchedJet();
 
@@ -1676,21 +1668,28 @@ Bool_t AliJetResponseMaker::FillHistograms()
       ce2 = jet2->ClosestJetDistance();
     }
 
+    Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
     Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
-    Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
 
     FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), 
                       jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), 
-                      jet1->NEF(), jet2->NEF(), jet1->MaxPartPt()/jet1->Pt(), jet2->MaxPartPt()/jet2->Pt());
+                      jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2);
   }
 
   jets1->ResetCurrentID();
-  while ((jet1 = jets1->GetNextAcceptJet())) {
+  while ((jet1 = jets1->GetNextJet())) {
+    if (!jets1->AcceptJet(jet1)) {
+      fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt());
+      continue;
+    }
     if (jet1->MCPt() < fMinJetMCPt) continue;
     AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
 
-    FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), jet1->MaxPartPt()/jet1->Pt(), 
-                jet1->Pt() - jets1->GetRhoVal() * jet1->Area(), jet1->MCPt(), 1);
+    Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+    Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
+
+    FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1, 
+                corrpt1, jet1->MCPt(), 1);
   }
   return kTRUE;
 }