Updated macros for PHOS alignment calculation
[u/mrichter/AliRoot.git] / 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 //_______________________________________________________________