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