update/fixes from Salvatore
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Feb 2013 14:41:04 +0000 (14:41 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Feb 2013 14:41:04 +0000 (14:41 +0000)
13 files changed:
PWG/EMCAL/AliEmcalMCTrackSelector.cxx
PWG/EMCAL/AliEmcalPicoTrackMaker.cxx
PWG/EMCAL/AliEmcalPicoTrackMaker.h
PWG/EMCAL/macros/AddTaskEmcalPicoTrackMaker.C
PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
PWGJE/EMCALJetTasks/AliJetRandomizerTask.cxx
PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
PWGJE/EMCALJetTasks/AliJetResponseMaker.h
PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C
PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C

index 6d2e6d3..d6e15f6 100644 (file)
@@ -93,13 +93,13 @@ void AliEmcalMCTrackSelector::UserExec(Option_t *)
     fInit = kTRUE;
   }
 
-  new (fTracksMap) TH1I(fTracksMapName, fTracksMapName, 1000, 0, 1);
-
   // clear container (normally a null operation as the event should clean it already)
   fTracksOut->Delete();
-  // loop over tracks
+
   const Int_t Ntracks = mcevent->GetNumberOfTracks();
+  new (fTracksMap) TH1I(fTracksMapName, fTracksMapName, Ntracks-2, 0, 1);  // Ntracks - 2, we use also over- and uner-flow bins
+
+  // loop over tracks
   for (Int_t iTracks = 0, nacc = 0; iTracks < Ntracks; ++iTracks) {
 
     fTracksMap->SetBinContent(iTracks, -1);
index 0bb128b..f344601 100644 (file)
@@ -33,6 +33,8 @@ AliEmcalPicoTrackMaker::AliEmcalPicoTrackMaker() :
   fMaxTrackPhi(10),
   fTrackEfficiency(1),
   fIncludeNoITS(kTRUE),
+  fUseNegativeLabels(kFALSE),
+  fIsMC(kFALSE),
   fTracksIn(0),
   fTracksOut(0)
 {
@@ -56,6 +58,8 @@ AliEmcalPicoTrackMaker::AliEmcalPicoTrackMaker(const char *name) :
   fMaxTrackPhi(10),
   fTrackEfficiency(1),
   fIncludeNoITS(kTRUE),
+  fUseNegativeLabels(kFALSE),
+  fIsMC(kFALSE),
   fTracksIn(0),
   fTracksOut(0)
 {
@@ -185,9 +189,14 @@ void AliEmcalPicoTrackMaker::UserExec(Option_t *)
        continue;
     }
 
-    Int_t label = -1;
-    if (track->GetLabel() >= 0)
+    Int_t label = 0;
+    if (track->GetLabel() > 0)
       label = track->GetLabel();
+    else if (fUseNegativeLabels)
+      label = -track->GetLabel();
+
+    if (fIsMC && label == 0) 
+      label = 99999;
 
     /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(track->Pt(), 
                                                                         track->Eta(), 
index 500377b..45d20d8 100644 (file)
@@ -26,8 +26,10 @@ class AliEmcalPicoTrackMaker : public AliAnalysisTaskSE {
   void SetTrackPtLimits(Double_t min, Double_t max)      { fMaxTrackPt       = max ; fMinTrackPt       = min ; }
   void SetTrackEtaLimits(Double_t min, Double_t max)     { fMaxTrackEta      = max ; fMinTrackEta      = min ; }
   void SetTrackPhiLimits(Double_t min, Double_t max)     { fMaxTrackPhi      = max ; fMinTrackPhi      = min ; }
-  void SetTrackEfficiency(Double_t eff = 0.95)           { fTrackEfficiency = eff  ; }
-  void SetIncludeNoITS(Bool_t f)                         { fIncludeNoITS     = f;    }
+  void SetTrackEfficiency(Double_t eff = 0.95)           { fTrackEfficiency  = eff ; }
+  void SetIncludeNoITS(Bool_t f)                         { fIncludeNoITS     = f   ; }
+  void SetUseNegativeLabes(Bool_t f)                     { fUseNegativeLabels= f   ; }
+  void SetMC(Bool_t a)                                   { fIsMC             = a   ; }
 
  protected:
   Int_t              fAODfilterBits[2];     // AOD track filter bit map
@@ -42,6 +44,8 @@ class AliEmcalPicoTrackMaker : public AliAnalysisTaskSE {
   Double_t           fMaxTrackPhi;          // cut on track phi
   Double_t           fTrackEfficiency;      // track efficiency
   Bool_t             fIncludeNoITS;         // includes tracks with failed ITS refit
+  Bool_t             fUseNegativeLabels;    // whether or not should use negative MC labels
+  Bool_t             fIsMC;                 // whether it is a MC event or not
   TClonesArray      *fTracksIn;             //!track array in
   TClonesArray      *fTracksOut;            //!track array out
 
@@ -49,6 +53,6 @@ class AliEmcalPicoTrackMaker : public AliAnalysisTaskSE {
   AliEmcalPicoTrackMaker(const AliEmcalPicoTrackMaker&);            // not implemented
   AliEmcalPicoTrackMaker &operator=(const AliEmcalPicoTrackMaker&); // not implemented
 
-  ClassDef(AliEmcalPicoTrackMaker, 4); // Task to make PicoTracks in AOD/ESD events
+  ClassDef(AliEmcalPicoTrackMaker, 5); // Task to make PicoTracks in AOD/ESD events
 };
 #endif
index 795baf4..02fdb5f 100644 (file)
@@ -37,9 +37,15 @@ AliEmcalPicoTrackMaker* AddTaskEmcalPicoTrackMaker(
   runPeriod.ToLower();
   if (runPeriod == "lhc11h") {
     eTask->SetAODfilterBits(256,512); // hybrid tracks for LHC11h
+    eTask->SetMC(kFALSE);
   }
-  else if (runPeriod == "lhc11a" || runPeriod == "lhc12a15a" || runPeriod == "lhc12a15e") {
-    eTask->SetAODfilterBits(256,16); // hybrid tracks for LHC11a, LHC12a15a and LHC12a15e
+  else if (runPeriod == "lhc11a") {
+    eTask->SetAODfilterBits(256,16); // hybrid tracks for LHC11a
+    eTask->SetMC(kFALSE);
+  }
+  else if (runPeriod == "lhc12a15a" || runPeriod == "lhc12a15e") {
+    eTask->SetAODfilterBits(256,16); // hybrid tracks for LHC12a15a and LHC12a15e
+    eTask->SetMC(kTRUE);
   }
   else {
     if (runPeriod.IsNull())
index c67c684..99ee21f 100644 (file)
@@ -456,7 +456,7 @@ void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
     if (vtrack && !AcceptTrack(vtrack)) 
       continue;
 
-    if (track->GetLabel() >= 0) {
+    if (track->GetLabel() > 0) {
       fEmbeddedTrackIds[fEmbeddedTrackNIds] = i;
       fEmbeddedTrackNIds++;
     }
@@ -485,7 +485,7 @@ void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
     if (!AcceptCluster(cluster)) 
       continue;
 
-    if (cluster->GetLabel() >= 0) {
+    if (cluster->GetLabel() > 0) {
       fEmbeddedClusterIds[fEmbeddedClusterNIds] = iClusters;
       fEmbeddedClusterNIds++;
     }
@@ -501,6 +501,8 @@ AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i)
 
   if (i >= 0)
     iJet = i;
+  else
+    iJet++;
 
   if (!fEmbJets)
     return 0;
@@ -607,7 +609,7 @@ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &p
        AliError(Form("Could not retrieve track %d",iTracks)); 
        continue; 
       }
-      
+
       if (!AcceptTrack(track)) 
        continue;
       
index 0725ee6..0c78aa2 100644 (file)
@@ -290,7 +290,7 @@ void AliEmcalJetTask::FindJets()
             maxCh = cPt;
         }
 
-        if (fIsMcPart || t->GetLabel() >= 0) // check if MC particle
+        if (fIsMcPart || t->GetLabel() > 0) // check if MC particle
           mcpt += cPt;
 
         if (cPhi<0) 
@@ -335,7 +335,7 @@ void AliEmcalJetTask::FindJets()
         if (cPt > maxNe)
           maxNe = cPt;
 
-        if (c->GetLabel() >= 0) // MC particle
+        if (c->GetLabel() > 0) // MC particle
           mcpt += cPt;
 
         if (cPhi<0) 
index ef5ba0d..23c424d 100644 (file)
@@ -62,7 +62,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
 {
   // Default constructor.
   SetSuffix("AODEmbedding");
-  fMarkMC = kFALSE;
+  SetMarkMC(0);
   fAODfilterBits[0] = -1;
   fAODfilterBits[1] = -1;
 }
@@ -97,7 +97,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
 {
   // Standard constructor.
   SetSuffix("AODEmbedding");
-  fMarkMC = kFALSE;
+  SetMarkMC(0);
   fAODfilterBits[0] = -1;
   fAODfilterBits[1] = -1;
 }
@@ -379,7 +379,7 @@ void AliJetEmbeddingFromAODTask::Run()
 
     if (fAODCaloCells) {
       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
-       Short_t mclabel = -1;
+       Short_t mclabel = 0;
        Double_t efrac = 0.;
        Double_t time = -1;
        Short_t cellNum = -1;
index 2377f73..bc201b0 100644 (file)
@@ -48,7 +48,7 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fNClusters(0),
   fNCells(0),
   fNTracks(0),
-  fMarkMC(100),
+  fMarkMC(99999),
   fPtSpectrum(0),
   fQAhistos(kFALSE),
   fIsInit(0),
@@ -87,7 +87,7 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fNClusters(0),
   fNCells(0),
   fNTracks(1),
-  fMarkMC(100),
+  fMarkMC(99999),
   fPtSpectrum(0),
   fQAhistos(drawqa),
   fIsInit(0),
@@ -298,12 +298,16 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     if (!fGeom) {
       if (fGeomName.Length() > 0) {
        fGeom = AliEMCALGeometry::GetInstance(fGeomName);
-       if (!fGeom)
-         AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
+       if (!fGeom) {
+         AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
+         return kFALSE;
+       }
       } else {
        fGeom = AliEMCALGeometry::GetInstance();
-       if (!fGeom) 
-         AliError("Could not get default geometry!");
+       if (!fGeom) {
+         AliFatal("Could not get default geometry!");
+         return kFALSE;
+       }
       }
     }
   
@@ -346,7 +350,7 @@ Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
 void AliJetModelBaseTask::CopyCells()
 {
   for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
-    Short_t mclabel = -1;
+    Short_t mclabel = 0;
     Double_t efrac = 0.;
     Double_t time = -1;
     Short_t cellNum = -1;
index 7305d6e..57230e1 100644 (file)
@@ -29,7 +29,7 @@ AliJetRandomizerTask::AliJetRandomizerTask() :
   // Default constructor.
 
   SetSuffix("Randomized");
-  SetMarkMC(-1);
+  SetMarkMC(0);
 }
 
 //________________________________________________________________________
@@ -40,7 +40,7 @@ AliJetRandomizerTask::AliJetRandomizerTask(const char *name) :
   // Standard constructor.
 
   SetSuffix("Randomized");
-  SetMarkMC(-1);
+  SetMarkMC(0);
 }
 
 //________________________________________________________________________
index f3c5fb3..d06d655 100644 (file)
@@ -34,7 +34,8 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fAreCollections1MC(kFALSE),  
   fAreCollections2MC(kTRUE),
   fMatching(kNoMatching),
-  fMatchingPar(0),
+  fMatchingPar1(0),
+  fMatchingPar2(0),
   fJet2MinEta(-999),
   fJet2MaxEta(-999),
   fJet2MinPhi(-999),
@@ -92,7 +93,8 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fAreCollections1MC(kFALSE),  
   fAreCollections2MC(kTRUE),
   fMatching(kNoMatching),
-  fMatchingPar(0.25),
+  fMatchingPar1(0),
+  fMatchingPar2(0),
   fJet2MinEta(-999),
   fJet2MaxEta(-999),
   fJet2MinPhi(-999),
@@ -390,9 +392,14 @@ void AliJetResponseMaker::ExecOnce()
 
     if (fAreCollections2MC) {
       fTracks2Map = dynamic_cast<TH1*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
+      // this is needed to map the MC labels with the indexes of the MC particle collection
+      // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
       if (!fTracks2Map) {
-       AliError(Form("%s: Could not retrieve map for tracks2 %s!", GetName(), fTracks2Name.Data())); 
-       return;
+       AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data())); 
+       fTracks2Map = new TH1I("tracksMap","tracksMap",9999,0,1);
+       for (Int_t i = 0; i < 9999; i++) {
+         fTracks2Map->SetBinContent(i,i);
+       }
       }
     }
   }
@@ -512,7 +519,7 @@ Bool_t AliJetResponseMaker::DoJetMatching()
       continue;
 
     if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 && 
-        jet1->ClosestJetDistance() < fMatchingPar && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar) {    // Matched jet found
+        jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) {    // Matched jet found
       jet1->SetMatchedToClosest(fMatching);
       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
     }
@@ -583,146 +590,241 @@ void AliJetResponseMaker::DoJetLoop(Bool_t order)
          continue;
       }
 
-      GetMatchingLevel(jet1, jet2, fMatching);
+      SetMatchingLevel(jet1, jet2, fMatching);
     }
   }
 }
 
 //________________________________________________________________________
-Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
+void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
 {
-  Double_t d1 = -1;
-  Double_t d2 = -1;
+  Double_t deta = jet2->Eta() - jet1->Eta();
+  Double_t dphi = jet2->Phi() - jet1->Phi();
+  d = TMath::Sqrt(deta * deta + dphi * dphi);
+}
 
-  switch (matching) {
-  case kGeometrical:
-    {
-      Double_t deta = jet2->Eta() - jet1->Eta();
-      Double_t dphi = jet2->Phi() - jet1->Phi();
-      d1 = TMath::Sqrt(deta * deta + dphi * dphi);
-      d2 = d1;
+//________________________________________________________________________
+void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{ 
+  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
+  d1 = jet1->Pt();
+  d2 = jet2->Pt();
+  Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
+
+  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
+    Bool_t track2Found = kFALSE;
+    Int_t index2 = jet2->TrackAt(iTrack2);
+    for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
+      AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
+      if (!track) {
+       AliWarning(Form("Could not find track %d!", iTrack));
+       continue;
+      }
+      Int_t MClabel = track->GetLabel();
+      Int_t index = -1;
+         
+      if (MClabel <= 0) {// this is not a MC particle; remove it completely
+       AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+       totalPt1 -= track->Pt();
+       d1 -= track->Pt();
+       continue;
+      }
+      else if (MClabel < fTracks2Map->GetNbinsX()-2) {
+       index = fTracks2Map->GetBinContent(MClabel);
+      }
+         
+      if (index < 0) {
+       AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+       continue;
+      }
+
+      if (index2 == index) { // found common particle
+       track2Found = kTRUE;
+       d1 -= track->Pt();
+       AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
+       AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                       iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+       d2 -= MCpart->Pt();
+       break;
+      }
     }
-    break;
-  case kMCLabel: // jet1 = detector level and jet2 = particle level!
-    { 
-      if (!fTracks2Map) {
-       fTracks2Map = new TH1I("tracksMap","tracksMap",1000,0,1);
-       for (Int_t i = 0; i < 1000; i++) {
-         fTracks2Map->SetBinContent(i,i);
+    for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+      AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
+      if (!clus) {
+       AliWarning(Form("Could not find cluster %d!", iClus));
+       continue;
+      }
+      TLorentzVector part;
+      clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+         
+      if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
+       for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
+         Int_t cellId = clus->GetCellAbsId(iCell);
+         Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
+
+         Int_t MClabel = fCaloCells->GetCellMCLabel(cellId);
+         Int_t index = -1;
+         
+         if (MClabel <= 0) {// this is not a MC particle; remove it completely
+           AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+           totalPt1 -= part.Pt() * cellFrac;
+           d1 -= part.Pt() * cellFrac;
+           continue;
+         }
+         else if (MClabel < fTracks2Map->GetNbinsX()-2) {
+           index = fTracks2Map->GetBinContent(MClabel);
+         }
+
+         if (index < 0) {
+           AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+           continue;
+         }
+         if (index2 == index) { // found common particle
+           d1 -= part.Pt() * cellFrac;
+               
+           if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
+             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
+             AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                             iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));             
+             d2 -= MCpart->Pt() * cellFrac;
+           }
+           break;
+         }
        }
       }
-      d1 = jet1->Pt();
-      d2 = jet2->Pt();
-      Double_t totalPt1 = d1;
-      for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
-       AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
-       if (!track) {
-         AliWarning(Form("Could not find track %d!", iTrack));
+      else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
+       Int_t MClabel = clus->GetLabel();
+       Int_t index = -1;
+           
+       if (MClabel <= 0) {// this is not a MC particle; remove it completely
+         AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
+         totalPt1 -= part.Pt();
+         d1 -= part.Pt();
          continue;
        }
-       Int_t MClabel = track->GetLabel();
-       if (MClabel < 0) {// this is not a MC particle; remove it completely
-         AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
-         totalPt1 -= track->Pt();
-         d1 -= track->Pt();
-         continue;
+       else if (MClabel < fTracks2Map->GetNbinsX()-2) {
+         index = fTracks2Map->GetBinContent(MClabel);
        }
-       Int_t index = fTracks2Map->GetBinContent(MClabel);
+        
        if (index < 0) {
-         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+         AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
          continue;
        }
-       for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
-         Int_t index2 = jet2->TrackAt(iTrack2);
-         if (index2 == index) { // found common particle
-           d1 -= track->Pt();
+       if (index2 == index) { // found common particle
+         d1 -= part.Pt();
+
+         if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
            AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
-           AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
-                           iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+           AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                           iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+               
            d2 -= MCpart->Pt();
-           break;
          }
+         break;
        }
       }
-      for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
-       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
-       if (!clus) {
-         AliWarning(Form("Could not find cluster %d!", iClus));
-         continue;
-       }
-       TLorentzVector part;
-       clus->GetMomentum(part, fVertex);
-
-       if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
-         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
-           Int_t cellId = clus->GetCellAbsId(iCell);
-           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
-           Int_t MClabel = fCaloCells->GetCellMCLabel(cellId);
-
-           if (MClabel < 0) {// this is not a MC particle; remove it competely
-             AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
-             totalPt1 -= part.Pt() * cellFrac;
-             d1 -= part.Pt() * cellFrac;
-             continue;
-           }
+    }
+  }
+  if (d1 <= 0 || totalPt1 < 1)
+    d1 = 0;
+  else
+    d1 /= totalPt1;
 
-           Int_t index = fTracks2Map->GetBinContent(MClabel);
-           if (index < 0) {
-             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
-             continue;
-           }
-           for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
-             Int_t index2 = jet2->TrackAt(iTrack2);
-             if (index2 == index) { // found common particle
-               d1 -= part.Pt() * cellFrac;
-               AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
-               AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
-                               iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
-               if (MCpart->Charge() != 0) // only if it is a neutral particle (charged particles are most likely already removed, to be fixed)
-                 d2 -= MCpart->Pt() * cellFrac;
-               break;
-             }
-           }
-         }
-       }
-       else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
-         Int_t MClabel = clus->GetLabel();
-         
-         if (MClabel < 0) {// this is not a MC particle; remove it competely
-           AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
-           totalPt1 -= part.Pt();
-           d1 -= part.Pt();
+  if (jet2->Pt() > 0 && d2 > 0)
+    d2 /= jet2->Pt();
+  else
+    d2 = 0;
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{ 
+  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
+  d1 = jet1->Pt();
+  d2 = jet2->Pt();
+
+  if (fTracks && fTracks2) {
+
+    for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
+      Int_t index2 = jet2->TrackAt(iTrack2);
+      for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
+       Int_t index = jet1->TrackAt(iTrack);
+       if (index2 == index) { // found common particle
+         AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
+         if (!part) {
+           AliWarning(Form("Could not find track %d!", index));
            continue;
          }
-         
-         Int_t index = fTracks2Map->GetBinContent(MClabel);
-         if (index < 0) {
-           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
+         AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
+         if (!part2) {
+           AliWarning(Form("Could not find track %d!", index2));
            continue;
          }
-         for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
-           Int_t index2 = jet2->TrackAt(iTrack2);
-           if (index2 == index) { // found common particle
-             d1 -= part.Pt();
-             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
-             AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
-                             iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
-             if (MCpart->Charge() != 0) // only if it is a neutral particle (charged particles are most likely already removed, to be fixed)
-               d2 -= MCpart->Pt();
-             break;
-           }
-         }
+
+         d1 -= part->Pt();
+         d2 -= part2->Pt();
+         break;
        }
       }
-      if (d1 <= 0 || totalPt1 < 1)
-       d1 = 0;
-      else
-       d1 /= totalPt1;
-      if (jet2->Pt() > 0)
-       d2 /= jet2->Pt();
-      else
-       d2 = 0;
     }
+
+  }
+
+  if (fCaloClusters && fCaloClusters2) {
+
+    for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
+      Int_t index2 = jet2->ClusterAt(iClus2);
+      for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+       Int_t index = jet1->ClusterAt(iClus);
+       AliVCluster *clus =  static_cast<AliVCluster*>(fCaloClusters->At(index));
+       if (!clus) {
+         AliWarning(Form("Could not find cluster %d!", index));
+         continue;
+       }
+       AliVCluster *clus2 =  static_cast<AliVCluster*>(fCaloClusters2->At(index2));
+       if (!clus2) {
+         AliWarning(Form("Could not find cluster %d!", index2));
+         continue;
+       }
+       TLorentzVector part, part2;
+       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+       clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+
+       d1 -= part.Pt();
+       d2 -= part2.Pt();
+       break;
+      }
+    }
+
+  }
+
+  if (jet1->Pt() > 0 && d1 > 0)
+    d1 /= jet1->Pt();
+  else
+    d1 = 0;
+
+  if (jet2->Pt() > 0 && d2 > 0)
+    d2 /= jet2->Pt();
+  else
+    d2 = 0;
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
+{
+  Double_t d1 = -1;
+  Double_t d2 = -1;
+
+  switch (matching) {
+  case kGeometrical:
+    GetGeometricalMatchingLevel(jet1,jet2,d1);
+    d2 = d1;
+    break;
+  case kMCLabel: // jet1 = detector level and jet2 = particle level!
+    GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
+    break;
+  case kSameCollections:
+    GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
     break;
   default:
     ;
@@ -749,8 +851,6 @@ Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *j
       jet2->SetSecondClosestJet(jet1, d2);
     }
   }
-  
-  return d1;
 }
 
 //________________________________________________________________________
@@ -808,12 +908,20 @@ Bool_t AliJetResponseMaker::FillHistograms()
        fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
       }
       else {
-       if (jet2->GetMatchingType() == kGeometrical)
-         fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
-       else if (jet2->GetMatchingType() == kMCLabel)
-         fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), jet2->ClosestJetDistance());
-       else
-         fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
+       Double_t d1=-1, d2=-1;
+       if (jet2->GetMatchingType() == kGeometrical) {
+
+         if (fAreCollections2MC && !fAreCollections1MC)
+           GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
+         else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
+           GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
+
+         fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), d2);
+       }
+       else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
+         GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
+         fHistDistancevsCommonEnergy->Fill(d1, jet2->ClosestJetDistance());
+       }
          
        fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
 
index 0165e54..c5a4e53 100644 (file)
@@ -21,23 +21,23 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   enum MatchingType{
     kNoMatching = 0,
     kGeometrical = 1,
-    kMCLabel = 2
+    kMCLabel = 2,
+    kSameCollections = 3
   };
 
   void                        UserCreateOutputObjects();
 
-  void                        SetJets2Name(const char *n)                            { fJets2Name         = n         ; }
-  void                        SetTracks2Name(const char *n)                          { fTracks2Name       = n         ; }
-  void                        SetClus2Name(const char *n)                            { fCalo2Name         = n         ; }
-  void                        SetJet2EtaLimits(Float_t min=-999, Float_t max=-999)   { fJet2MinEta = min, fJet2MaxEta = max ; }
-  void                        SetJet2PhiLimits(Float_t min=-999, Float_t max=-999)   { fJet2MinPhi = min, fJet2MaxPhi = max ; }
-  void                        SetRho2Name(const char *n)                             { fRho2Name          = n         ; }
-  void                        SetPtBiasJet2Clus(Float_t b)                           { fPtBiasJet2Clus    = b         ; }
-  void                        SetPtBiasJet2Track(Float_t b)                          { fPtBiasJet2Track   = b         ; }
-  void                        SetMatching(MatchingType t, Double_t p=1)              { fMatching = t; fMatchingPar = p; }
-  void                        SetPtHardBin(Int_t b)                                  { fSelectPtHardBin   = b         ; }
-  void                        SetMCflag1(Bool_t f)                                   { fAreCollections1MC = f         ; }
-  void                        SetMCflag2(Bool_t f)                                   { fAreCollections2MC = f         ; }
+  void                        SetJets2Name(const char *n)                                     { fJets2Name         = n         ; }
+  void                        SetTracks2Name(const char *n)                                   { fTracks2Name       = n         ; }
+  void                        SetClus2Name(const char *n)                                     { fCalo2Name         = n         ; }
+  void                        SetJet2EtaLimits(Float_t min=-999, Float_t max=-999)            { fJet2MinEta = min, fJet2MaxEta = max ; }
+  void                        SetJet2PhiLimits(Float_t min=-999, Float_t max=-999)            { fJet2MinPhi = min, fJet2MaxPhi = max ; }
+  void                        SetRho2Name(const char *n)                                      { fRho2Name          = n         ; }
+  void                        SetPtBiasJet2Clus(Float_t b)                                    { fPtBiasJet2Clus    = b         ; }
+  void                        SetPtBiasJet2Track(Float_t b)                                   { fPtBiasJet2Track   = b         ; }
+  void                        SetMatching(MatchingType t, Double_t p1=1, Double_t p2=1)       { fMatching = t; fMatchingPar1 = p1; fMatchingPar2 = p2; }
+  void                        SetPtHardBin(Int_t b)                                           { fSelectPtHardBin   = b         ; }
+  void                        SetAreMCCollections(Bool_t f1, Bool_t f2)                       { fAreCollections1MC = f1; fAreCollections2MC = f2; }
 
  protected:
   Bool_t                      IsEventSelected();
@@ -49,7 +49,10 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   Bool_t                      RetrieveEventObjects();
   Bool_t                      Run();
   Bool_t                      DoJetMatching();
-  Double_t                    GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching);
+  void                        SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching);
+  void                        GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const;
+  void                        GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const;
+  void                        GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const;
 
   TString                     fTracks2Name;                   // name of second track collection
   TString                     fCalo2Name;                     // name of second cluster collection
@@ -60,7 +63,8 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   Bool_t                      fAreCollections1MC;             // collections 1 MC
   Bool_t                      fAreCollections2MC;             // collections 1 MC
   MatchingType                fMatching;                      // matching type
-  Double_t                    fMatchingPar;                   // maximum distance between matched particle and detector level jets
+  Double_t                    fMatchingPar1;                  // matching parameter for generated-reconstructed matching
+  Double_t                    fMatchingPar2;                  // matching parameter for reconstructed-generated matching
   Float_t                     fJet2MinEta;                    // minimum eta jet 2 acceptance
   Float_t                     fJet2MaxEta;                    // maximum eta jet 2 acceptance
   Float_t                     fJet2MinPhi;                    // minimum phi jet 2 acceptance
@@ -110,6 +114,6 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   AliJetResponseMaker(const AliJetResponseMaker&);            // not implemented
   AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented
 
-  ClassDef(AliJetResponseMaker, 10) // Jet response matrix producing task
+  ClassDef(AliJetResponseMaker, 11) // Jet response matrix producing task
 };
 #endif
index 8f3add9..5048946 100644 (file)
@@ -15,7 +15,8 @@ AliJetResponseMaker* AddTaskJetRespPtHard(
   Double_t    jetBiasTrack       = 5,
   Double_t    jetBiasClus        = 5,
   UInt_t      matching           = AliJetResponseMaker::kGeometrical,
-  Double_t    maxDistance        = 0.25,
+  Double_t    maxDistance1       = 0.25,
+  Double_t    maxDistance2       = 0.25,
   UInt_t      type               = AliAnalysisTaskEmcal::kTPC,
   Int_t       maxPtHardBin       = -999,
   Int_t       minPtHardBin       = -999,
@@ -31,7 +32,7 @@ AliJetResponseMaker* AddTaskJetRespPtHard(
   for (Int_t i = minPtBin; i <= maxPtBin; i++) {
     AddTaskJetResponseMaker(ntracks1, nclusters1, njets1, nrho1, ntracks2, nclusters2, njets2, nrho2,
                            jetradius, jetptcut, jetareacut, jetBiasTrack, jetBiasClus, 
-                           matching, maxDistance, type, i, taskname, biggerMatrix, jetTask + i - minPtBin);
+                           matching, maxDistance1, maxDistance2, type, i, taskname, biggerMatrix, jetTask + i - minPtBin);
   }
   
   return jetTask;
index 7ea7288..3458981 100644 (file)
@@ -15,7 +15,8 @@ AliJetResponseMaker* AddTaskJetResponseMaker(
   Double_t    jetBiasTrack       = 5,
   Double_t    jetBiasClus        = 5,
   UInt_t      matching           = AliJetResponseMaker::kGeometrical,
-  Double_t    maxDistance        = 0.25,
+  Double_t    maxDistance1       = 0.25,
+  Double_t    maxDistance2       = 0.25,
   UInt_t      type               = AliAnalysisTaskEmcal::kTPC,
   Int_t       ptHardBin          = -999,
   const char *taskname           = "AliJetResponseMaker",
@@ -76,7 +77,7 @@ AliJetResponseMaker* AddTaskJetResponseMaker(
   jetTask->SetPercAreaCut(jetareacut);
   jetTask->SetPtBiasJetTrack(jetBiasTrack);
   jetTask->SetPtBiasJetClus(jetBiasClus);
-  jetTask->SetMatching(matching, maxDistance);
+  jetTask->SetMatching(matching, maxDistance1, maxDistance2);
   jetTask->SetVzRange(-10,10);
   jetTask->SetPtHardBin(ptHardBin);