PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx

   1 // $Id: AliJetEmbeddingFromAODTask.cxx $
   2 //
   3 // Jet embedding from AOD task.
   4 //
   5 // Author: S.Aiola, C.Loizides
   6
   7 #include "AliJetEmbeddingFromAODTask.h"
   8
   9 //#include <iostream>
  10
  11 #include <TFile.h>
  12 #include <TTree.h>
  13 #include <TClonesArray.h>
  14 #include <TObjArray.h>
  15 #include <TObjString.h>
  16 #include <TGrid.h>
  17 #include <TH2C.h>
  18
  19 #include "AliVEvent.h"
  20 #include "AliAODTrack.h"
  21 #include "AliESDtrack.h"
  22 #include "AliPicoTrack.h"
  23 #include "AliVTrack.h"
  24 #include "AliVCluster.h"
  25 #include "AliVCaloCells.h"
  26 #include "AliCentrality.h"
  27 #include "AliVHeader.h"
  28 #include "AliVVertex.h"
  29 #include "AliAODHeader.h"
  30 #include "AliLog.h"
  31 #include "AliInputEventHandler.h"
  32
  33 ClassImp(AliJetEmbeddingFromAODTask)
  34
  35 //________________________________________________________________________
  36 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
  37   AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
  38   fFileList(0),
  39   fAODTreeName(),
  40   fAODHeaderName(),
  41   fAODVertexName(),
  42   fAODTrackName(),
  43   fAODClusName(),
  44   fAODCellsName(),
  45   fMinCentrality(0),
  46   fMaxCentrality(10),
  47   fTriggerMask(AliVEvent::kAny),
  48   fZVertexCut(10),
  49   fIncludeNoITS(kTRUE),
  50   fTotalFiles(2200),
  51   fEsdTreeMode(kFALSE),
  52   fCurrentFileID(0),
  53   fCurrentAODFileID(-1),
  54   fCurrentAODFile(0),
  55   fAODHeader(0),
  56   fAODVertex(0),
  57   fAODTracks(0),
  58   fAODClusters(0),
  59   fAODCaloCells(0),
  60   fHistFileIDs(0)
  61 {
  62   // Default constructor.
  63   SetSuffix("AODEmbedding");
  64   fMarkMC = kFALSE;
  65   fAODfilterBits[0] = -1;
  66   fAODfilterBits[1] = -1;
  67 }
  68
  69 //________________________________________________________________________
  70 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
  71   AliJetModelBaseTask(name, drawqa),
  72   fFileList(0),
  73   fAODTreeName("aodTree"),
  74   fAODHeaderName("header"),
  75   fAODVertexName("vertices"),
  76   fAODTrackName("tracks"),
  77   fAODClusName("caloClusters"),
  78   fAODCellsName("emcalCells"),
  79   fMinCentrality(0),
  80   fMaxCentrality(10),
  81   fTriggerMask(AliVEvent::kAny),
  82   fZVertexCut(10),
  83   fIncludeNoITS(kTRUE),
  84   fTotalFiles(2200),
  85   fEsdTreeMode(kFALSE),
  86   fCurrentFileID(0),
  87   fCurrentAODFileID(-1),
  88   fCurrentAODFile(0),
  89   fAODHeader(0),
  90   fAODVertex(0),
  91   fAODTracks(0),
  92   fAODClusters(0),
  93   fAODCaloCells(0),
  94   fHistFileIDs(0)
  95 {
  96   // Standard constructor.
  97   SetSuffix("AODEmbedding");
  98   fMarkMC = kFALSE;
  99   fAODfilterBits[0] = -1;
 100   fAODfilterBits[1] = -1;
 101 }
 102
 103 //________________________________________________________________________
 104 AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
 105 {
 106   // Destructor
 107
 108   if (fCurrentAODFile) {
 109     fCurrentAODFile->Close();
 110     delete fCurrentAODFile;
 111   }
 112 }
 113
 114 //________________________________________________________________________
 115 void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
 116 {
 117   if (!fQAhistos)
 118     return;
 119
 120   AliJetModelBaseTask::UserCreateOutputObjects();
 121
 122   fHistFileIDs = new TH2C("fHistFileIDs", "fHistFileIDs", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
 123   fHistFileIDs->GetXaxis()->SetTitle("File no. (PYTHIA)");
 124   fHistFileIDs->GetYaxis()->SetTitle("File no. (Embedded AOD)");
 125   fOutput->Add(fHistFileIDs);
 126
 127   PostData(1, fOutput);
 128 }
 129
 130 //________________________________________________________________________
 131 Bool_t AliJetEmbeddingFromAODTask::UserNotify()
 132 {
 133   TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
 134   if (path.EndsWith("/"))
 135     path.Remove(path.Length()-1);
 136   path.Remove(path.Last('/'));
 137   path.Remove(0, path.Last('/')+1);
 138   if (!path.IsDec()) {
 139     AliWarning(Form("Could not extract file number from path %s", path.Data()));
 140     return kTRUE;
 141   }
 142   fCurrentFileID = path.Atoi();
 143   fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles;
 144   AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
 145   return kTRUE;
 146 }
 147
 148 //________________________________________________________________________
 149 Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
 150 {
 151   if (fAODTreeName.Contains("aod"))
 152     fEsdTreeMode = kFALSE;
 153   else
 154     fEsdTreeMode = kTRUE;
 155
 156   return AliJetModelBaseTask::ExecOnce();
 157 }
 158
 159 //________________________________________________________________________
 160 Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
 161 {
 162   if (fCurrentAODFile) {
 163     fCurrentAODFile->Close();
 164     delete fCurrentAODFile;
 165     fCurrentAODFile = 0;
 166   }
 167
 168   if (fCurrentAODFileID >= fFileList->GetEntriesFast())
 169     return kFALSE;
 170
 171   TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
 172   TString fileName = objFileName->GetString();
 173
 174   if (fileName.BeginsWith("alien://") && !gGrid) {
 175       AliInfo("Trying to connect to AliEn ...");
 176       TGrid::Connect("alien://");
 177   }
 178
 179   fCurrentAODFile = TFile::Open(fileName);
 180
 181   if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
 182     return kFALSE;
 183   }
 184
 185   if (fQAhistos)
 186     fHistFileIDs->Fill(fCurrentFileID, fCurrentAODFileID);
 187
 188   fCurrentAODFileID++;
 189   return kTRUE;
 190 }
 191
 192 //________________________________________________________________________
 193 Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
 194 {
 195   static TTree *tree = 0;
 196   static Int_t entry = 0;
 197
 198   Int_t attempts = 0;
 199
 200   do {
 201
 202     if (!fCurrentAODFile || !tree || entry >= tree->GetEntries()) {
 203       if (!OpenNextFile())
 204         return kFALSE;
 205
 206       tree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
 207       if (!tree)
 208         return kFALSE;
 209
 210       if (!fAODHeaderName.IsNull())
 211         tree->SetBranchAddress(fAODHeaderName, &fAODHeader);
 212
 213       if (!fAODVertexName.IsNull())
 214         tree->SetBranchAddress(fAODVertexName, &fAODVertex);
 215
 216       if (!fAODTrackName.IsNull())
 217         tree->SetBranchAddress(fAODTrackName, &fAODTracks);
 218
 219       if (!fAODClusName.IsNull())
 220         tree->SetBranchAddress(fAODClusName, &fAODClusters);
 221
 222       if (!fAODCellsName.IsNull())
 223         tree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
 224
 225       entry = 0;
 226     }
 227
 228     tree->GetEntry(entry);
 229     entry++;
 230     attempts++;
 231
 232     if (attempts == 1000)
 233       AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
 234
 235   } while (!IsAODEventSelected() && tree);
 236
 237   return (tree!=0);
 238 }
 239
 240 //________________________________________________________________________
 241 Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
 242 {
 243   if (!fEsdTreeMode && fAODHeader) {
 244     AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
 245
 246     UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
 247
 248     if ((offlineTrigger & fTriggerMask) == 0) {
 249       AliDebug(2, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
 250       return kFALSE;
 251     }
 252
 253     AliCentrality *cent = aodHeader->GetCentralityP();
 254     Float_t centVal = cent->GetCentralityPercentile("V0M");
 255     if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
 256       AliDebug(2, Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
 257       return kFALSE;
 258     }
 259   }
 260
 261   if (fAODVertex) {
 262     Double_t vert[3]={0};
 263     ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
 264     if (TMath::Abs(vert[2]) > fZVertexCut) {
 265       AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
 266       return kFALSE;
 267     }
 268   }
 269
 270   return kTRUE;
 271 }
 272
 273 //________________________________________________________________________
 274 void AliJetEmbeddingFromAODTask::Run()
 275 {
 276   if (!GetNextEntry()) {
 277     AliError("Unable to get AOD event to embed. Nothing will be embedded.");
 278     return;
 279   }
 280
 281   if (fTracks && fOutTracks) {
 282
 283     if (fCopyArray)
 284       CopyTracks();
 285
 286     if (fAODTracks) {
 287       AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
 288       for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
 289         AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
 290         if (!track) {
 291           AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
 292           continue;
 293         }
 294
 295         Int_t type = 0;
 296         Bool_t isEmc = kFALSE;
 297         if (!fEsdTreeMode) {
 298           AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
 299           if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
 300             type = 0;
 301           }
 302           else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
 303             if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
 304               if (fIncludeNoITS)
 305                 type = 2;
 306               else
 307                 continue;
 308             }
 309             else {
 310               type = 1;
 311             }
 312           }
 313           else { /*not a good track*/
 314             continue;
 315           }
 316
 317           if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
 318               aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
 319               aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
 320             isEmc = kTRUE;
 321         }
 322         else { /*not AOD mode, let's see if it is PicoTrack*/
 323           AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
 324           if (ptrack) {
 325             type = ptrack->GetTrackType();
 326             isEmc = ptrack->IsEMCAL();
 327           }
 328         }
 329
 330         AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f", track->Pt(), track->Eta(), track->Phi()));
 331         AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc);
 332       }
 333     }
 334   }
 335
 336   if (fClusters && fOutClusters) {
 337
 338     if (fCopyArray)
 339       CopyClusters();
 340
 341     if (fAODClusters) {
 342       for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
 343         AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
 344         if (!clus) {
 345           AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
 346           continue;
 347         }
 348         TLorentzVector vect;
 349         Double_t vert[3] = {0,0,0};
 350         clus->GetMomentum(vect,vert);
 351         AddCluster(clus->E(), vect.Eta(), vect.Phi());
 352       }
 353     }
 354   }
 355
 356   if (fCaloCells && fOutCaloCells) {
 357
 358     Int_t totalCells = fCaloCells->GetNumberOfCells();
 359     if (fAODCaloCells)
 360       totalCells += fAODCaloCells->GetNumberOfCells();
 361
 362     SetNumberOfOutCells(totalCells);
 363
 364     CopyCells();
 365
 366     if (fAODCaloCells) {
 367       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
 368         Short_t mclabel = -1;
 369         Double_t efrac = 0.;
 370         Double_t time = -1;
 371         Short_t cellNum = -1;
 372         Double_t amp = -1;
 373
 374         fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 375         AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
 376         AddCell(amp, cellNum, time);
 377       }
 378     }
 379   }
 380 }