PWGGA/EMCALJetTasks/AliJetModelBaseTask.cxx

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