PWGGA/EMCALJetTasks/AliHadCorrTask.cxx

   1 // $Id$
   2 //
   3 // Hadronic correction task.
   4 //
   5 // Author: R.Reed, C.Loizides
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TParticle.h>
  14 #include "AliAODCaloCluster.h"
  15 #include "AliAODEvent.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliESDCaloCluster.h"
  19 #include "AliESDtrack.h"
  20 #include "AliFJWrapper.h"
  21 #include "AliPicoTrack.h"
  22 #include "AliVEventHandler.h"
  23 #include "AliHadCorrTask.h"
  24
  25 ClassImp(AliHadCorrTask)
  26
  27 //________________________________________________________________________
  28 AliHadCorrTask::AliHadCorrTask() :
  29   AliAnalysisTaskSE("AliHadCorrTask"),
  30   fTracksName(),
  31   fCaloName(),
  32   fOutCaloName(),
  33   fPhiMatch(0.05),
  34   fEtaMatch(0.025),
  35   fHadCorr(0),
  36   fMinPt(0.15),
  37   fOutClusters(0),
  38   fOutputList(0),
  39   fHistNclusvsCent(0),
  40   fHistNclusMatchvsCent(0),
  41   fHistEbefore(0),
  42   fHistEafter(0),
  43   fHistEoPCent(0),
  44   fHistNMatchCent(0)
  45 {
  46   // Default constructor.
  47
  48   for(Int_t i=0; i<8; i++) {
  49     if (i<4) {
  50       fHistMatchEvsP[i]   = 0;
  51       fHistMatchdRvsEP[i] = 0;
  52       fHistEsubPch[i]     = 0;
  53     }
  54     for(Int_t j=0; j<5; j++) {
  55       fHistMatchEtaPhi[i][j] = 0;
  56     }
  57   }
  58 }
  59
  60 //________________________________________________________________________
  61 AliHadCorrTask::AliHadCorrTask(const char *name) :
  62   AliAnalysisTaskSE(name),
  63   fTracksName("Tracks"),
  64   fCaloName("CaloClusters"),
  65   fOutCaloName("CaloClustersCorr"),
  66   fPhiMatch(0.05),
  67   fEtaMatch(0.025),
  68   fHadCorr(0),
  69   fMinPt(0.15),
  70   fOutClusters(0),
  71   fOutputList(0),
  72   fHistNclusvsCent(0),
  73   fHistNclusMatchvsCent(0),
  74   fHistEbefore(0),
  75   fHistEafter(0),
  76   fHistEoPCent(0),
  77   fHistNMatchCent(0)
  78 {
  79   // Standard constructor.
  80
  81    for(Int_t i=0; i<8; i++) {
  82     if (i<4) {
  83       fHistMatchEvsP[i]=0;
  84       fHistMatchdRvsEP[i]=0;
  85       fHistEsubPch[i]=0;
  86     }
  87     for(Int_t j=0; j<5; j++) {
  88       fHistMatchEtaPhi[i][j] = 0;
  89     }
  90   }
  91
  92   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  93
  94   DefineInput(0,TChain::Class());
  95   DefineOutput(1,TList::Class());
  96 }
  97
  98 //________________________________________________________________________
  99 AliHadCorrTask::~AliHadCorrTask()
 100 {
 101   // Destructor
 102 }
 103
 104 //________________________________________________________________________
 105 Int_t AliHadCorrTask::GetCentBin(Double_t cent) const
 106 {
 107   // Get centrality bin.
 108
 109   Int_t centbin = -1;
 110   if (cent>=0 && cent<10)
 111     centbin=0;
 112   else if (cent>=10 && cent<30)
 113     centbin=1;
 114   else if (cent>=30 && cent<50)
 115     centbin=2;
 116   else if (cent>=50 && cent<=100)
 117     centbin=3;
 118
 119   return centbin;
 120 }
 121
 122 //________________________________________________________________________
 123 Int_t AliHadCorrTask::GetMomBin(Double_t p) const
 124 {
 125   // Get momenum bin.
 126
 127   Int_t pbin=-1;
 128   if (p>=0 && p<0.5)
 129     pbin=0;
 130   else if (p>=0.5 && p<2.)
 131     pbin=1;
 132   else if (p>=2. && p<3.)
 133     pbin=2;
 134   else if (p>=3. && p<5.)
 135     pbin=3;
 136   else if (p>=5.)
 137     pbin=4;
 138
 139   return pbin;
 140 }
 141
 142 //________________________________________________________________________
 143 void AliHadCorrTask::UserCreateOutputObjects()
 144 {
 145   // Create my user objects.
 146
 147   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 148   if (!handler) {
 149     AliError("Input handler not available!");
 150     return;
 151   }
 152
 153   if (handler->InheritsFrom("AliESDInputHandler")) {
 154     fOutClusters = new TClonesArray("AliESDCaloCluster");
 155   } else if (handler->InheritsFrom("AliAODInputHandler")) {
 156     fOutClusters = new TClonesArray("AliAODCaloCluster");
 157   } else {
 158     AliError("Input handler not recognized!");
 159     return;
 160   }
 161   fOutClusters->SetName(fOutCaloName);
 162
 163   OpenFile(1);
 164   fOutputList = new TList();
 165   fOutputList->SetOwner();
 166
 167   for(Int_t icent=0; icent<8; ++icent) {
 168     for(Int_t ipt=0; ipt<5; ++ipt) {
 169       TString name(Form("fHistMatchEtaPhi_%i_%i",icent,ipt));
 170       fHistMatchEtaPhi[icent][ipt] = new TH2F(name, name, 400, -0.2, 0.2, 1600, -0.8, 0.8);
 171       fOutputList->Add(fHistMatchEtaPhi[icent][ipt]);
 172     }
 173
 174     if(icent<4){
 175       TString name(Form("fHistMatchEvsP_%i",icent));
 176       fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 177       fOutputList->Add(fHistMatchEvsP[icent]);
 178
 179       name = Form("fHistMatchdRvsEP_%i",icent);
 180       fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
 181       fOutputList->Add(fHistMatchdRvsEP[icent]);
 182
 183       name = Form("fHistEsubPch_%i",icent);
 184       fHistEsubPch[icent]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 185       fOutputList->Add(fHistEsubPch[icent]);
 186     }
 187   }
 188
 189   fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
 190   fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
 191   fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
 192   fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
 193   fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
 194   fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 101, -0.5, 100.5);
 195   fOutputList->Add(fHistNclusMatchvsCent);
 196   fOutputList->Add(fHistNclusvsCent);
 197   fOutputList->Add(fHistEbefore);
 198   fOutputList->Add(fHistEafter);
 199   fOutputList->Add(fHistEoPCent);
 200   fOutputList->Add(fHistNMatchCent);
 201
 202   PostData(1, fOutputList);
 203 }
 204
 205 //________________________________________________________________________
 206 void AliHadCorrTask::UserExec(Option_t *)
 207 {
 208   // Execute per event.
 209
 210   // post output in event if not yet present
 211   if (!(InputEvent()->FindListObject(fOutCaloName)))
 212     InputEvent()->AddObject(fOutClusters);
 213
 214   // delete output
 215   fOutClusters->Delete();
 216
 217   // esd or aod mode
 218   Bool_t esdMode = kTRUE;
 219   if (dynamic_cast<AliAODEvent*>(InputEvent()))
 220       esdMode = kFALSE;
 221
 222   // optimization in case autobranch loading is off
 223   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 224   if (fCaloName == "CaloClusters")
 225     am->LoadBranch("CaloClusters");
 226   if (fTracksName == "Tracks")
 227     am->LoadBranch("Tracks");
 228
 229   // get centrality
 230   Double_t cent = -1;
 231   AliCentrality *centrality = InputEvent()->GetCentrality() ;
 232   if (centrality)
 233     cent = centrality->GetCentralityPercentile("V0M");
 234   else
 235     cent=99; // probably pp data
 236   if (cent<0) {
 237     AliError(Form("Centrality negative: %f", cent));
 238     return;
 239   }
 240   Int_t centbin = GetCentBin(cent);
 241
 242   // get input collections
 243   TClonesArray *tracks = 0;
 244   TClonesArray *clus   = 0;
 245   TList *l = InputEvent()->GetList();
 246
 247   tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
 248   if (!tracks) {
 249     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 250     return;
 251   }
 252   clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
 253   if (!clus) {
 254     AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
 255     return;
 256   }
 257
 258   // get primary vertex
 259   Double_t vertex[3] = {0, 0, 0};
 260   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 261
 262   // loop over clusters and tracks
 263   const Int_t Nclus = clus->GetEntries();
 264   const Int_t Ntrks = tracks->GetEntries();
 265
 266   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
 267     AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
 268     if (!c)
 269       continue;
 270     if (!c->IsEMCAL())
 271       continue;
 272
 273     // make primary particle out of cluster
 274     TLorentzVector nPart;
 275     c->GetMomentum(nPart, vertex);
 276
 277     Double_t etclus = nPart.Pt();
 278     if (etclus<fMinPt)
 279       continue;
 280
 281     Double_t energyclus = nPart.P();
 282
 283     // reset cluster/track matching
 284     c->SetEmcCpvDistance(-1);
 285     c->SetTrackDistance(999,999);
 286
 287     // run cluster-track matching
 288     Int_t    imin       = -1;
 289     Double_t dEtaMin    = 1e9;
 290     Double_t dPhiMin    = 1e9;
 291     Double_t dRmin      = 1e9;
 292     Double_t totalTrkP  = 0.0; // count total track momentum
 293     Int_t    Nmatches   = 0;   // count total number of matches
 294     for (Int_t t = 0; t<Ntrks; ++t) {
 295       AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(t));
 296       if (!track)
 297         continue;
 298       if (!track->IsEMCAL())
 299         continue;
 300       if (track->Pt()<fMinPt)
 301         continue;
 302       Double_t etadiff=999;
 303       Double_t phidiff=999;
 304       AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
 305       Double_t dR = TMath::Sqrt(etadiff*etadiff+phidiff*phidiff);
 306       if(dR<dRmin){
 307         dEtaMin = etadiff;
 308         dPhiMin = phidiff;
 309         dRmin = dR;
 310         imin = t;
 311       }
 312
 313       if (fHadCorr>1) {
 314         Double_t mom = track->P();
 315         Int_t mombin = GetMomBin(mom);
 316         Int_t centbinch = centbin;
 317         if (track->Charge()==-1)
 318           centbinch += 4;
 319         if (mombin>-1) {
 320           fHistMatchEtaPhi[centbinch][mombin]->Fill(etadiff,phidiff);
 321           fHistMatchdRvsEP[centbin]->Fill(dR,energyclus/mom);
 322         }
 323       }
 324       if (TMath::Abs(phidiff)<fPhiMatch && TMath::Abs(etadiff)<fEtaMatch) {
 325         ++Nmatches;
 326         totalTrkP += track->P();
 327       }
 328     }
 329
 330     // store closest track
 331     c->SetEmcCpvDistance(imin);
 332     c->SetTrackDistance(dPhiMin, dEtaMin);
 333
 334     fHistNclusvsCent->Fill(cent);
 335     fHistEbefore->Fill(cent,energyclus);
 336     fHistNMatchCent->Fill(cent,Nmatches);
 337     if(Nmatches>0)
 338       fHistNclusMatchvsCent->Fill(cent);
 339
 340     // apply correction / subtraction
 341     if (fHadCorr>0) {
 342       // to subtract only the closest track set fHadCor to a %
 343       // to subtract all tracks within the cut set fHadCor to %+1
 344       if (fHadCorr>1) {
 345         if (totalTrkP>0) {
 346           Double_t EoP  = energyclus/totalTrkP;
 347           Double_t Esub = (fHadCorr-1)*totalTrkP;
 348           if (Esub>energyclus)
 349             Esub = energyclus;
 350           fHistEoPCent->Fill(cent,EoP);
 351           fHistMatchEvsP[centbin]->Fill(energyclus,EoP);
 352           fHistEsubPch[centbin]->Fill(totalTrkP,Esub/totalTrkP);
 353         }
 354         energyclus -= (fHadCorr-1)*totalTrkP;
 355       } else if (imin>=0) {
 356         AliVTrack *t = dynamic_cast<AliVTrack*>(tracks->At(imin));
 357         if (t) {
 358           Double_t mom = t->P();
 359           Int_t mombin = GetMomBin(mom);
 360           Int_t centbinch = centbin;
 361           if (t->Charge()==-1)
 362             centbinch += 4;
 363           fHistMatchEtaPhi[centbinch][mombin]->Fill(dEtaMin,dPhiMin);
 364           if (mom>0){
 365             fHistEoPCent->Fill(cent,energyclus/mom);
 366             fHistMatchEvsP[centbin]->Fill(energyclus,energyclus/mom);
 367             fHistMatchdRvsEP[centbin]->Fill(dRmin,energyclus/mom);
 368           }
 369           if (TMath::Abs(dPhiMin)<fPhiMatch && TMath::Abs(dEtaMin)<fEtaMatch) {
 370             energyclus -= fHadCorr*t->P();
 371           }
 372         }
 373       }
 374     }
 375     if (energyclus<0)
 376       energyclus = 0;
 377     fHistEafter->Fill(cent,energyclus);
 378
 379     if (energyclus>0) { // create corrected cluster
 380       AliVCluster *oc;
 381       if (esdMode) {
 382         oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*>(c)));
 383       } else {
 384         oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*>(c)));
 385       }
 386       oc->SetE(energyclus);
 387       ++clusCount;
 388     }
 389   }
 390 }
 391
 392 //________________________________________________________________________
 393 void AliHadCorrTask::Terminate(Option_t *)
 394 {
 395   // Nothing to be done.
 396 }
 397