PWGGA/EMCALJetTasks/AliHadCorrTask.cxx

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