]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/trigger/AliHLTD0Trigger.cxx
using ITS Tracks and Vertex from global vertexer directly, some further development...
[u/mrichter/AliRoot.git] / HLT / trigger / AliHLTD0Trigger.cxx
CommitLineData
5f4502cc 1// $Id: AliHLTD0Trigger.cxx
2//**************************************************************************
3//* This file is property of and copyright by the ALICE HLT Project *
4//* ALICE Experiment at CERN, All rights reserved. *
5//* *
6//* Primary Authors: Gaute Ovrebekk *
7//* for The ALICE HLT Project. *
8//* *
9//* Permission to use, copy, modify and distribute this software and its *
10//* documentation strictly for non-commercial purposes is hereby granted *
11//* without fee, provided that the above copyright notice appears in all *
12//* copies and that both the copyright notice and this permission notice *
13//* appear in the supporting documentation. The authors make no claims *
14//* about the suitability of this software for any purpose. It is *
15//* provided "as is" without express or implied warranty. *
16//**************************************************************************
17
18/// @file AliHLTD0Trigger.cxx
19/// @author Gaute Ovrebekk
20/// @date 2009-10-28
21/// @brief HLT trigger component for D0->Kpi
22
23// see header file for class documentation
24// or
25// refer to README to build package
26// or
27// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
28
29#include "AliHLTD0Trigger.h"
30#include "AliESDEvent.h"
5f4502cc 31#include "AliESDv0.h"
32#include "AliHLTTriggerDecision.h"
33#include "AliHLTDomainEntry.h"
34#include "AliHLTGlobalBarrelTrack.h"
35#include "TObjArray.h"
36#include "TObjString.h"
37#include "TDatabasePDG.h"
38#include "AliESDVertex.h"
39#include "TH1F.h"
40#include "AliHLTD0toKpi.h"
9929f8f5 41#include "AliAODVertex.h"
629b904b 42#include "AliESDVertex.h"
3a626a12 43#include "AliAODRecoDecay.h"
5f4502cc 44
45/** ROOT macro for the implementation of ROOT specific class methods */
46ClassImp(AliHLTD0Trigger)
47
48AliHLTD0Trigger::AliHLTD0Trigger()
49 : AliHLTTrigger()
50 , fPtMin(0.0)
51 , fdca(0.0)
52 , finvMass(0.0)
53 , fcosThetaStar(0.0)
54 , fd0(0.0)
55 , fd0d0(0.0)
56 , fcosPoint(0.0)
9929f8f5 57 , fplothisto(false)
58 , fUseV0(false)
59 , mD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
5f4502cc 60 , fD0mass(NULL)
9929f8f5 61 , fPos()
62 , fNeg()
63 , fd0calc(NULL)
64 , ftwoTrackArray(NULL)
65 , fTotalD0(0)
629b904b 66 , fVertex(NULL)
67 , fField(0)
5f4502cc 68{
9929f8f5 69
5f4502cc 70 // see header file for class documentation
71 // or
72 // refer to README to build package
73 // or
74 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
75
76}
77
78const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
79
80AliHLTD0Trigger::~AliHLTD0Trigger()
81{
9929f8f5 82 //if(fd0calc){delete fd0calc;}
83 //if(fD0mass){delete fD0mass;}
84 //if(ftwoTrackArray){delete ftwoTrackArray;}
5f4502cc 85 // see header file for class documentation
5f4502cc 86}
9929f8f5 87
5f4502cc 88const char* AliHLTD0Trigger::GetTriggerName() const
89{
90 // see header file for class documentation
91 return "D0Trigger";
92}
93
94AliHLTComponent* AliHLTD0Trigger::Spawn()
95{
96 // see header file for class documentation
97 return new AliHLTD0Trigger;
98}
99
100int AliHLTD0Trigger::DoTrigger()
101{
9929f8f5 102 // -- Iterator over Data Blocks --
103 //const AliHLTComponentBlockData* iter = NULL;
5f4502cc 104
9929f8f5 105 if (!IsDataEvent()) return 0;
5f4502cc 106
9929f8f5 107 Int_t nD0=0;
108 TString description;
5f4502cc 109
9929f8f5 110 for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
111 if(fUseV0){
112 nD0=RecV0(iter);
113 }
114 else{
629b904b 115 AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
116 event->GetStdContent();
117 fField = event->GetMagneticField();
118 const AliESDVertex* pv = event->GetPrimaryVertexTracks();
119 fVertex = new AliESDVertex(*pv);
3a626a12 120 if(fVertex->GetNContributors()<2){
121 HLTWarning("Contributors in ESD vertex to low or not been set");
122 continue;
123 }
629b904b 124 for(Int_t it=0;it<event->GetNumberOfTracks();it++){
125 SingleTrackSelect(event->GetTrack(it));
126 }
127
128 nD0=RecD0();
5f4502cc 129 }
5f4502cc 130 }
9929f8f5 131
3a626a12 132 for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut);
629b904b 133 iter != NULL; iter = GetNextInputObject() ) {
134 fVertex = dynamic_cast<AliESDVertex*>(const_cast<TObject*>( iter ));
135 if(!fVertex){
3a626a12 136 HLTError("Vertex object is corrupted");
629b904b 137 //iResult = -EINVAL;
138 }
139 }
140
141 for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
142 iter != NULL; iter = GetNextInputBlock() ) {
143 vector<AliHLTGlobalBarrelTrack> tracksVector;
144 AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
145
146 fField = GetBz();
147 for(UInt_t i=0;i<tracksVector.size();i++){
148 SingleTrackSelect(&tracksVector[i]);
149 }
150 nD0=RecD0();
9929f8f5 151 }
152
153 fTotalD0+=nD0;
154
155 ftwoTrackArray->Clear();
156 fPos.clear();
157 fNeg.clear();
629b904b 158
159 HLTInfo("Number of D0 found: %d",nD0);
160 HLTInfo("Total Number of D0 found: %d",fTotalD0);
5f4502cc 161
162 if(fplothisto){PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);}
163
9929f8f5 164 //if (iResult>=0) {
165 if (1) {
5f4502cc 166
167 if (nD0>=1) {
168 description.Form("Event contains %d D0(s)", nD0);
169 SetDescription(description.Data());
170 // Enable the central detectors for readout.
171 GetReadoutList().Enable(
172 AliHLTReadoutList::kITSSPD |
173 AliHLTReadoutList::kITSSDD |
174 AliHLTReadoutList::kITSSSD |
175 AliHLTReadoutList::kTPC |
176 AliHLTReadoutList::kTRD |
177 AliHLTReadoutList::kTOF |
178 AliHLTReadoutList::kHMPID |
179 AliHLTReadoutList::kPHOS
180 );
181 // Add the available HLT information for readout too.
182 GetTriggerDomain().Add("CLUSTERS", "TPC ");
183 TriggerEvent(true);
184 return 0;
185 }
186 description.Form("No D0");
187 } else {
188 description.Form("No input blocks found");
189 }
190 SetDescription(description.Data());
191 TriggerEvent(false);
9929f8f5 192 //return iResult;
193 return 0;
5f4502cc 194}
195
196int AliHLTD0Trigger::DoInit(int argc, const char** argv)
197{
9929f8f5 198
199 fd0calc = new AliHLTD0toKpi();
200 ftwoTrackArray = new TObjArray(2);
201
5f4502cc 202 fplothisto=false;
203 // see header file for class documentation
9929f8f5 204 fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
5f4502cc 205 // first configure the default
206 int iResult=0;
207 if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
208
209 // configure from the command line parameters if specified
210 if (iResult>=0 && argc>0)
211 iResult=ConfigureFromArgumentString(argc, argv);
212 return iResult;
213}
214
215int AliHLTD0Trigger::DoDeinit()
216{
217 // see header file for class documentation
9929f8f5 218 if(fd0calc){delete fd0calc;}
219 if(fD0mass){delete fD0mass;}
220 if(ftwoTrackArray){delete ftwoTrackArray;}
629b904b 221 if(fVertex){delete fVertex;}
5f4502cc 222 return 0;
223}
224
225int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
226{
227 // see header file for class documentation
228
229 // configure from the specified antry or the default one
230 const char* entry=cdbEntry;
231 if (!entry || entry[0]==0) {
232 entry=fgkOCDBEntry;
233 }
234
235 return ConfigureFromCDBTObjString(entry);
236}
237
238int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
239{
240 // see header file for class documentation
241 if (argc<=0) return 0;
242 int i=0;
243 TString argument=argv[i];
244 // -minpt for decay
245 if (argument.CompareTo("-pt")==0) {
246 if (++i>=argc) return -EPROTO;
247 argument=argv[i];
248 fPtMin=argument.Atof();
249 return 2;
250 }
251 // minimum dca for decay tracks
252 if (argument.CompareTo("-dca")==0) {
253 if (++i>=argc) return -EPROTO;
254 argument=argv[i];
255 fdca=argument.Atof();
256 return 2;
257 }
258 // inv. mass half width.
259 if (argument.CompareTo("-invmass")==0) {
260 if (++i>=argc) return -EPROTO;
261 argument=argv[i];
262 finvMass=argument.Atof();
263 return 2;
264 }
265
266// cos theta for decay angle
267 if (argument.CompareTo("-costhetastar")==0) {
268 if (++i>=argc) return -EPROTO;
269 argument=argv[i];
270 fcosThetaStar=argument.Atof();
271 return 2;
272 }
273
274 // impact parameter for decay
275 if (argument.CompareTo("-d0")==0) {
276 if (++i>=argc) return -EPROTO;
277 argument=argv[i];
278 fd0=argument.Atof();
279 return 2;
280 }
281 // product of impact parameter
282 if (argument.CompareTo("-d0d0")==0) {
283 if (++i>=argc) return -EPROTO;
284 argument=argv[i];
285 fd0d0=argument.Atof();
286 return 2;
287 }
288 // product of impact parameter
289 if (argument.CompareTo("-cospoint")==0) {
290 if (++i>=argc) return -EPROTO;
291 argument=argv[i];
292 fcosPoint=argument.Atof();
293 return 2;
294 }
295 if (argument.CompareTo("-plothistogram")==0) {
5f4502cc 296 fplothisto=true;
9929f8f5 297 return 1;
298 }
299 if (argument.CompareTo("-usev0")==0) {
300 fUseV0=true;
301 return 1;
5f4502cc 302 }
9929f8f5 303
5f4502cc 304 // unknown argument
305 return -EINVAL;
306}
9929f8f5 307
629b904b 308void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t){
9929f8f5 309 // Offline har || på disse kuttene på de to henfallsproduktene
629b904b 310 Double_t pv[3];
311 fVertex->GetXYZ(pv);
312
9929f8f5 313 if(t->Pt()<fPtMin){return;}
629b904b 314 if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > fd0){return;}
9929f8f5 315
316 if(t->Charge()>0){
317 fPos.push_back(t);
318 }
319 else{
320 fNeg.push_back(t);
321 }
322}
323
629b904b 324Int_t AliHLTD0Trigger::RecD0(){
9929f8f5 325
326 int nD0=0;
327 Double_t D0,D0bar,xdummy,ydummy;
328 Double_t d0[2];
329 Double_t svpos[3];
330 Double_t pvpos[3];
331
629b904b 332 if(!fVertex){
333 HLTError("No Vertex is set");
334 return 0;
9929f8f5 335 }
629b904b 336 fVertex->GetXYZ(pvpos);
337
338 for(UInt_t ip=0;ip<fPos.size();ip++){
339 AliExternalTrackParam *tP=fPos[ip];
340 for(UInt_t in=0;in<fNeg.size();in++){
341 AliExternalTrackParam *tN=fNeg[in];
342
343 tP->PropagateToDCA(fVertex,fField,kVeryBig); //do I need this??????
344 tN->PropagateToDCA(fVertex,fField,kVeryBig);
9929f8f5 345
629b904b 346 Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
9929f8f5 347 if(dcatPtN>fdca) { continue; }
348
629b904b 349 ftwoTrackArray->AddAt(tP,0);
350 ftwoTrackArray->AddAt(tN,1);
351 AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,fField,fVertex);
9929f8f5 352 if(!vertexp1n1) {
353 ftwoTrackArray->Clear();
354 continue;
355 }
356
357 vertexp1n1->GetXYZ(svpos);
358
629b904b 359 tP->PropagateToDCA(vertexp1n1,fField,kVeryBig);
360 tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
9929f8f5 361
3a626a12 362 /*
363 Double_t px[2],py[2],pz[2];
364 Double_t momentum[3];
365 tP->GetPxPyPz(momentum);
366 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
367 tN->GetPxPyPz(momentum);
368 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
369
370 Short_t dummycharge=0;
371 Double_t *dummyd0 = new Double_t[2];
372 Int_t nprongs = 2;
373
374 for(Int_t ipr=0;ipr<nprongs;ipr++) dummyd0[ipr]=0.;
375 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
376 delete [] dummyd0; dummyd0=NULL;
377
378 UInt_t pdg2[2],pdg2bar[2];
379 Double_t mPDG,minv,minvbar;
380
381 pdg2[0]=211; pdg2[1]=321; pdg2bar[0]=321; pdg2bar[1]=211;
382 minv = rd->InvMass(nprongs,pdg2);
383 minvbar = rd->InvMass(nprongs,pdg2bar);
384 if(TMath::Abs(minv-mD0PDG)>finvMass && TMath::Abs(minv-mD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
385 */
9929f8f5 386 if((TMath::Abs(fd0calc->InvMass(tN,tP)-mD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-mD0PDG) > finvMass){continue;}
387 fd0calc->cosThetaStar(tN,tP,D0,D0bar);
388 if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){continue;}
629b904b 389 d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
390 d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
9929f8f5 391 if((d0[0]*d0[1]) > fd0d0){continue;}
392 if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
393
394 if(fplothisto){
3a626a12 395 //fD0mass->Fill(minv);
396 //fD0mass->Fill(minvbar);
397 fD0mass->Fill(fd0calc->InvMass(tN,tP));
398 fD0mass->Fill(fd0calc->InvMass(tP,tN));
399 /*
9929f8f5 400 if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
401 fD0mass->Fill(fd0calc->InvMass(tN,tP));
402 }
403 else{
404 fD0mass->Fill(fd0calc->InvMass(tP,tN));
3a626a12 405 }
406 */
9929f8f5 407 }
3a626a12 408
9929f8f5 409 nD0++;
410 delete vertexp1n1;
411 }
412 }
413 return nD0;
414}
415
416Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
417 int nD0=0;
418 Double_t d0[2];
419 Double_t D0,D0bar;
420 Double_t svpos[3];
421 Double_t pvpos[3];
422
423 AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
424 event->GetStdContent();
425 Int_t nV0 = event->GetNumberOfV0s();
426 Double_t field = event->GetMagneticField();
427 const AliESDVertex* pv = event->GetPrimaryVertexTracks();
428 pv->GetXYZ(pvpos);
429
430 for (Int_t iv=0; iv<nV0; iv++) {
431
432 AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
433 AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());
434
435 if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
436
437 d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
438 d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
439
440 if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
441
442 event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
443
444 if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
445 HLTInfo("Tracks could not be propagated to secondary vertex");
446 continue;
447 }
448
449 Double_t tmp1, tmp2;
450 if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
451
452 if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - mD0PDG) > finvMass){continue;}
453 fd0calc->cosThetaStar(tN,tP,D0,D0bar);
454 if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
455 if((d0[0]*d0[1]) > fd0d0){continue;}
456 if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
457
458 nD0++;
459 if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
460 fD0mass->Fill(fd0calc->InvMass(tN,tP));
461 }
462 else{
463 fD0mass->Fill(fd0calc->InvMass(tP,tN));
464 }
465 }
466 return nD0;
467
468}