PWGGA/EMCALJetTasks/AliJetModelBaseTask.cxx

   1 // $Id$
   2 //
   3 // Jet modelling task.
   4 //
   5 // Author: S.Aiola, C.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 "AliJetModelBaseTask.h"
  21
  22 ClassImp(AliJetModelBaseTask)
  23
  24 //________________________________________________________________________
  25 AliJetModelBaseTask::AliJetModelBaseTask() :
  26   AliAnalysisTaskSE("AliJetModelBaseTask"),
  27   fGeomName(),
  28   fTracksName(),
  29   fOutTracksName(),
  30   fCaloName(),
  31   fOutCaloName(),
  32   fSuffix(),
  33   fEtaMin(-1),
  34   fEtaMax(1),
  35   fPhiMin(0),
  36   fPhiMax(TMath::Pi() * 2),
  37   fPtMin(0),
  38   fPtMax(0),
  39   fCopyArray(kTRUE),
  40   fNClusters(0),
  41   fNTracks(0),
  42   fGeom(0),
  43   fClusters(0),
  44   fOutClusters(0),
  45   fTracks(0),
  46   fOutTracks(0)
  47 {
  48   // Default constructor.
  49 }
  50
  51 //________________________________________________________________________
  52 AliJetModelBaseTask::AliJetModelBaseTask(const char *name) :
  53   AliAnalysisTaskSE(name),
  54   fGeomName(""),
  55   fTracksName("PicoTracks"),
  56   fOutTracksName("PicoTracksEmbedded"),
  57   fCaloName("CaloClustersCorr"),
  58   fOutCaloName("CaloClustersCorrEmbedded"),
  59   fSuffix("Processed"),
  60   fEtaMin(-1),
  61   fEtaMax(1),
  62   fPhiMin(0),
  63   fPhiMax(TMath::Pi() * 2),
  64   fPtMin(50),
  65   fPtMax(60),
  66   fCopyArray(kTRUE),
  67   fNClusters(0),
  68   fNTracks(1),
  69   fGeom(0),
  70   fClusters(0),
  71   fOutClusters(0),
  72   fTracks(0),
  73   fOutTracks(0)
  74 {
  75   // Standard constructor.
  76 }
  77
  78 //________________________________________________________________________
  79 AliJetModelBaseTask::~AliJetModelBaseTask()
  80 {
  81   // Destructor
  82 }
  83
  84 //________________________________________________________________________
  85 void AliJetModelBaseTask::Init()
  86 {
  87   // Init task.
  88
  89   if (fPtMax < fPtMin) {
  90     AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
  91     fPtMax = fPtMin;
  92   }
  93
  94   if (fEtaMax < fEtaMin) {
  95     AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
  96     fEtaMax = fEtaMin;
  97   }
  98
  99   if (fPhiMax < fPhiMin) {
 100     AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
 101     fPhiMax = fPhiMin;
 102   }
 103
 104   if (fNTracks > 0) {
 105     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 106     if (!fTracks) {
 107       AliError(Form("Couldn't retrieve tracks with name %s!", fTracksName.Data()));
 108       return;
 109     }
 110
 111     if (strcmp(fTracks->GetClass()->GetName(), "AliPicoTrack")) {
 112       AliError("Can only embed PicoTracks!");
 113       return;
 114     }
 115
 116     if (!fOutTracks) {
 117       fOutTracksName = fTracksName;
 118       if (fCopyArray) {
 119         fOutTracksName += fSuffix;
 120         fOutTracks = new TClonesArray(*fTracks);
 121         fOutTracks->SetName(fOutTracksName);
 122       }
 123       else {
 124         fOutTracks = fTracks;
 125       }
 126     }
 127
 128     if (fCopyArray) {
 129       if (!(InputEvent()->FindListObject(fOutTracksName)))
 130         InputEvent()->AddObject(fOutTracks);
 131       fOutTracks->Clear();
 132     }
 133   }
 134
 135   if (fNClusters > 0) {
 136     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 137
 138     if (!fClusters) {
 139       AliError(Form("Couldn't retrieve clusters with name %s!", fCaloName.Data()));
 140       return;
 141     }
 142
 143     if (!fOutClusters) {
 144       fOutCaloName = fCaloName;
 145       if (fCopyArray) {
 146         fOutCaloName += fSuffix;
 147         fOutClusters = new TClonesArray(*fClusters);
 148         fOutClusters->SetName(fOutCaloName);
 149       }
 150       else {
 151         fOutClusters = fClusters;
 152       }
 153     }
 154
 155     if (fCopyArray) {
 156       if (!(InputEvent()->FindListObject(fOutCaloName)))
 157         InputEvent()->AddObject(fOutClusters);
 158       fOutClusters->Clear();
 159     }
 160
 161     if (!fGeom) {
 162       if (fGeomName.Length() > 0) {
 163         fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 164         if (!fGeom)
 165           AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
 166       } else {
 167         fGeom = AliEMCALGeometry::GetInstance();
 168         if (!fGeom)
 169           AliError("Could not get default geometry!");
 170       }
 171     }
 172
 173     const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 174     const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 175     const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 176     const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 177
 178     if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
 179     if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
 180     if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
 181     if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
 182
 183     if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
 184     if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
 185     if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
 186     if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
 187   }
 188 }
 189
 190 //________________________________________________________________________
 191 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 192 {
 193   // Get random cell.
 194
 195   Int_t repeats = 0;
 196   Double_t rndEta = eta;
 197   Double_t rndPhi = phi;
 198   do {
 199     if (eta < -100)
 200       rndEta = GetRandomEta();
 201     if (phi < 0)
 202       rndPhi = GetRandomPhi();
 203     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
 204     repeats++;
 205   } while (absId == -1 && repeats < 100);
 206
 207   if (!(absId > -1)) {
 208     AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
 209                     "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 210   }
 211   else {
 212     eta = rndEta;
 213     phi = rndPhi;
 214   }
 215 }
 216
 217 //________________________________________________________________________
 218 Double_t AliJetModelBaseTask::GetRandomEta()
 219 {
 220   // Get random eta.
 221
 222   return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
 223 }
 224
 225 //________________________________________________________________________
 226 Double_t AliJetModelBaseTask::GetRandomPhi()
 227 {
 228   // Get random phi.
 229
 230   return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
 231 }
 232
 233 //________________________________________________________________________
 234 Double_t AliJetModelBaseTask::GetRandomPt()
 235 {
 236   // Get random pt.
 237
 238   return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 239 }
 240
 241 //________________________________________________________________________
 242 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
 243 {
 244   // Add a cluster to the event.
 245
 246   Int_t absId = 0;
 247   if (eta < -100 || phi < 0) {
 248     GetRandomCell(eta, phi, absId);
 249   }
 250   else {
 251     fGeom->EtaPhiFromIndex(absId, eta, phi);
 252   }
 253
 254   if (absId == -1) {
 255     AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
 256                     " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 257                     eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 258     return 0;
 259   }
 260
 261   if (e < 0) {
 262     Double_t pt = GetRandomPt();
 263     TLorentzVector nPart;
 264     nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 265     e = nPart.E();
 266   }
 267
 268   return AddCluster(e, absId);
 269 }
 270
 271 //________________________________________________________________________
 272 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
 273 {
 274   // Add a cluster to the event.
 275
 276   const Int_t nClusters = fOutClusters->GetEntriesFast();
 277
 278   TClonesArray digits("AliEMCALDigit", 1);
 279
 280   AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
 281   digit->SetId(absId);
 282   digit->SetIndexInList(0);
 283   digit->SetType(AliEMCALDigit::kHG);
 284   digit->SetAmplitude(e);
 285
 286   AliEMCALRecPoint recPoint("");
 287   recPoint.AddDigit(*digit, e, kFALSE);
 288   recPoint.EvalGlobalPosition(0, &digits);
 289
 290   TVector3 gpos;
 291   recPoint.GetGlobalPosition(gpos);
 292   Float_t g[3];
 293   gpos.GetXYZ(g);
 294
 295   AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
 296   cluster->SetType(AliVCluster::kEMCALClusterv1);
 297   cluster->SetE(recPoint.GetEnergy());
 298   cluster->SetPosition(g);
 299   cluster->SetNCells(1);
 300   UShort_t shortAbsId = absId;
 301   cluster->SetCellsAbsId(&shortAbsId);
 302   Double32_t fract = 1;
 303   cluster->SetCellsAmplitudeFraction(&fract);
 304   cluster->SetID(nClusters);
 305   cluster->SetEmcCpvDistance(-1);
 306   cluster->SetChi2(100); // MC flag!
 307
 308   return cluster;
 309 }
 310
 311 //________________________________________________________________________
 312 AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi)
 313 {
 314   // Add a track to the event.
 315
 316   const Int_t nTracks = fOutTracks->GetEntriesFast();
 317
 318   if (pt < 0)
 319     pt = GetRandomPt();
 320   if (eta < -100)
 321     eta = GetRandomEta();
 322   if (phi < 0)
 323     phi = GetRandomPhi();
 324
 325   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
 326                                                 eta,
 327                                                 phi,
 328                                                 1,
 329                                                 100,    // MC flag!
 330                                                 0,
 331                                                 0,
 332                                                 kFALSE);
 333   return track;
 334 }
 335
 336 //________________________________________________________________________
 337 void AliJetModelBaseTask::Run()
 338 {
 339   // Run.
 340 }
 341
 342 //________________________________________________________________________
 343 void AliJetModelBaseTask::UserExec(Option_t *)
 344 {
 345   // Execute per event.
 346
 347   Init();
 348   Run();
 349 }