PWGGA/EMCALJetTasks/AliJetEmbeddingTask.cxx

   1 // $Id$
   2 //
   3 // Jet embedding task.
   4 //
   5 // Author: Salvatore Aiola, Constantin Loizides
   6
   7 #include <TClonesArray.h>
   8 #include <TLorentzVector.h>
   9 #include <TRandom3.h>
  10
  11 #include "AliAnalysisManager.h"
  12 #include "AliVEvent.h"
  13 #include "AliVCluster.h"
  14 #include "AliEMCALDigit.h"
  15 #include "AliEMCALRecPoint.h"
  16 #include "AliPicoTrack.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliLog.h"
  19
  20 #include "AliJetEmbeddingTask.h"
  21
  22 ClassImp(AliJetEmbeddingTask)
  23
  24 //________________________________________________________________________
  25 AliJetEmbeddingTask::AliJetEmbeddingTask() :
  26   AliAnalysisTaskSE("AliJetEmbeddingTask"),
  27   fGeomName(),
  28   fTracksName(),
  29   fOutTracksName(),
  30   fCaloName(),
  31   fOutCaloName(),
  32   fEtaMin(-1),
  33   fEtaMax(1),
  34   fPhiMin(0),
  35   fPhiMax(TMath::Pi() * 2),
  36   fPtMin(0),
  37   fPtMax(0),
  38   fCopyArray(kTRUE),
  39   fNEmbClusters(0),
  40   fNEmbTracks(0),
  41   fGeom(0),
  42   fClusters(0),
  43   fOutClusters(0),
  44   fTracks(0),
  45   fOutTracks(0)
  46 {
  47   // Default constructor.
  48 }
  49
  50 //________________________________________________________________________
  51 AliJetEmbeddingTask::AliJetEmbeddingTask(const char *name) :
  52   AliAnalysisTaskSE(name),
  53   fGeomName(""),
  54   fTracksName("PicoTracks"),
  55   fOutTracksName("PicoTracksEmbedded"),
  56   fCaloName("CaloClustersCorr"),
  57   fOutCaloName("CaloClustersCorrEmbedded"),
  58   fEtaMin(-1),
  59   fEtaMax(1),
  60   fPhiMin(0),
  61   fPhiMax(TMath::Pi() * 2),
  62   fPtMin(50),
  63   fPtMax(60),
  64   fCopyArray(kTRUE),
  65   fNEmbClusters(0),
  66   fNEmbTracks(1),
  67   fGeom(0),
  68   fClusters(0),
  69   fOutClusters(0),
  70   fTracks(0),
  71   fOutTracks(0)
  72 {
  73   // Standard constructor.
  74
  75 }
  76
  77 //________________________________________________________________________
  78 AliJetEmbeddingTask::~AliJetEmbeddingTask()
  79 {
  80   // Destructor
  81 }
  82
  83 //________________________________________________________________________
  84 void AliJetEmbeddingTask::Init()
  85 {
  86   // Init task.
  87
  88   if (fNEmbTracks > 0) {
  89     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
  90
  91     if (strcmp(fTracks->GetClass()->GetName(), "AliPicoTrack")) {
  92       AliError("Can only embed PicoTracks!");
  93       return;
  94     }
  95
  96     if (!fTracks) {
  97       AliError(Form("Couldn't retrieve tracks with name %s!", fTracksName.Data()));
  98       return;
  99     }
 100
 101     if (!fOutTracks) {
 102       fOutTracksName = fTracksName;
 103       if (fCopyArray) {
 104         fOutTracksName += "Embedded";
 105         fOutTracks = new TClonesArray(*fTracks);
 106         fOutTracks->SetName(fOutTracksName);
 107       }
 108       else {
 109         fOutTracks = fTracks;
 110       }
 111     }
 112
 113     if (fCopyArray) {
 114       if (!(InputEvent()->FindListObject(fOutTracksName)))
 115         InputEvent()->AddObject(fOutTracks);
 116       fOutTracks->Clear();
 117     }
 118   }
 119
 120   if (fNEmbClusters > 0) {
 121     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 122
 123     if (!fClusters) {
 124       AliError(Form("Couldn't retrieve clusters with name %s!", fCaloName.Data()));
 125       return;
 126     }
 127
 128     if (!fOutClusters) {
 129       fOutCaloName = fCaloName;
 130       if (fCopyArray) {
 131         fOutCaloName += "Embedded";
 132         fOutClusters = new TClonesArray(*fClusters);
 133         fOutClusters->SetName(fOutCaloName);
 134       }
 135       else {
 136         fOutClusters = fClusters;
 137       }
 138     }
 139
 140     if (fCopyArray) {
 141       if (!(InputEvent()->FindListObject(fOutCaloName)))
 142         InputEvent()->AddObject(fOutClusters);
 143       fOutClusters->Clear();
 144     }
 145   }
 146
 147   if (!fGeom) {
 148     if (fGeomName.Length()>0) {
 149       fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 150       if (!fGeom)
 151         AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
 152     } else {
 153       fGeom = AliEMCALGeometry::GetInstance();
 154       if (!fGeom)
 155         AliError("Could not get default geometry!");
 156     }
 157   }
 158 }
 159
 160 //________________________________________________________________________
 161 void AliJetEmbeddingTask::Embed()
 162 {
 163   // Embed particles.
 164
 165   if (fNEmbClusters > 0 && fOutClusters) {
 166     const Int_t nClusters = fOutClusters->GetEntriesFast();
 167     TClonesArray digits("AliEMCALDigit", 1);
 168     for (Int_t i = 0; i < fNEmbClusters; ++i) {
 169       Double_t pt  = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 170       Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
 171       Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
 172
 173       TLorentzVector nPart;
 174       nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 175       Double_t e = nPart.E();
 176
 177       Int_t absId = 0;
 178       fGeom->GetAbsCellIdFromEtaPhi(eta, phi, absId);
 179
 180       if (absId == -1) {
 181         AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
 182                         " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 183                         eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 184         continue;
 185       }
 186
 187       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
 188       digit->SetId(absId);
 189       digit->SetIndexInList(0);
 190       digit->SetType(AliEMCALDigit::kHG);
 191       digit->SetAmplitude(e);
 192
 193       AliEMCALRecPoint recPoint("");
 194       recPoint.AddDigit(*digit, e, kFALSE);
 195       recPoint.EvalGlobalPosition(0, &digits);
 196
 197       TVector3 gpos;
 198       recPoint.GetGlobalPosition(gpos);
 199       Float_t g[3];
 200       gpos.GetXYZ(g);
 201
 202       AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters + i));
 203       cluster->SetType(AliVCluster::kEMCALClusterv1);
 204       cluster->SetE(recPoint.GetEnergy());
 205       cluster->SetPosition(g);
 206       cluster->SetNCells(1);
 207       UShort_t shortAbsId = absId;
 208       cluster->SetCellsAbsId(&shortAbsId);
 209       Double32_t fract = 1;
 210       cluster->SetCellsAmplitudeFraction(&fract);
 211       cluster->SetID(nClusters + i);
 212       cluster->SetEmcCpvDistance(-1);
 213       cluster->SetChi2(100); // MC flag!
 214     }
 215   }
 216
 217   if (fNEmbTracks > 0 && fOutTracks) {
 218     Int_t nTracks = fOutTracks->GetEntriesFast();
 219     for (Int_t i = 0; i < fNEmbTracks; ++i) {
 220       Double_t pt  = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 221       Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
 222       Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
 223
 224       new ((*fOutTracks)[nTracks + i]) AliPicoTrack(pt,
 225                                                     eta,
 226                                                     phi,
 227                                                     1,
 228                                                     100,    // MC flag!
 229                                                     0,
 230                                                     0,
 231                                                     kFALSE);
 232     }
 233   }
 234 }
 235
 236 //________________________________________________________________________
 237 void AliJetEmbeddingTask::UserExec(Option_t *)
 238 {
 239   // Execute per event.
 240
 241   Init();
 242
 243   Embed();
 244 }