HLT/trigger/AliHLTTriggerFastJet.cxx

   1 #include "AliHLTTriggerFastJet.h"
   2 #include "AliESDEvent.h"
   3 #include "AliESDtrackCuts.h"
   4 #include "AliESDCaloCluster.h"
   5 #include "AliHLTCaloClusterReader.h"
   6 #include "AliHLTCaloClusterDataStruct.h"
   7 #include "AliESDVZERO.h"
   8 #include "AliHLTTriggerDecision.h"
   9 #include "AliHLTDomainEntry.h"
  10 #include "TLorentzVector.h"
  11 #include "TRefArray.h"
  12 #include "TString.h"
  13 #include "TMap.h"
  14 #include "TObjArray.h"
  15 #include "TArrayI.h"
  16
  17 #include "AliFJWrapper.h"
  18
  19
  20 ClassImp(AliHLTTriggerFastJet)
  21
  22 AliHLTTriggerFastJet::AliHLTTriggerFastJet() :
  23 AliHLTTrigger(),
  24   fEThreshold(0.0),
  25   fDetector("EMCAL"),
  26   fFastJetWrapper(NULL),
  27   EsdTrackCuts(NULL),
  28   fOCDBEntry("HLT/ConfigHLT/EmcalJetTrigger")
  29 {
  30   fFastJetWrapper = new AliFJWrapper("FastJet","FastJet");
  31
  32 }
  33 //_____________________________________________________________
  34
  35 AliHLTTriggerFastJet::~AliHLTTriggerFastJet() {
  36
  37   if (fFastJetWrapper)
  38     delete fFastJetWrapper;
  39   fFastJetWrapper = NULL;
  40
  41   if (EsdTrackCuts)
  42     delete EsdTrackCuts;
  43   EsdTrackCuts = NULL;
  44
  45 }
  46 //_____________________________________________________________
  47
  48 int AliHLTTriggerFastJet::DoInit(int argc, const char** argv) {
  49
  50   int iResult = ConfigureFromCDBTObjString(fOCDBEntry);
  51
  52   if ( iResult>=0 && argc>0 ) {
  53     iResult = ConfigureFromArgumentString(argc, argv);
  54     HLTImportant("Trigger threshold set from argument string :  %.02f GeV:", fEThreshold );
  55   } else if ( iResult>=0 ) {
  56     HLTImportant("Trigger threshold set from OCDB database entry : %.02f Gev:", fEThreshold );
  57   }
  58   return iResult;
  59 }
  60 //______________________________________________________________
  61
  62 int AliHLTTriggerFastJet::DoDeInit() {
  63
  64   return 0;
  65
  66 }
  67 //______________________________________________________________
  68
  69
  70
  71 AliHLTComponent* AliHLTTriggerFastJet::Spawn() {
  72   // see header file for class documentation
  73   return new AliHLTTriggerFastJet;
  74 }
  75
  76
  77 const char* AliHLTTriggerFastJet::GetTriggerName() const {
  78   // see header file for class documentation
  79   return "EmcalJetTrigger";
  80 }
  81
  82 Int_t AliHLTTriggerFastJet::DoTrigger() {
  83
  84   Int_t iResult = 0;
  85
  86   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
  87     return 0;
  88
  89   const TObject *obj = GetFirstInputObject( kAliHLTAllDataTypes, "AliESDEvent" );
  90   AliESDEvent   *esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
  91
  92   // check for MatchedTrack to avoid double counting
  93
  94
  95   if ( esd != NULL ) {
  96     esd->GetStdContent();
  97
  98     //Double_t vertex[3] = {0,0,0};
  99     //esd->GetVertex()->GetXYZ(vertex);
 100     //cout << "Vertex " << vertex[0] << vertex[1] << vertex [2] << endl;
 101     //TLorentzVector gamma;
 102
 103     // -- add VZERO
 104     // uncomment for usage (to avoid warning)
 105
 106     /*
 107     AliESDVZERO *v0  = esd->GetVZEROData();
 108
 109     Float_t MtotVOA = v0->GetMTotV0A();
 110     Float_t MtotVOC = v0->GetMTotV0A();
 111     */
 112     Int_t nTracks = esd->GetNumberOfTracks();
 113
 114     TObjArray *tracks = new TObjArray(nTracks);
 115     tracks->SetOwner(1);
 116
 117
 118     EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 119     EsdTrackCuts->SetMinNClustersTPC(70);
 120
 121     // add the tracks to FastJet Wrapper
 122
 123     for ( Int_t i = 0; i < nTracks; i++ ) {
 124
 125       const AliESDVertex *vtxTPC = esd->GetPrimaryVertexTPC();
 126       if (!vtxTPC) continue;
 127       cout << "Got primary vertex " << endl;
 128
 129       AliESDtrack *esdtrack = esd->GetTrack(i);
 130       if (!esdtrack) continue;
 131       tracks->Add(esdtrack);
 132       cout << "Got esd track " << endl;
 133
 134       AliESDtrack *tpctrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdtrack->GetID());
 135       if (!tpctrack) continue;
 136       cout << "Got tpc track " << endl;
 137
 138       AliExternalTrackParam exParam;
 139       Bool_t relate = tpctrack->RelateToVertexTPC(vtxTPC,esd->GetMagneticField(),kVeryBig,&exParam);
 140       if (!relate) {
 141         delete tpctrack;
 142         continue;
 143       }
 144
 145       tpctrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 146       cout << "Set param to tpctrack " << endl;
 147
 148       fFastJetWrapper->AddInputVector(tpctrack->Px(),tpctrack->Py(),tpctrack->Pz(),tpctrack->P());
 149       cout << "Added track with P: " << tpctrack->P() << endl;
 150     }
 151
 152
 153 //    // add the Calo Clusters to FastJet Wrapper
 154 //    TRefArray *caloClustersArr = new TRefArray();
 155 //    esd->GetEMCALClusters(caloClustersArr);
 156 //    const Int_t nEMCALClusters = caloClustersArr->GetEntries();
 157 //
 158 //    for ( UInt_t ic = 0; ic < nEMCALClusters; ic++ )
 159 //      {
 160 //      AliESDCaloCluster *cluster = dynamic_cast<AliESDCaloCluster*>( caloClustersArr->At(ic) );
 161 //        Double_t clusterEn = cluster->E();
 162 //      Double_t subE = 0;
 163 //
 164 //        TArrayI *matched = cluster->GetTracksMatched();
 165 //      for ( UInt_t jk = 0; jk < matched->GetSize(); jk++ )
 166 //          {
 167 //            if ( matched->At(jk) > -1 && matched->At(jk) < esd ->GetNumberOfTracks() )
 168 //              {
 169 //                AliESDtrack *track = dynamic_cast<AliESDtrack*>( tracks->At( (matched->At(jk)) ) );
 170 //                if ( !track ) continue;
 171 //                subE += track->P();
 172 //              }
 173 //          }
 174 //        clusterEn -= subE;
 175 //        if ( clusterEn < 0)
 176 //          clusterEn = 0;
 177 //        cluster->SetE( clusterEn );
 178 //        if ( cluster->E() == 0 ) continue;
 179 //
 180 //        cluster->GetMomentum( gamma, vertex );
 181 //        fFastJetWrapper->AddInputVector( gamma.Px(), gamma.Py(), gamma.Pz(), clusterEn );
 182 //        gamma.Clear();
 183 //      }
 184 //  }
 185     for ( Int_t j = 0; j< esd->GetNumberOfCaloClusters(); j++) {
 186       //AliESDCaloCluster *cluster = static_cast<AliESDCaloCluster*>(esd->GetCaloCluster(j));
 187       AliVCluster *cluster = dynamic_cast<AliVCluster*>(esd->GetCaloCluster(j));
 188       if ( cluster->IsEMCAL() ) {
 189         Int_t trackindex = cluster->GetTrackMatchedIndex();
 190         if (trackindex<0) continue;
 191         AliESDtrack *track = esd->GetTrack(trackindex);
 192         if (!track) continue;
 193         fFastJetWrapper->AddInputVector(track->Px(),track->Py(),track->Pz(),cluster->E());
 194         cout << "Added cluster with E: " << cluster->E() << endl;
 195       }
 196     }
 197   }
 198
 199   fFastJetWrapper->SetupStrategyfromOpt("Best");
 200   fFastJetWrapper->SetupAlgorithmfromOpt("antikt");
 201   fFastJetWrapper->SetupSchemefromOpt("BIpt");
 202   fFastJetWrapper->SetupAreaTypefromOpt("active");
 203   fFastJetWrapper->SetNRepeats(1);
 204   fFastJetWrapper->SetGhostArea(0.01);
 205   fFastJetWrapper->SetMaxRap(0.9);
 206   fFastJetWrapper->SetR(0.4);
 207
 208   double median = 0.0;
 209   double sigma  = 0.0;
 210
 211   fFastJetWrapper->Run();
 212   fFastJetWrapper->GetMedianAndSigma(median,sigma);
 213
 214   if ( TriggerOnJet(fFastJetWrapper->GetSubtractedJetsPts(median)) ) {
 215     return iResult;
 216   }
 217
 218   TString description;
 219   description.Form(" No jets with energy >  %.02f GeV found! ", fDetector.Data(), fEThreshold);
 220   SetDescription(description.Data());
 221   TriggerEvent(kFALSE);
 222
 223   return iResult;
 224 }
 225
 226 //______________________________________________________________
 227
 228 template <class T>
 229 Bool_t AliHLTTriggerFastJet::TriggerOnJet(T Jet) {
 230
 231   for ( unsigned int ij = 0; ij<Jet.size(); ij++ ) {
 232     if ( Jet[ij] > fEThreshold ) {
 233       TString description;
 234       description.Form(" Event contains at least one %s jet with energy > %.02f GeV! ", fDetector.Data(), fEThreshold);
 235       SetDescription(description.Data());
 236
 237       // Enable the detectors for readout
 238       GetReadoutList().Enable(AliHLTReadoutList::kEMCAL |
 239                               AliHLTReadoutList::kTPC);
 240
 241       // Add the available HLT info for readout
 242       GetTriggerDomain().Add(kAliHLTAnyDataTypeID, fDetector.Data());
 243
 244       // Set trigger desicion
 245       TriggerEvent(kTRUE);
 246
 247       return(kTRUE);
 248     }
 249   }
 250
 251   return kFALSE;
 252
 253 }
 254 //______________________________________________________________
 255
 256 int AliHLTTriggerFastJet::Reconfigure(const char* cdbEntry, const char* /*chainId*/) {
 257
 258   const char* entry = cdbEntry;
 259   if ( !entry || entry[0]==0 ) entry=fOCDBEntry;
 260
 261   return ConfigureFromCDBTObjString(entry);
 262
 263 }
 264
 265 //______________________________________________________________
 266
 267 int AliHLTTriggerFastJet::ScanConfigurationArgument(int argc, const char** argv) {
 268
 269   if ( argc<=0 ) return 0;
 270   int i = 0;
 271   TString argument = argv[i];
 272
 273   if ( argument.CompareTo("-energy") == 0 ) {
 274     if (++i>=argc) return -EPROTO;
 275     argument = argv[i];
 276     fEThreshold=argument.Atof();
 277     return 2;
 278   }
 279
 280   return -EINVAL;
 281
 282 }
 283 //______________________________________________________________
 284
 285 void AliHLTTriggerFastJet::GetOutputDataSize(unsigned long &constBase, double &inputMultiplier) {
 286
 287   constBase = sizeof(AliHLTTriggerDecision) + sizeof(AliHLTDomainEntry)*14;
 288   inputMultiplier = 1;
 289
 290 }
 291 //______________________________________________________________
 292
 293 void AliHLTTriggerFastJet::GetOCDBObjectDescription(TMap* const targetMap) {
 294
 295   if ( !targetMap ) return;
 296   targetMap->Add(new TObjString(fOCDBEntry),
 297                  new TObjString(Form("%s threshold trigger OCDB object", fDetector.Data()) )
 298                  );
 299
 300 }
 301 //_______________________________________________________________