PWG/EMCAL/AliEmcalClusTrackMatcherTask.cxx

   1 // $Id$
   2 //
   3 // Track/cluster matcher
   4 //
   5 // Author: C.Loizides, S.Aiola
   6
   7 #include "AliEmcalClusTrackMatcherTask.h"
   8
   9 #include <TClonesArray.h>
  10
  11 #include "AliAODCaloCluster.h"
  12 #include "AliESDCaloCluster.h"
  13 #include "AliEmcalParticle.h"
  14 #include "AliLog.h"
  15 #include "AliParticleContainer.h"
  16 #include "AliPicoTrack.h"
  17 #include "AliVCluster.h"
  18 #include "AliVTrack.h"
  19
  20 ClassImp(AliEmcalClusTrackMatcherTask)
  21
  22 //________________________________________________________________________
  23 AliEmcalClusTrackMatcherTask::AliEmcalClusTrackMatcherTask() :
  24   AliAnalysisTaskEmcal("AliEmcalClusTrackMatcherTask",kFALSE),
  25   fMaxDistance(0.1),
  26   fModifyObjs(kFALSE),
  27   fOrigTracks(0),
  28   fOrigClus(0),
  29   fHistMatchEtaAll(0),
  30   fHistMatchPhiAll(0)
  31 {
  32   // Constructor.
  33
  34   for(Int_t icent=0; icent<8; ++icent) {
  35     for(Int_t ipt=0; ipt<9; ++ipt) {
  36       for(Int_t ieta=0; ieta<2; ++ieta) {
  37         fHistMatchEta[icent][ipt][ieta] = 0;
  38         fHistMatchPhi[icent][ipt][ieta] = 0;
  39       }
  40     }
  41   }
  42 }
  43
  44 //________________________________________________________________________
  45 AliEmcalClusTrackMatcherTask::AliEmcalClusTrackMatcherTask(const char *name, Bool_t histo) :
  46   AliAnalysisTaskEmcal(name,histo),
  47   fMaxDistance(0.1),
  48   fModifyObjs(kFALSE),
  49   fOrigTracks(0),
  50   fOrigClus(0),
  51   fHistMatchEtaAll(0),
  52   fHistMatchPhiAll(0)
  53 {
  54   // Standard constructor.
  55
  56   for(Int_t icent=0; icent<8; ++icent) {
  57     for(Int_t ipt=0; ipt<9; ++ipt) {
  58       for(Int_t ieta=0; ieta<2; ++ieta) {
  59         fHistMatchEta[icent][ipt][ieta] = 0;
  60         fHistMatchPhi[icent][ipt][ieta] = 0;
  61       }
  62     }
  63   }
  64 }
  65
  66 //________________________________________________________________________
  67 AliEmcalClusTrackMatcherTask::~AliEmcalClusTrackMatcherTask()
  68 {
  69   // Destructor.
  70 }
  71
  72 //________________________________________________________________________
  73 void AliEmcalClusTrackMatcherTask::ExecOnce()
  74 {
  75   // Initialize the analysis.
  76
  77   if (fParticleCollArray.GetEntriesFast()<2) {
  78     AliError(Form("Wrong number of particle collections (%d), required 2",fParticleCollArray.GetEntriesFast()));
  79     return;
  80   }
  81
  82   for (Int_t i = 0; i < 2; i++) {
  83     AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
  84     cont->SetClassName("AliEmcalParticle");
  85     // make sure objects are not double matched
  86     TClonesArray *dummy = new TClonesArray("TObject",0);
  87     dummy->SetName(Form("%s_matched", cont->GetArrayName().Data()));
  88     AddObjectToEvent(dummy);
  89     // get pointer to original collections
  90     TString tmp(cont->GetArrayName());
  91     TObjArray *arr = tmp.Tokenize("_");
  92     if (arr) {
  93       const Int_t aid = arr->GetEntries()-1;
  94       if (aid>0) {
  95         TString tname(arr->At(aid)->GetName());
  96         TClonesArray *oarr =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(tname));
  97         if (oarr && (i==0)) {
  98           AliInfo(Form("Setting orig tracks to %s", tname.Data()));
  99           fOrigTracks = oarr;
 100         } else if (oarr && (i==1)) {
 101           AliInfo(Form("Setting orig clusters to %s", tname.Data()));
 102           fOrigClus = oarr;
 103         }
 104       }
 105     }
 106     delete arr;
 107   }
 108
 109   AliAnalysisTaskEmcal::ExecOnce();
 110 }
 111
 112 //________________________________________________________________________
 113 void AliEmcalClusTrackMatcherTask::UserCreateOutputObjects()
 114 {
 115   // Create my user objects.
 116
 117   if (!fCreateHisto)
 118     return;
 119
 120   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 121
 122   const Int_t nCentChBins = fNcentBins * 2;
 123
 124   fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
 125   fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
 126   fOutput->Add(fHistMatchEtaAll);
 127   fOutput->Add(fHistMatchPhiAll);
 128
 129   for(Int_t icent=0; icent<nCentChBins; ++icent) {
 130     for(Int_t ipt=0; ipt<9; ++ipt) {
 131       for(Int_t ieta=0; ieta<2; ++ieta) {
 132         TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
 133         fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
 134         fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
 135         TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
 136         fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
 137         fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
 138         fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
 139         fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
 140       }
 141     }
 142   }
 143
 144   PostData(1, fOutput);
 145 }
 146
 147 //________________________________________________________________________
 148 Int_t AliEmcalClusTrackMatcherTask::GetMomBin(Double_t p) const
 149 {
 150   // Get momenum bin.
 151
 152   Int_t pbin=-1;
 153   if (p<0.5)
 154     pbin=0;
 155   else if (p>=0.5 && p<1.0)
 156     pbin=1;
 157   else if (p>=1.0 && p<1.5)
 158     pbin=2;
 159   else if (p>=1.5 && p<2.)
 160     pbin=3;
 161   else if (p>=2. && p<3.)
 162     pbin=4;
 163   else if (p>=3. && p<4.)
 164     pbin=5;
 165   else if (p>=4. && p<5.)
 166     pbin=6;
 167   else if (p>=5. && p<8.)
 168     pbin=7;
 169   else if (p>=8.)
 170     pbin=8;
 171
 172   return pbin;
 173 }
 174
 175 //________________________________________________________________________
 176 Bool_t AliEmcalClusTrackMatcherTask::Run()
 177 {
 178   // Run the matching for the selected options.
 179
 180   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 181   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 182
 183   AliEmcalParticle *partC = 0;
 184   AliEmcalParticle *partT = 0;
 185
 186   const Double_t maxd2 = fMaxDistance*fMaxDistance;
 187
 188   // set the links between tracks and clusters
 189   clusters->ResetCurrentID();
 190   while ((partC = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle()))) {
 191     AliVCluster *clust = partC->GetCluster();
 192
 193     tracks->ResetCurrentID();
 194     while ((partT = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle()))) {
 195       AliVTrack *track = partT->GetTrack();
 196       Double_t deta = 999;
 197       Double_t dphi = 999;
 198       AliPicoTrack::GetEtaPhiDiff(track, clust, dphi, deta);
 199       Double_t d2 = deta * deta + dphi * dphi;
 200       if (d2 > maxd2)
 201         continue;
 202
 203       Double_t d = TMath::Sqrt(d2);
 204       partC->AddMatchedObj(tracks->GetCurrentID(), d);
 205       partC->SetMatchedPtr(fOrigTracks);
 206       partT->AddMatchedObj(clusters->GetCurrentID(), d);
 207       partT->SetMatchedPtr(fOrigClus);
 208
 209       if (fCreateHisto) {
 210         Int_t mombin = GetMomBin(track->P());
 211         Int_t centbinch = fCentBin;
 212         if (track->Charge()<0)
 213           centbinch += fNcentBins;
 214         Int_t etabin = 0;
 215         if(track->Eta() > 0)
 216           etabin = 1;
 217
 218         fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
 219         fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
 220         fHistMatchEtaAll->Fill(deta);
 221         fHistMatchPhiAll->Fill(dphi);
 222       }
 223     }
 224   }
 225
 226   if (!fModifyObjs)
 227     return kTRUE;
 228
 229   clusters->ResetCurrentID();
 230   while ((partC = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle()))) {
 231     AliVCluster *clust = partC->GetCluster();
 232     clust->SetEmcCpvDistance(-1);
 233     clust->SetTrackDistance(1024, 1024);
 234     AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(clust);
 235     AliESDCaloCluster *ec = 0;
 236     if (ac) {
 237       const Int_t N = ac->GetNTracksMatched();
 238       for (Int_t i=N-1; i>=0; --i) {
 239         TObject *ptr = ac->GetTrackMatched(i);
 240         ac->RemoveTrackMatched(ptr);
 241       }
 242     } else {
 243       ec = dynamic_cast<AliESDCaloCluster*>(clust);
 244       TArrayI *arr = ec->GetTracksMatched();
 245       arr->Set(0);
 246     }
 247     const Int_t N = partC->GetNumberOfMatchedObj();
 248     if (N <= 0)
 249       continue;
 250     const UInt_t matchedId = partC->GetMatchedObjId();
 251     partT = static_cast<AliEmcalParticle*>(tracks->GetParticle(matchedId));
 252     AliVTrack   *track = partT->GetTrack();
 253     Double_t deta = 999;
 254     Double_t dphi = 999;
 255     AliPicoTrack::GetEtaPhiDiff(track, clust, dphi, deta);
 256     clust->SetEmcCpvDistance(matchedId);
 257     clust->SetTrackDistance(deta, dphi);
 258     if (ac) {
 259       for (Int_t i=0; i<N; ++i) {
 260         Int_t id = partC->GetMatchedObjId(i);
 261         partT = static_cast<AliEmcalParticle*>(tracks->GetParticle(id));
 262         TObject *obj = partT->GetTrack();
 263         ac->AddTrackMatched(obj);
 264       }
 265     } else {
 266       TArrayI arr(N);
 267       for (Int_t i=0; i<N; ++i) {
 268         Int_t id = partC->GetMatchedObjId(i);
 269         partT = static_cast<AliEmcalParticle*>(tracks->GetParticle(id));
 270         arr.AddAt(partT->IdInCollection(),i);
 271       }
 272       ec->AddTracksMatched(arr);
 273     }
 274   }
 275
 276   tracks->ResetCurrentID();
 277   while ((partT = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle()))) {
 278     AliVTrack *track = partT->GetTrack();
 279     track->ResetStatus(AliVTrack::kEMCALmatch);
 280     if (partT->GetNumberOfMatchedObj() <= 0)
 281       continue;
 282     partC = static_cast<AliEmcalParticle*>(clusters->GetParticle(partT->GetMatchedObjId()));
 283     track->SetEMCALcluster(partC->IdInCollection());
 284     track->SetStatus(AliVTrack::kEMCALmatch);
 285   }
 286
 287   return kTRUE;
 288 }