Updated macros for PHOS alignment calculation
[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"
10d746a6 44#include "TMap.h"
0f4240c3 45#include "AliHLTD0Candidate.h"
46#include "TFile.h"
47#include "TClonesArray.h"
5f4502cc 48
49/** ROOT macro for the implementation of ROOT specific class methods */
50ClassImp(AliHLTD0Trigger)
51
52AliHLTD0Trigger::AliHLTD0Trigger()
53 : AliHLTTrigger()
54 , fPtMin(0.0)
55 , fdca(0.0)
56 , finvMass(0.0)
57 , fcosThetaStar(0.0)
58 , fd0(0.0)
59 , fd0d0(0.0)
60 , fcosPoint(0.0)
9929f8f5 61 , fplothisto(false)
62 , fUseV0(false)
052f676e 63 , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
5f4502cc 64 , fD0mass(NULL)
252bba71 65 , fD0pt(NULL)
9929f8f5 66 , fPos()
67 , fNeg()
68 , fd0calc(NULL)
69 , ftwoTrackArray(NULL)
70 , fTotalD0(0)
15f6cee9 71 , fuseKF(false)
0f4240c3 72 , fSendCandidate(false)
73 , fCandidateTree(NULL)
74 , fCandidateArray(NULL)
75 , fWriteFile(false)
5f4502cc 76{
9929f8f5 77
5f4502cc 78 // see header file for class documentation
79 // or
80 // refer to README to build package
81 // or
82 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
83
84}
85
86const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
87
88AliHLTD0Trigger::~AliHLTD0Trigger()
89{
0f4240c3 90 // destructor
5f4502cc 91}
9929f8f5 92
5f4502cc 93const char* AliHLTD0Trigger::GetTriggerName() const
94{
95 // see header file for class documentation
96 return "D0Trigger";
97}
98
99AliHLTComponent* AliHLTD0Trigger::Spawn()
100{
101 // see header file for class documentation
102 return new AliHLTD0Trigger;
103}
104
105int AliHLTD0Trigger::DoTrigger()
106{
9929f8f5 107 // -- Iterator over Data Blocks --
108 //const AliHLTComponentBlockData* iter = NULL;
5f4502cc 109
9929f8f5 110 if (!IsDataEvent()) return 0;
5f4502cc 111
9f96f527 112 HLTDebug("Used Cuts: pT: %f , DCA: %f , InvMass: %f , CosThetaStar: %f , d0: %f , d0xd0: %f , CosPointingAngle: %f",fPtMin,fdca,finvMass,fcosThetaStar,fd0,fd0d0,fcosPoint);
113
9929f8f5 114 Int_t nD0=0;
115 TString description;
5f4502cc 116
a7cf6e02 117 const AliESDVertex *Vertex=NULL;
0f4240c3 118 if(fCandidateArray){fCandidateArray->Clear();}
119
9929f8f5 120 for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
121 if(fUseV0){
122 nD0=RecV0(iter);
123 }
124 else{
629b904b 125 AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
126 event->GetStdContent();
14978adc 127 Double_t field = event->GetMagneticField();
a7cf6e02 128 Vertex = event->GetPrimaryVertexTracks();
129 if(!Vertex || Vertex->GetNContributors()<2){
9f96f527 130 HLTDebug("Contributors in ESD vertex to low or not been set");
3a626a12 131 continue;
132 }
629b904b 133 for(Int_t it=0;it<event->GetNumberOfTracks();it++){
14978adc 134 SingleTrackSelect(event->GetTrack(it),Vertex,field);
629b904b 135 }
0f4240c3 136 RecD0(nD0,Vertex,field);
5f4502cc 137 }
5f4502cc 138 }
9929f8f5 139
3a626a12 140 for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut);
629b904b 141 iter != NULL; iter = GetNextInputObject() ) {
a7cf6e02 142 Vertex = dynamic_cast<const AliESDVertex*>(iter);
143 if(!Vertex){
3a626a12 144 HLTError("Vertex object is corrupted");
629b904b 145 //iResult = -EINVAL;
146 }
147 }
148
149 for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
150 iter != NULL; iter = GetNextInputBlock() ) {
151 vector<AliHLTGlobalBarrelTrack> tracksVector;
152 AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
153
14978adc 154 Double_t field = GetBz();
a7cf6e02 155 if (Vertex){
156 for(UInt_t i=0;i<tracksVector.size();i++){
14978adc 157 SingleTrackSelect(&tracksVector[i],Vertex,field);
a7cf6e02 158 }
0f4240c3 159 RecD0(nD0,Vertex,field);
a7cf6e02 160 }
161 }
162
9929f8f5 163 fTotalD0+=nD0;
2a743e38 164
9929f8f5 165 ftwoTrackArray->Clear();
166 fPos.clear();
167 fNeg.clear();
629b904b 168
2a743e38 169 HLTDebug("Number of D0 found: %d",nD0);
2a743e38 170 HLTDebug("Total Number of D0 found: %d",fTotalD0);
0f4240c3 171
252bba71 172 if(fplothisto){
173 PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
174 PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
175 }
0f4240c3 176 if(fSendCandidate){
177 PushBack( (TObject*) fCandidateArray, kAliHLTDataTypeTObjArray,0);
178 }
179 if(fWriteFile && fCandidateTree!=NULL){
180 fCandidateTree->Fill();
181 }
182 if(fCandidateArray){
183 fCandidateArray->Clear();
184 }
185
9929f8f5 186 //if (iResult>=0) {
187 if (1) {
5f4502cc 188
189 if (nD0>=1) {
190 description.Form("Event contains %d D0(s)", nD0);
191 SetDescription(description.Data());
192 // Enable the central detectors for readout.
193 GetReadoutList().Enable(
194 AliHLTReadoutList::kITSSPD |
195 AliHLTReadoutList::kITSSDD |
196 AliHLTReadoutList::kITSSSD |
197 AliHLTReadoutList::kTPC |
198 AliHLTReadoutList::kTRD |
199 AliHLTReadoutList::kTOF |
200 AliHLTReadoutList::kHMPID |
201 AliHLTReadoutList::kPHOS
202 );
203 // Add the available HLT information for readout too.
204 GetTriggerDomain().Add("CLUSTERS", "TPC ");
205 TriggerEvent(true);
206 return 0;
207 }
208 description.Form("No D0");
209 } else {
210 description.Form("No input blocks found");
211 }
212 SetDescription(description.Data());
213 TriggerEvent(false);
9929f8f5 214 //return iResult;
215 return 0;
5f4502cc 216}
217
218int AliHLTD0Trigger::DoInit(int argc, const char** argv)
219{
0f4240c3 220 // component initialization
221
9929f8f5 222 fd0calc = new AliHLTD0toKpi();
223 ftwoTrackArray = new TObjArray(2);
5f4502cc 224 // see header file for class documentation
9929f8f5 225 fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
252bba71 226 fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
5f4502cc 227 // first configure the default
228 int iResult=0;
229 if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
230
231 // configure from the command line parameters if specified
232 if (iResult>=0 && argc>0)
233 iResult=ConfigureFromArgumentString(argc, argv);
0f4240c3 234
235 if (iResult<0) return iResult;
236
237 if(fSendCandidate || fWriteFile){
238 fCandidateArray = new TClonesArray("AliHLTD0Candidate",1000);
239 }
240 if(fWriteFile){
241 fCandidateTree = new TTree("D0Candidates","Tree with D0 Candidates");
242 fCandidateTree->Branch("Candidates","TClonesArray",&fCandidateArray,32000,99);
243 }
244
5f4502cc 245 return iResult;
246}
247
248int AliHLTD0Trigger::DoDeinit()
249{
250 // see header file for class documentation
a7cf6e02 251 if(fd0calc){delete fd0calc; fd0calc = NULL;}
252 if(fD0mass){delete fD0mass; fD0mass = NULL;}
925bbd08 253 if(fD0pt){delete fD0pt; fD0pt=NULL;}
254 if(ftwoTrackArray){delete ftwoTrackArray; ftwoTrackArray=NULL;}
0f4240c3 255 if(fWriteFile){
256 TFile *mcFile = new TFile("D0_Candidates.root","RECREATE");
257 fCandidateTree->Write();
258 mcFile->Close();
259 }
260 if(fCandidateTree){delete fCandidateTree; fCandidateTree=NULL;}
261 if(fCandidateArray){delete fCandidateArray; fCandidateArray=NULL;}
5f4502cc 262 return 0;
263}
264
265int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
266{
267 // see header file for class documentation
268
269 // configure from the specified antry or the default one
270 const char* entry=cdbEntry;
271 if (!entry || entry[0]==0) {
272 entry=fgkOCDBEntry;
273 }
274
275 return ConfigureFromCDBTObjString(entry);
276}
277
278int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
279{
280 // see header file for class documentation
281 if (argc<=0) return 0;
282 int i=0;
283 TString argument=argv[i];
284 // -minpt for decay
285 if (argument.CompareTo("-pt")==0) {
286 if (++i>=argc) return -EPROTO;
287 argument=argv[i];
288 fPtMin=argument.Atof();
289 return 2;
290 }
291 // minimum dca for decay tracks
292 if (argument.CompareTo("-dca")==0) {
293 if (++i>=argc) return -EPROTO;
294 argument=argv[i];
295 fdca=argument.Atof();
296 return 2;
297 }
298 // inv. mass half width.
299 if (argument.CompareTo("-invmass")==0) {
300 if (++i>=argc) return -EPROTO;
301 argument=argv[i];
302 finvMass=argument.Atof();
303 return 2;
304 }
305
306// cos theta for decay angle
307 if (argument.CompareTo("-costhetastar")==0) {
308 if (++i>=argc) return -EPROTO;
309 argument=argv[i];
310 fcosThetaStar=argument.Atof();
311 return 2;
312 }
313
314 // impact parameter for decay
315 if (argument.CompareTo("-d0")==0) {
316 if (++i>=argc) return -EPROTO;
317 argument=argv[i];
318 fd0=argument.Atof();
319 return 2;
320 }
321 // product of impact parameter
322 if (argument.CompareTo("-d0d0")==0) {
323 if (++i>=argc) return -EPROTO;
324 argument=argv[i];
325 fd0d0=argument.Atof();
326 return 2;
327 }
328 // product of impact parameter
329 if (argument.CompareTo("-cospoint")==0) {
330 if (++i>=argc) return -EPROTO;
331 argument=argv[i];
332 fcosPoint=argument.Atof();
333 return 2;
334 }
335 if (argument.CompareTo("-plothistogram")==0) {
5f4502cc 336 fplothisto=true;
9929f8f5 337 return 1;
338 }
339 if (argument.CompareTo("-usev0")==0) {
340 fUseV0=true;
341 return 1;
5f4502cc 342 }
15f6cee9 343 if (argument.CompareTo("-useKF")==0) {
344 fuseKF=true;
345 return 1;
346 }
0f4240c3 347 if (argument.CompareTo("-send-candidates")==0) {
348 fSendCandidate=true;
349 return 1;
350 }
351 if (argument.CompareTo("-write-file")==0) {
352 fWriteFile=true;
353 return 1;
354 }
355
5f4502cc 356 // unknown argument
357 return -EINVAL;
358}
9929f8f5 359
0f4240c3 360void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, const AliESDVertex* pv, Double_t field){
a7cf6e02 361 // Offline uses || on the cuts of decay products
362 if (!pv) return;
8d9d674e 363
a7cf6e02 364 Double_t pvpos[3];
365 pv->GetXYZ(pvpos);
629b904b 366
9929f8f5 367 if(t->Pt()<fPtMin){return;}
14978adc 368 if(TMath::Abs(t->GetD(pvpos[0],pvpos[1],field)) > fd0){return;}
9929f8f5 369
370 if(t->Charge()>0){
371 fPos.push_back(t);
372 }
373 else{
374 fNeg.push_back(t);
375 }
376}
377
0f4240c3 378void AliHLTD0Trigger::RecD0(Int_t& nD0,const AliESDVertex* pv, Double_t field){
a7cf6e02 379 // Default way of reconstructing D0's. Both from ESD and Barreltracks
9929f8f5 380 Double_t D0,D0bar,xdummy,ydummy;
381 Double_t d0[2];
382 Double_t svpos[3];
383 Double_t pvpos[3];
384
a7cf6e02 385 if(!pv){
629b904b 386 HLTError("No Vertex is set");
2a743e38 387 return;
9929f8f5 388 }
a7cf6e02 389 pv->GetXYZ(pvpos);
629b904b 390
391 for(UInt_t ip=0;ip<fPos.size();ip++){
392 AliExternalTrackParam *tP=fPos[ip];
393 for(UInt_t in=0;in<fNeg.size();in++){
394 AliExternalTrackParam *tN=fNeg[in];
395
14978adc 396 tP->PropagateToDCA(pv,field,kVeryBig); //do I need this??????
397 tN->PropagateToDCA(pv,field,kVeryBig);
9929f8f5 398
14978adc 399 Double_t dcatPtN = tP->GetDCA(tN,field,xdummy,ydummy);
9929f8f5 400 if(dcatPtN>fdca) { continue; }
401
629b904b 402 ftwoTrackArray->AddAt(tP,0);
403 ftwoTrackArray->AddAt(tN,1);
14978adc 404 AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,field,pv,fuseKF);
9929f8f5 405 if(!vertexp1n1) {
406 ftwoTrackArray->Clear();
407 continue;
408 }
409
410 vertexp1n1->GetXYZ(svpos);
411
14978adc 412 tP->PropagateToDCA(vertexp1n1,field,kVeryBig);
413 tN->PropagateToDCA(vertexp1n1,field,kVeryBig);
9929f8f5 414
3a626a12 415 /*
416 Double_t px[2],py[2],pz[2];
417 Double_t momentum[3];
418 tP->GetPxPyPz(momentum);
419 px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
420 tN->GetPxPyPz(momentum);
421 px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
422
423 Short_t dummycharge=0;
424 Double_t *dummyd0 = new Double_t[2];
425 Int_t nprongs = 2;
426
427 for(Int_t ipr=0;ipr<nprongs;ipr++) dummyd0[ipr]=0.;
428 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
429 delete [] dummyd0; dummyd0=NULL;
430
431 UInt_t pdg2[2],pdg2bar[2];
432 Double_t mPDG,minv,minvbar;
433
434 pdg2[0]=211; pdg2[1]=321; pdg2bar[0]=321; pdg2bar[1]=211;
435 minv = rd->InvMass(nprongs,pdg2);
436 minvbar = rd->InvMass(nprongs,pdg2bar);
052f676e 437 if(TMath::Abs(minv-fD0PDG)>finvMass && TMath::Abs(minv-fD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
3a626a12 438 */
0f4240c3 439 if((TMath::Abs(fd0calc->InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-fD0PDG) > finvMass){
440 vertexp1n1->RemoveCovMatrix();
441 delete vertexp1n1;
442 ftwoTrackArray->Clear();
443 continue;
444 }
445
36e4e16c 446 fd0calc->CosThetaStar(tN,tP,D0,D0bar);
0f4240c3 447 if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){
448 vertexp1n1->RemoveCovMatrix();
449 delete vertexp1n1;
450 ftwoTrackArray->Clear();
451 continue;
452 }
453
14978adc 454 d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
455 d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
0f4240c3 456 if((d0[0]*d0[1]) > fd0d0){
457 vertexp1n1->RemoveCovMatrix();
458 delete vertexp1n1;
459 ftwoTrackArray->Clear();
460 continue;
461 }
462
36e4e16c 463 if(fd0calc->PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){
0f4240c3 464 vertexp1n1->RemoveCovMatrix();
465 delete vertexp1n1;
466 ftwoTrackArray->Clear();
467 continue;
468 }
9929f8f5 469
470 if(fplothisto){
3a626a12 471 //fD0mass->Fill(minv);
472 //fD0mass->Fill(minvbar);
473 fD0mass->Fill(fd0calc->InvMass(tN,tP));
474 fD0mass->Fill(fd0calc->InvMass(tP,tN));
252bba71 475 fD0pt->Fill(fd0calc->Pt(tP,tN));
3a626a12 476 /*
052f676e 477 if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
9929f8f5 478 fD0mass->Fill(fd0calc->InvMass(tN,tP));
479 }
480 else{
481 fD0mass->Fill(fd0calc->InvMass(tP,tN));
3a626a12 482 }
483 */
9929f8f5 484 }
0f4240c3 485
486 if(fSendCandidate && fCandidateArray){
487 new(fCandidateArray->AddrAt(nD0)) AliHLTD0Candidate(fd0calc->InvMass(tN,tP),fd0calc->InvMass(tP,tN),fd0calc->Pt(tP,tN),tP->GetLabel(),tN->GetLabel(),tP->Pt(),tN->Pt());
9f96f527 488 }
2a743e38 489
9929f8f5 490 nD0++;
052f676e 491 vertexp1n1->RemoveCovMatrix();
9929f8f5 492 delete vertexp1n1;
0f4240c3 493 ftwoTrackArray->Clear();
9929f8f5 494 }
495 }
0f4240c3 496 ftwoTrackArray->Clear();
9929f8f5 497}
498
499Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
a7cf6e02 500 // Reconstruction D0's from the V0 found by the V0 finder
9929f8f5 501 int nD0=0;
502 Double_t d0[2];
503 Double_t D0,D0bar;
504 Double_t svpos[3];
505 Double_t pvpos[3];
506
507 AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
55f4e3b5 508 if (!event) {
509 HLTError("input object is not of type AliESDEvent");
510 return 0;
511 }
512
9929f8f5 513 event->GetStdContent();
514 Int_t nV0 = event->GetNumberOfV0s();
515 Double_t field = event->GetMagneticField();
516 const AliESDVertex* pv = event->GetPrimaryVertexTracks();
517 pv->GetXYZ(pvpos);
518
519 for (Int_t iv=0; iv<nV0; iv++) {
520
521 AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
522 AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());
523
524 if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
525
526 d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
527 d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
528
529 if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
530
531 event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
532
533 if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
534 HLTInfo("Tracks could not be propagated to secondary vertex");
535 continue;
536 }
537
538 Double_t tmp1, tmp2;
539 if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
540
052f676e 541 if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - fD0PDG) > finvMass){continue;}
36e4e16c 542 fd0calc->CosThetaStar(tN,tP,D0,D0bar);
9929f8f5 543 if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
544 if((d0[0]*d0[1]) > fd0d0){continue;}
36e4e16c 545 if(fd0calc->PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
9929f8f5 546
547 nD0++;
052f676e 548 if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
9929f8f5 549 fD0mass->Fill(fd0calc->InvMass(tN,tP));
550 }
551 else{
552 fD0mass->Fill(fd0calc->InvMass(tP,tN));
553 }
554 }
555 return nD0;
556
557}
10d746a6 558void AliHLTD0Trigger::GetOCDBObjectDescription( TMap* const targetMap)
559{
560 // Get a list of OCDB object description.
561 if (!targetMap) return;
562 targetMap->Add(new TObjString("HLT/ConfigHLT/D0Trigger"),new TObjString("Object with default cuts for D0 reconstruction" ));
563}