PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx

   1 // $Id$
   2 //
   3 // Jet modelling task.
   4 //
   5 // Author: S.Aiola, C.Loizides
   6
   7 #include "AliJetModelBaseTask.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TH1.h>
  11 #include <TLorentzVector.h>
  12 #include <TRandom3.h>
  13
  14 #include "AliAODCaloCluster.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliEMCALDigit.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliEMCALRecPoint.h"
  19 #include "AliESDCaloCluster.h"
  20 #include "AliLog.h"
  21 #include "AliPicoTrack.h"
  22 #include "AliVCluster.h"
  23 #include "AliVEvent.h"
  24
  25 ClassImp(AliJetModelBaseTask)
  26
  27 //________________________________________________________________________
  28 AliJetModelBaseTask::AliJetModelBaseTask() :
  29   AliAnalysisTaskSE("AliJetModelBaseTask"),
  30   fGeomName(),
  31   fTracksName(),
  32   fOutTracksName(),
  33   fCaloName(),
  34   fOutCaloName(),
  35   fSuffix(),
  36   fEtaMin(-1),
  37   fEtaMax(1),
  38   fPhiMin(0),
  39   fPhiMax(TMath::Pi() * 2),
  40   fPtMin(0),
  41   fPtMax(0),
  42   fCopyArray(kTRUE),
  43   fNClusters(0),
  44   fNTracks(0),
  45   fMarkMC(kTRUE),
  46   fPtSpectrum(0),
  47   fIsInit(0),
  48   fGeom(0),
  49   fClusters(0),
  50   fOutClusters(0),
  51   fTracks(0),
  52   fOutTracks(0)
  53 {
  54   // Default constructor.
  55 }
  56
  57 //________________________________________________________________________
  58 AliJetModelBaseTask::AliJetModelBaseTask(const char *name) :
  59   AliAnalysisTaskSE(name),
  60   fGeomName(""),
  61   fTracksName("PicoTracks"),
  62   fOutTracksName("PicoTracksEmbedded"),
  63   fCaloName("CaloClustersCorr"),
  64   fOutCaloName("CaloClustersCorrEmbedded"),
  65   fSuffix("Processed"),
  66   fEtaMin(-1),
  67   fEtaMax(1),
  68   fPhiMin(0),
  69   fPhiMax(TMath::Pi() * 2),
  70   fPtMin(50),
  71   fPtMax(60),
  72   fCopyArray(kTRUE),
  73   fNClusters(0),
  74   fNTracks(1),
  75   fMarkMC(kTRUE),
  76   fPtSpectrum(0),
  77   fIsInit(0),
  78   fGeom(0),
  79   fClusters(0),
  80   fOutClusters(0),
  81   fTracks(0),
  82   fOutTracks(0)
  83 {
  84   // Standard constructor.
  85 }
  86
  87 //________________________________________________________________________
  88 AliJetModelBaseTask::~AliJetModelBaseTask()
  89 {
  90   // Destructor
  91 }
  92
  93 //________________________________________________________________________
  94 void AliJetModelBaseTask::UserExec(Option_t *)
  95 {
  96   // Execute per event.
  97
  98   if (!fIsInit)
  99     fIsInit = ExecOnce();
 100
 101   if (!fIsInit)
 102     return;
 103
 104   if (fCopyArray) {
 105     if (fOutTracks)
 106       fOutTracks->Delete();
 107     if (fOutClusters)
 108       fOutClusters->Delete();
 109   }
 110
 111   Run();
 112 }
 113
 114 //________________________________________________________________________
 115 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
 116 {
 117   // Add a cluster to the event.
 118
 119   Int_t absId = 0;
 120   if (eta < -100 || phi < 0) {
 121     GetRandomCell(eta, phi, absId);
 122   }
 123   else {
 124     fGeom->EtaPhiFromIndex(absId, eta, phi);
 125   }
 126
 127   if (absId == -1) {
 128     AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
 129                     " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 130                     eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 131     return 0;
 132   }
 133
 134   if (e < 0) {
 135     Double_t pt = GetRandomPt();
 136     TLorentzVector nPart;
 137     nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 138     e = nPart.E();
 139   }
 140
 141   return AddCluster(e, absId);
 142 }
 143
 144 //________________________________________________________________________
 145 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
 146 {
 147   // Add a cluster to the event.
 148
 149   const Int_t nClusters = fOutClusters->GetEntriesFast();
 150
 151   TClonesArray digits("AliEMCALDigit", 1);
 152
 153   AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
 154   digit->SetId(absId);
 155   digit->SetIndexInList(0);
 156   digit->SetType(AliEMCALDigit::kHG);
 157   digit->SetAmplitude(e);
 158
 159   AliEMCALRecPoint recPoint("");
 160   recPoint.AddDigit(*digit, e, kFALSE);
 161   recPoint.EvalGlobalPosition(0, &digits);
 162
 163   TVector3 gpos;
 164   recPoint.GetGlobalPosition(gpos);
 165   Float_t g[3];
 166   gpos.GetXYZ(g);
 167
 168   AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
 169   cluster->SetType(AliVCluster::kEMCALClusterv1);
 170   cluster->SetE(recPoint.GetEnergy());
 171   cluster->SetPosition(g);
 172   cluster->SetNCells(1);
 173   UShort_t shortAbsId = absId;
 174   cluster->SetCellsAbsId(&shortAbsId);
 175   Double32_t fract = 1;
 176   cluster->SetCellsAmplitudeFraction(&fract);
 177   cluster->SetID(nClusters);
 178   cluster->SetEmcCpvDistance(-1);
 179   if (fMarkMC)
 180     cluster->SetChi2(100); // MC flag!
 181
 182   return cluster;
 183 }
 184
 185 //________________________________________________________________________
 186 AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi)
 187 {
 188   // Add a track to the event.
 189
 190   const Int_t nTracks = fOutTracks->GetEntriesFast();
 191
 192   if (pt < 0)
 193     pt = GetRandomPt();
 194   if (eta < -100)
 195     eta = GetRandomEta();
 196   if (phi < 0)
 197     phi = GetRandomPhi();
 198
 199   Int_t label = fMarkMC ? 100 : 0;
 200
 201   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
 202                                                                   eta,
 203                                                                   phi,
 204                                                                   1,
 205                                                                   label,    // MC flag!
 206                                                                   0,
 207                                                                   0,
 208                                                                   kFALSE);
 209   return track;
 210 }
 211
 212 //________________________________________________________________________
 213 void AliJetModelBaseTask::CopyClusters()
 214 {
 215   // Copy all the clusters in the new collection
 216
 217   Bool_t esdMode = (Bool_t)(fClusters->GetClass()->GetBaseClass("AliESDCaloCluster") != 0);
 218   const Int_t nClusters = fClusters->GetEntriesFast();
 219   Int_t nCopiedClusters = 0;
 220
 221   if (esdMode) {
 222     for (Int_t i = 0; i < nClusters; ++i) {
 223       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
 224       if (!esdcluster || !esdcluster->IsEMCAL())
 225         continue;
 226       new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
 227       nCopiedClusters++;
 228     }
 229   }
 230   else {
 231     for (Int_t i = 0; i < nClusters; ++i) {
 232       AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
 233       if (!aodcluster || !aodcluster->IsEMCAL())
 234         continue;
 235       new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
 236       nCopiedClusters++;
 237     }
 238   }
 239 }
 240
 241 //________________________________________________________________________
 242 void AliJetModelBaseTask::CopyTracks()
 243 {
 244   const Int_t nTracks = fTracks->GetEntriesFast();
 245   Int_t nCopiedTracks = 0;
 246   for (Int_t i = 0; i < nTracks; ++i) {
 247     AliPicoTrack *track = static_cast<AliPicoTrack*>(fTracks->At(i));
 248     if (!track)
 249       continue;
 250     new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*track);
 251     nCopiedTracks++;
 252   }
 253 }
 254
 255 //________________________________________________________________________
 256 Bool_t AliJetModelBaseTask::ExecOnce()
 257 {
 258   // Init task.
 259
 260   if (fPtMax < fPtMin) {
 261     AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
 262     fPtMax = fPtMin;
 263   }
 264
 265   if (fEtaMax < fEtaMin) {
 266     AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
 267     fEtaMax = fEtaMin;
 268   }
 269
 270   if (fPhiMax < fPhiMin) {
 271     AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
 272     fPhiMax = fPhiMin;
 273   }
 274
 275   if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) {
 276     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 277     if (!fTracks) {
 278       AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
 279       return kFALSE;
 280     }
 281
 282     if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
 283       AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data()));
 284       fTracks = 0;
 285       return kFALSE;
 286     }
 287
 288     if (!fOutTracks) {
 289       fOutTracksName = fTracksName;
 290       if (fCopyArray) {
 291         fOutTracksName += fSuffix;
 292         fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
 293         fOutTracks->SetName(fOutTracksName);
 294         if (InputEvent()->FindListObject(fOutTracksName)) {
 295           AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data()));
 296           return kFALSE;
 297         }
 298         else {
 299           InputEvent()->AddObject(fOutTracks);
 300         }
 301       }
 302       else {
 303         fOutTracks = fTracks;
 304       }
 305     }
 306   }
 307
 308   if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) {
 309     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 310
 311     if (!fClusters) {
 312       AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
 313       return kFALSE;
 314     }
 315
 316     if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
 317       AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data()));
 318       fClusters = 0;
 319       return kFALSE;
 320     }
 321
 322     if (!fOutClusters) {
 323       fOutCaloName = fCaloName;
 324       if (fCopyArray) {
 325         fOutCaloName += fSuffix;
 326         fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
 327         fOutClusters->SetName(fOutCaloName);
 328         if (InputEvent()->FindListObject(fOutCaloName)) {
 329           AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data()));
 330           return kFALSE;
 331         }
 332         else {
 333           InputEvent()->AddObject(fOutClusters);
 334         }
 335       }
 336       else {
 337         fOutClusters = fClusters;
 338       }
 339     }
 340
 341     if (!fGeom) {
 342       if (fGeomName.Length() > 0) {
 343         fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 344         if (!fGeom)
 345           AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
 346       } else {
 347         fGeom = AliEMCALGeometry::GetInstance();
 348         if (!fGeom)
 349           AliError("Could not get default geometry!");
 350       }
 351     }
 352
 353     const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 354     const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 355     const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 356     const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 357
 358     if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
 359     if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
 360     if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
 361     if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
 362
 363     if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
 364     if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
 365     if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
 366     if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
 367   }
 368
 369   return kTRUE;
 370 }
 371
 372 //________________________________________________________________________
 373 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 374 {
 375   // Get random cell.
 376
 377   Int_t repeats = 0;
 378   Double_t rndEta = eta;
 379   Double_t rndPhi = phi;
 380   do {
 381     if (eta < -100)
 382       rndEta = GetRandomEta();
 383     if (phi < 0)
 384       rndPhi = GetRandomPhi();
 385     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
 386     repeats++;
 387   } while (absId == -1 && repeats < 100);
 388
 389   if (!(absId > -1)) {
 390     AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
 391                     "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 392   }
 393   else {
 394     eta = rndEta;
 395     phi = rndPhi;
 396   }
 397 }
 398
 399 //________________________________________________________________________
 400 Double_t AliJetModelBaseTask::GetRandomEta()
 401 {
 402   // Get random eta.
 403
 404   return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
 405 }
 406
 407 //________________________________________________________________________
 408 Double_t AliJetModelBaseTask::GetRandomPhi()
 409 {
 410   // Get random phi.
 411
 412   return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
 413 }
 414
 415 //________________________________________________________________________
 416 Double_t AliJetModelBaseTask::GetRandomPt()
 417 {
 418   // Get random pt.
 419
 420   if (fPtSpectrum)
 421     return fPtSpectrum->GetRandom();
 422   else
 423     return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 424 }
 425
 426 //________________________________________________________________________
 427 void AliJetModelBaseTask::Run()
 428 {
 429   // Run.
 430 }