]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/EMCALJetTasks/AliHadCorrTask.cxx
including task for user mconnors
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliHadCorrTask.cxx
index 2866b4fb023064e7220ebadeed042e49bbfffa4c..5ace53e4aad1c75231810859e7c8129ca2ef697a 100644 (file)
@@ -39,9 +39,7 @@ AliHadCorrTask::AliHadCorrTask() :
   fHistEafter(0),
   fHistEoPCent(0),
   fHistNMatchCent(0),
-  fHistNMatchCent_trk(0),
-  fHistCentrality(0),
-  fHistNoMatchEtaPhi(0)
+  fHistCentrality(0)
 {
   // Default constructor.
 
@@ -64,56 +62,13 @@ AliHadCorrTask::AliHadCorrTask() :
   } 
 }
 
-//________________________________________________________________________
-AliHadCorrTask::AliHadCorrTask(const char *name) : 
-  AliAnalysisTaskEmcal(name, kFALSE),
-  fOutCaloName("CaloClustersCorr"),
-  fPhiMatch(0.05),
-  fEtaMatch(0.025),
-  fDoTrackClus(0),
-  fHadCorr(0),
-  fEexclCell(0),
-  fOutClusters(0),
-  fHistNclusvsCent(0),
-  fHistNclusMatchvsCent(0),
-  fHistEbefore(0),
-  fHistEafter(0),
-  fHistEoPCent(0),
-  fHistNMatchCent(0),
-  fHistNMatchCent_trk(0),
-  fHistCentrality(0),
-  fHistNoMatchEtaPhi(0)
-
-{
-  // Standard constructor.
-
-  for(Int_t i=0; i<8; i++) {
-      fHistEsubPch[i]    = 0;
-      fHistEsubPchRat[i] = 0;
-    for(Int_t j=0; j<3; j++) {
-    }
-    
-    if (i<4) {
-      fHistMatchEvsP[i]   = 0;
-      fHistMatchdRvsEP[i] = 0;
-    }
-    
-    for(Int_t j=0; j<9; j++) {
-      for(Int_t k=0; k<2; k++) {
-       fHistMatchEtaPhi[i][j][k] = 0;
-      }
-    } 
-  }
-  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
-}
-
 //________________________________________________________________________
 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) : 
   AliAnalysisTaskEmcal(name, histo),
   fOutCaloName("CaloClustersCorr"),
   fPhiMatch(0.05),
   fEtaMatch(0.025),
-  fDoTrackClus(0),
+  fDoTrackClus(1),
   fHadCorr(0),
   fEexclCell(0),
   fOutClusters(0),
@@ -123,10 +78,7 @@ AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
   fHistEafter(0),
   fHistEoPCent(0),
   fHistNMatchCent(0),
-  fHistNMatchCent_trk(0),
-  fHistCentrality(0),
-  fHistNoMatchEtaPhi(0)
-
+  fHistCentrality(0)
 {
   // Standard constructor.
 
@@ -396,7 +348,7 @@ void AliHadCorrTask::UserCreateOutputObjects()
 
   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
   if (!handler) {
-    AliError("Input handler not available!");
+    AliError(Form("%s: Input handler not available!", GetName()));
     return;
   }
  
@@ -405,7 +357,7 @@ void AliHadCorrTask::UserCreateOutputObjects()
   } else if (handler->InheritsFrom("AliAODInputHandler")) {
     fOutClusters = new TClonesArray("AliAODCaloCluster");
   } else {
-    AliError("Input handler not recognized!");
+    AliError(Form("%s: Input handler not recognized!", GetName()));
     return;
   }
   fOutClusters->SetName(fOutCaloName);
@@ -460,12 +412,10 @@ void AliHadCorrTask::UserCreateOutputObjects()
   fHistCentrality       = new TH1F("fHistCentrality",  "Centrality",       100, 0, 100);
   fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
   fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
-  fHistEbefore          = new TH2F("Ebefore",          "Ebefore",          100, 0, 100, 101, -0.5, 100.5);
-  fHistEafter           = new TH2F("Eafter",           "Eafter",           100, 0, 100, 101, -0.5, 100.5);
+  fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
+  fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
   fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
   fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 11, -0.5, 10.5);
-  fHistNMatchCent_trk   = new TH2F("NMatchesCent_trk", "NMatchesCent_trk", 100, 0, 100, 11, -0.5, 10.5);
-  fHistNoMatchEtaPhi    = new TH2F("NoMatchEtaPhi",    "NoMatchEtaPhi",    200, -1, 1, 90, 1, 4);
 
   fOutput->Add(fHistNclusMatchvsCent);
   fOutput->Add(fHistNclusvsCent);
@@ -473,61 +423,11 @@ void AliHadCorrTask::UserCreateOutputObjects()
   fOutput->Add(fHistEafter);
   fOutput->Add(fHistEoPCent);
   fOutput->Add(fHistNMatchCent);
-  fOutput->Add(fHistNMatchCent_trk);
   fOutput->Add(fHistCentrality);
-  fOutput->Add(fHistNoMatchEtaPhi);
 
   PostData(1, fOutput);
 }
 
-//________________________________________________________________________
-void AliHadCorrTask::DoTrackClusLoop() 
-{
-  // Do the track cluster loop.
-
-  const Int_t Ntrks = fTracks->GetEntries();
-
-  // loop over all tracks
-  for (Int_t t = 0; t < Ntrks; ++t) {
-    AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(t));
-    if (!emctrack)
-      continue;
-
-    AliVTrack *track = emctrack->GetTrack();
-    if (!track)
-      continue;
-    if (!AcceptTrack(track))
-      continue;
-
-    Int_t Nclus    = emctrack->GetNumberOfMatchedObj();
-    Int_t Nmatches = 0;
-
-    // loop over matched clusters
-    for (Int_t i = 0; i < Nclus; ++i) {
-      Int_t c = emctrack->GetMatchedObjId(i);
-      AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(c));
-      if (!emccluster)
-       continue;
-
-      AliVCluster *cluster = emccluster->GetCluster();
-      if (!cluster)
-       continue;
-
-      Double_t etadiff = 999;
-      Double_t phidiff = 999;
-      AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
-
-      if (TMath::Abs(phidiff) < 0.050 && TMath::Abs(etadiff) < 0.025) // pp cuts!!!
-       Nmatches++;
-    }
-
-    fHistNMatchCent_trk->Fill(fCent, Nmatches);
-
-    if(Nmatches == 0 && track->Pt() > 2.0)
-      fHistNoMatchEtaPhi->Fill(track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal());
-  }
-}
-
 //________________________________________________________________________
 void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches) 
 {
@@ -546,13 +446,16 @@ void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t
     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
     if (!emctrack)
       continue;
-
     AliVTrack *track = emctrack->GetTrack();
     if (!track)
       continue;
     if (!AcceptTrack(track))
       continue;
 
+    // check if track also points to cluster
+    if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
+      continue;
+
     Double_t etadiff = 999;
     Double_t phidiff = 999;
     AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
@@ -566,7 +469,6 @@ void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t
     if (fCreateHisto) {
       Int_t etabin = 0;
       if(track->Eta() > 0) etabin=1;
-       
       fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
     }
     
@@ -576,30 +478,26 @@ void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t
 
     if (fPhiMatch > 0) {
       phiCutlo = -fPhiMatch;
-      phiCuthi = fPhiMatch;
-    }
-    else {
+      phiCuthi = +fPhiMatch;
+    } else {
       phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
       phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
     }
 
     if (fEtaMatch > 0) {
       etaCut = fEtaMatch;
-    }
-    else {
+    } else {
       etaCut = GetEtaSigma(mombin);
     }
 
     if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
-      if ((fDoTrackClus && (track->GetEMCALcluster()) == iClus) || !fDoTrackClus) {
-       ++Nmatches;
-       totalTrkP += track->P();
-
-       if (fCreateHisto) {
-         if (fHadCorr > 1 && mombin > -1) {
-           fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
-         }
-       }
+      ++Nmatches;
+      totalTrkP += track->P();
+
+      if (fCreateHisto) {
+        if (fHadCorr > 1 && mombin > -1) {
+          fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
+        }
       }
     }
   }
@@ -628,18 +526,15 @@ Bool_t AliHadCorrTask::Run()
       am->LoadBranch("CaloClusters");
     if (fTracksName == "Tracks")
       am->LoadBranch("Tracks");
-    am->LoadBranch("Centrality");      
+    am->LoadBranch("Centrality.");      
   }
 
-  if (fCreateHisto)
+  if (fCreateHisto) {
     fHistCentrality->Fill(fCent);
+  }
 
+   // loop over all clusters
   const Int_t Nclus = fCaloClusters->GetEntries();
-  if (fDoTrackClus && fCreateHisto)
-    DoTrackClusLoop();
-
-  // loop over all clusters
   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
     AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
     if (!emccluster)
@@ -652,7 +547,8 @@ Bool_t AliHadCorrTask::Run()
       continue;
 
     Double_t energyclus = 0;
-    fHistEbefore->Fill(fCent, cluster->E());
+    if (fCreateHisto)
+      fHistEbefore->Fill(fCent, cluster->E());
   
     // apply correction / subtraction
     if (fHadCorr > 0) {
@@ -698,10 +594,8 @@ Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Doub
   // Apply the hadronic correction with one track only.
 
   AliVCluster *cluster = emccluster->GetCluster();
-  
-  Double_t energyclus = cluster->E();
-  Int_t    iMin       = emccluster->GetMatchedObjId();
-
+  Double_t energyclus  = cluster->E();
+  Int_t    iMin        = emccluster->GetMatchedObjId();
   if (iMin < 0)
     return energyclus;
 
@@ -709,8 +603,9 @@ Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Doub
   if (!emctrack)
     return energyclus;
 
-  Int_t cid      = emctrack->GetMatchedObjId();
-  if (cid!=iMin) 
+  // check if track also points to cluster
+  Int_t cid = emctrack->GetMatchedObjId();
+  if (fDoTrackClus && (cid!=emccluster->IdInCollection())) 
     return energyclus;
 
   AliVTrack *track = emctrack->GetTrack();
@@ -731,7 +626,7 @@ Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Doub
   if (track->Charge()<0) 
     centbinch += 4;
 
-  // Plot some histograms if switched on
+  // plot some histograms if switched on
   if (fCreateHisto) {
     Int_t etabin = 0;
     if(track->Eta() > 0) 
@@ -746,13 +641,13 @@ Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Doub
     }
   }
           
-  // Define eta/phi cuts
+  // define eta/phi cuts
   Double_t etaCut   = 0.0;
   Double_t phiCutlo = 0.0;
   Double_t phiCuthi = 0.0;
   if (fPhiMatch > 0) {
     phiCutlo = -fPhiMatch;
-    phiCuthi = fPhiMatch;
+    phiCuthi = +fPhiMatch;
   }
   else {
     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
@@ -765,12 +660,9 @@ Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Doub
     etaCut = GetEtaSigma(mombin);
   }
   
-  // Apply the correction if the track is in the eta/phi window
+  // apply the correction if the track is in the eta/phi window
   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
-
-    if ((fDoTrackClus && (track->GetEMCALcluster()) == emccluster->IdInCollection()) || !fDoTrackClus) {
-      energyclus -= hadCorr * mom;
-    }
+    energyclus -= hadCorr * mom;
   }
 
   return energyclus;
@@ -789,32 +681,47 @@ Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Dou
   Double_t totalTrkP  = 0.0; // count total track momentum
   Int_t    Nmatches   = 0;   // count total number of matches
   
-  // Do the loop over the matched tracks and get the number of matches and the total momentum
+  // do the loop over the matched tracks and get the number of matches and the total momentum
   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
 
+  if(fCreateHisto && Nmatches == 0) {
+    fHistNclusvsCent->Fill(fCent);
+    fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
+  }
+
   if (totalTrkP <= 0)
     return energyclus;
 
   Double_t Esub = hadCorr * totalTrkP;
 
+  Double_t clusEexcl = fEexclCell * cNcells;
+
+  if (energyclus < clusEexcl) 
+    clusEexcl = energyclus;
+       
+  if (Esub > energyclus) 
+    Esub = energyclus;
+       
+  // applying Peter's proposed algorithm
+  // never subtract the full energy of the cluster 
+  if ((energyclus - Esub) < clusEexcl) 
+    Esub = (energyclus - clusEexcl);
+
   Double_t EoP = energyclus / totalTrkP;
 
-  // Plot some histograms if switched on
-  if(fCreateHisto) {
-    fHistNclusvsCent->Fill(fCent);
+  // plot some histograms if switched on
+  if (fCreateHisto) {
     fHistNMatchCent->Fill(fCent, Nmatches);
     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
     
     if(Nmatches > 0) 
       fHistNclusMatchvsCent->Fill(fCent);
       
-    if(Nmatches == 0)
-      fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
-    else if(Nmatches == 1)
+    if(Nmatches == 1)
       fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);       
     else if(Nmatches == 2)
       fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
-    else
+    else if(Nmatches > 2)
       fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
 
     if (EoP > 0) {
@@ -834,27 +741,10 @@ Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Dou
       fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
     } 
   }
-  Double_t clusEexcl = fEexclCell * cNcells;
 
-  if (energyclus < clusEexcl) 
-    clusEexcl = energyclus;
-       
-  if (Esub > energyclus) 
-    Esub = energyclus;
-       
-  //applying Peter's proposed algorithm
-  //Never subtract the full energy of the cluster 
-  if ((energyclus - Esub) < clusEexcl) 
-    Esub = (energyclus - clusEexcl);
-
-  //apply the correction
+  // apply the correction
   energyclus -= Esub;
 
   return energyclus;
 }
 
-void AliHadCorrTask::Terminate(Option_t *)
-{
-  // Nothing to be done.
-}