HLT/trigger/AliHLTTriggerFastJet.cxx

   1 //**************************************************************************
   2 //* This file is property of and copyright by the ALICE HLT Project        *
   3 //* ALICE Experiment at CERN, All rights reserved.                         *
   4 //*                                                                        *
   5 //* Primary Authors: leonidas.xaplanteris@gmail.com                        *
   6 //*                  for The ALICE HLT Project.                            *
   7 //*                                                                        *
   8 //* Permission to use, copy, modify and distribute this software and its   *
   9 //* documentation strictly for non-commercial purposes is hereby granted   *
  10 //* without fee, provided that the above copyright notice appears in all   *
  11 //* copies and that both the copyright notice and this permission notice   *
  12 //* appear in the supporting documentation. The authors make no claims     *
  13 //* about the suitability of this software for any purpose. It is          *
  14 //* provided "as is" without express or implied warranty.                  *
  15 //**************************************************************************
  16
  17 #include "AliHLTTriggerFastJet.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDtrackCuts.h"
  20 #include "AliESDCaloCluster.h"
  21 #include "AliHLTCaloClusterReader.h"
  22 #include "AliHLTCaloClusterDataStruct.h"
  23 #include "AliESDVZERO.h"
  24 #include "AliHLTTriggerDecision.h"
  25 #include "AliHLTDomainEntry.h"
  26 #include "TLorentzVector.h"
  27 #include "TRefArray.h"
  28 #include "TString.h"
  29 #include "TMap.h"
  30 #include "TObjArray.h"
  31 #include "TArrayI.h"
  32 #include "AliHLTScalars.h"
  33
  34 #include "AliFJWrapper.h"
  35
  36
  37 ClassImp(AliHLTTriggerFastJet)
  38
  39 AliHLTTriggerFastJet::AliHLTTriggerFastJet() :
  40 AliHLTTrigger(),
  41   fEThreshold(0.0),
  42   fDetector("EMCAL"),
  43   fFastJetWrapper(NULL),
  44   fMakeStats(kFALSE),
  45   EsdTrackCuts(NULL),
  46   fOCDBEntry("HLT/ConfigHLT/EmcalJetTrigger")
  47 {
  48   fFastJetWrapper = new AliFJWrapper("FastJet","FastJet");
  49
  50   // check if the fMakeStats flag is set to use AliHLTScalars
  51   // to Add a new scalar:
  52   // scalars.Add(const char *name, const char *description, Double_t value)
  53
  54   //  if ( fMakeStats ) AliHLTScalars scalars;
  55
  56 }
  57 //_____________________________________________________________
  58
  59 AliHLTTriggerFastJet::~AliHLTTriggerFastJet() {
  60
  61   if (fFastJetWrapper)
  62     delete fFastJetWrapper;
  63   fFastJetWrapper = NULL;
  64
  65   if (EsdTrackCuts)
  66     delete EsdTrackCuts;
  67   EsdTrackCuts = NULL;
  68
  69 }
  70 //_____________________________________________________________
  71
  72 int AliHLTTriggerFastJet::DoInit(int argc, const char** argv) {
  73
  74   int iResult = ConfigureFromCDBTObjString(fOCDBEntry);
  75
  76   if ( iResult>=0 && argc>0 ) {
  77     iResult = ConfigureFromArgumentString(argc, argv);
  78     HLTImportant("Trigger threshold set from argument string :  %.02f GeV:", fEThreshold );
  79   } else if ( iResult>=0 ) {
  80     HLTImportant("Trigger threshold set from OCDB database entry : %.02f Gev:", fEThreshold );
  81   }
  82   return iResult;
  83 }
  84 //______________________________________________________________
  85
  86 int AliHLTTriggerFastJet::DoDeInit() {
  87
  88   return 0;
  89
  90 }
  91 //______________________________________________________________
  92
  93
  94
  95 AliHLTComponent* AliHLTTriggerFastJet::Spawn() {
  96   // see header file for class documentation
  97   return new AliHLTTriggerFastJet;
  98 }
  99
 100
 101 const char* AliHLTTriggerFastJet::GetTriggerName() const {
 102   // see header file for class documentation
 103   return "EmcalJetTrigger";
 104 }
 105
 106 Int_t AliHLTTriggerFastJet::DoTrigger() {
 107
 108   Int_t iResult = 0;
 109   AliHLTScalars scalars;
 110
 111   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
 112     return 0;
 113
 114   const TObject *obj = GetFirstInputObject( kAliHLTAllDataTypes, "AliESDEvent" );
 115   AliESDEvent   *esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
 116
 117   // check for MatchedTrack to avoid double counting
 118
 119
 120
 121   if ( esd != NULL ) {
 122     esd->GetStdContent();
 123
 124     Double_t vertex[3] = {0,0,0};
 125     vertex[0] = esd->GetVertex()->GetX();
 126     vertex[1] = esd->GetVertex()->GetY();
 127     vertex[2] = esd->GetVertex()->GetZ();
 128     //cout << "Vertex " << vertex[0] << vertex[1] << vertex [2] << endl;
 129     //TLorentzVector gamma;
 130
 131     // -- add VZERO
 132     // uncomment for usage (to avoid warning)
 133
 134     /*
 135     AliESDVZERO *v0  = esd->GetVZEROData();
 136
 137     Float_t MtotVOA = v0->GetMTotV0A();
 138     Float_t MtotVOC = v0->GetMTotV0A();
 139     */
 140     Int_t nTracks = esd->GetNumberOfTracks();
 141
 142     TObjArray *tracks = new TObjArray(nTracks);
 143     tracks->SetOwner(1);
 144
 145
 146     EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 147     EsdTrackCuts->SetMinNClustersTPC(70);
 148
 149     // add the tracks to FastJet Wrapper
 150
 151     for ( Int_t i = 0; i < nTracks; i++ ) {
 152
 153       const AliESDVertex *vtxTPC = esd->GetPrimaryVertexTPC();
 154       if (!vtxTPC) continue;
 155       //cout << "Got primary vertex " << endl;
 156
 157       AliESDtrack *esdtrack = esd->GetTrack(i);
 158       if (!esdtrack) continue;
 159       tracks->Add(esdtrack);
 160       //cout << "Got esd track " << endl;
 161
 162       AliESDtrack *tpctrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdtrack->GetID());
 163       if (!tpctrack) continue;
 164       //cout << "Got tpc track " << endl;
 165
 166       AliExternalTrackParam exParam;
 167       Bool_t relate = tpctrack->RelateToVertexTPC(vtxTPC,esd->GetMagneticField(),kVeryBig,&exParam);
 168       if (!relate) {
 169         delete tpctrack;
 170         continue;
 171       }
 172
 173       tpctrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 174       //cout << "Set param to tpctrack " << endl;
 175
 176       fFastJetWrapper->AddInputVector(tpctrack->Px(),tpctrack->Py(),tpctrack->Pz(),tpctrack->P());
 177       if ( fMakeStats ) scalars.Add("TracksPt","TPC tracks pT", tpctrack->Pt());
 178       //cout << "Added track with P: " << tpctrack->P() << endl;
 179     }
 180
 181
 182 //    // add the Calo Clusters to FastJet Wrapper
 183 //    TRefArray *caloClustersArr = new TRefArray();
 184 //    esd->GetEMCALClusters(caloClustersArr);
 185 //    const Int_t nEMCALClusters = caloClustersArr->GetEntries();
 186 //
 187 //    for ( UInt_t ic = 0; ic < nEMCALClusters; ic++ )
 188 //      {
 189 //      AliESDCaloCluster *cluster = dynamic_cast<AliESDCaloCluster*>( caloClustersArr->At(ic) );
 190 //        Double_t clusterEn = cluster->E();
 191 //      Double_t subE = 0;
 192 //
 193 //        TArrayI *matched = cluster->GetTracksMatched();
 194 //      for ( UInt_t jk = 0; jk < matched->GetSize(); jk++ )
 195 //          {
 196 //            if ( matched->At(jk) > -1 && matched->At(jk) < esd ->GetNumberOfTracks() )
 197 //              {
 198 //                AliESDtrack *track = dynamic_cast<AliESDtrack*>( tracks->At( (matched->At(jk)) ) );
 199 //                if ( !track ) continue;
 200 //                subE += track->P();
 201 //              }
 202 //          }
 203 //        clusterEn -= subE;
 204 //        if ( clusterEn < 0)
 205 //          clusterEn = 0;
 206 //        cluster->SetE( clusterEn );
 207 //        if ( cluster->E() == 0 ) continue;
 208 //
 209 //        cluster->GetMomentum( gamma, vertex );
 210 //        fFastJetWrapper->AddInputVector( gamma.Px(), gamma.Py(), gamma.Pz(), clusterEn );
 211 //        if ( fMakeStats ) {
 212 //           scalars.Add("ClusterEn","Clusters Energy",clusterEn);
 213 //           scalars.Add("ClusterEta","Clusters Eta", gamma.Eta());
 214 //           scalars.Add("ClusterPhi","Clusters Phi", gamma.Phi());
 215 //        }
 216 //        gamma.Clear();
 217 //      }
 218 //  }
 219     for ( Int_t j = 0; j< esd->GetNumberOfCaloClusters(); j++) {
 220       //AliESDCaloCluster *cluster = static_cast<AliESDCaloCluster*>(esd->GetCaloCluster(j));
 221       AliVCluster *cluster = dynamic_cast<AliVCluster*>(esd->GetCaloCluster(j));
 222       if ( cluster->IsEMCAL() ) {
 223         Int_t trackindex = cluster->GetTrackMatchedIndex();
 224         if (trackindex<0) continue;
 225         AliESDtrack *track = esd->GetTrack(trackindex);
 226         if (!track) continue;
 227         fFastJetWrapper->AddInputVector(track->Px(),track->Py(),track->Pz(),cluster->E());
 228         //cout << "Added cluster with E: " << cluster->E() << endl;
 229         if ( fMakeStats ) {
 230           scalars.Add("ClusterEn" , "Clusters Energy", cluster->E());
 231           scalars.Add("ClusterEta", "Clusters Eta"   , track->Eta());
 232           scalars.Add("ClusterPhi", "Clusters Phi"   , track->Phi());
 233         }
 234       }
 235     }
 236   }
 237
 238   fFastJetWrapper->SetupStrategyfromOpt("Best");
 239   fFastJetWrapper->SetupAlgorithmfromOpt("antikt");
 240   fFastJetWrapper->SetupSchemefromOpt("BIpt");
 241   fFastJetWrapper->SetupAreaTypefromOpt("active");
 242   fFastJetWrapper->SetNRepeats(1);
 243   fFastJetWrapper->SetGhostArea(0.01);
 244   fFastJetWrapper->SetMaxRap(0.9);
 245   fFastJetWrapper->SetR(0.4);
 246
 247   double median = 0.0;
 248   double sigma  = 0.0;
 249
 250   fFastJetWrapper->Run();
 251   fFastJetWrapper->GetMedianAndSigma(median,sigma);
 252
 253   // checks if trigger criteria are fulfilled
 254   if ( TriggerOnJet(fFastJetWrapper->GetSubtractedJetsPts(median)) ) {
 255     if ( fMakeStats ) iResult = PushBack(&scalars, kAliHLTDataTypeEventStatistics|kAliHLTDataOriginHLT);
 256     return iResult;
 257   }
 258
 259   // if we got to this point it meens the trigger criteria were not fulfilled
 260   TString description;
 261   description.Form(" No jets with energy >  %.02f GeV found! ", fDetector.Data(), fEThreshold);
 262   SetDescription(description.Data());
 263   TriggerEvent(kFALSE);
 264
 265   if ( fMakeStats ) iResult = PushBack(&scalars, kAliHLTDataTypeEventStatistics|kAliHLTDataOriginHLT);
 266   return iResult;
 267 }
 268
 269 //______________________________________________________________
 270
 271 template <class T>
 272 Bool_t AliHLTTriggerFastJet::TriggerOnJet(T Jet) {
 273
 274   AliHLTScalars scalars;
 275
 276   for ( unsigned int ij = 0; ij<Jet.size(); ij++ ) {
 277     if ( Jet[ij] > fEThreshold ) {
 278       TString description;
 279       description.Form(" Event contains at least one %s jet with energy > %.02f GeV! ", fDetector.Data(), fEThreshold);
 280       SetDescription(description.Data());
 281
 282       // Enable the detectors for readout
 283       GetReadoutList().Enable(AliHLTReadoutList::kEMCAL |
 284                               AliHLTReadoutList::kTPC);
 285
 286       // Add the available HLT info for readout
 287       GetTriggerDomain().Add(kAliHLTAnyDataTypeID, fDetector.Data());
 288
 289       // Set trigger desicion
 290       TriggerEvent(kTRUE);
 291
 292       // if fMakeStats true fill stats
 293       if ( fMakeStats ) scalars.Add("JetsPt","Jets pT", Jet[ij]);
 294
 295       return(kTRUE);
 296     }
 297   }
 298
 299   return kFALSE;
 300
 301 }
 302 //______________________________________________________________
 303
 304 int AliHLTTriggerFastJet::Reconfigure(const char* cdbEntry, const char* /*chainId*/) {
 305
 306   const char* entry = cdbEntry;
 307   if ( !entry || entry[0]==0 ) entry=fOCDBEntry;
 308
 309   return ConfigureFromCDBTObjString(entry);
 310
 311 }
 312
 313 //______________________________________________________________
 314
 315 int AliHLTTriggerFastJet::ScanConfigurationArgument(int argc, const char** argv) {
 316
 317   if ( argc<=0 ) return 0;
 318   int i = 0;
 319   TString argument = argv[i];
 320
 321   if ( argument.CompareTo("-energy") == 0 ) {
 322     if (++i>=argc) return -EPROTO;
 323     argument = argv[i];
 324     fEThreshold=argument.Atof();
 325     return 2;
 326   }
 327
 328   if ( argument.CompareTo("-makestats") == 0) {
 329     fMakeStats = kTRUE;
 330     return 2;
 331   }
 332
 333   return -EINVAL;
 334
 335 }
 336 //______________________________________________________________
 337
 338 void AliHLTTriggerFastJet::GetOutputDataTypes(AliHLTComponentDataTypeList &list) const {
 339   // return the output data types generated
 340
 341   list.push_back(kAliHLTDataTypeTriggerDecision);
 342   list.push_back(kAliHLTDataTypeEventStatistics|kAliHLTDataOriginHLT);
 343 }
 344 //_______________________________________________________________
 345
 346 void AliHLTTriggerFastJet::GetOutputDataSize(unsigned long &constBase, double &inputMultiplier) {
 347
 348   constBase = sizeof(AliHLTTriggerDecision) + sizeof(AliHLTDomainEntry)*14;
 349   inputMultiplier = 1;
 350
 351 }
 352 //______________________________________________________________
 353
 354 void AliHLTTriggerFastJet::GetOCDBObjectDescription(TMap* const targetMap) {
 355
 356   if ( !targetMap ) return;
 357   targetMap->Add(new TObjString(fOCDBEntry),
 358                  new TObjString(Form("%s threshold trigger OCDB object", fDetector.Data()) )
 359                  );
 360
 361 }
 362 //_______________________________________________________________