8c91efa21f124a6139c93f5ec4f6cff462526d73
[u/mrichter/AliRoot.git] / HLT / trigger / AliHLTD0Trigger.cxx
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"
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"
41 #include "AliAODVertex.h"
42 #include "AliESDVertex.h"
43 #include "AliAODRecoDecay.h"
44 #include "TMap.h"
45 #include "AliHLTD0Candidate.h"
46 #include "TFile.h"
47 #include "TClonesArray.h"
48
49 /** ROOT macro for the implementation of ROOT specific class methods */
50 ClassImp(AliHLTD0Trigger)
51
52 AliHLTD0Trigger::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)
61   , fplothisto(false)
62   , fUseV0(false)
63   , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
64   , fD0mass(NULL)
65   , fD0pt(NULL)
66   , fPos()
67   , fNeg()
68   , fd0calc(NULL)                  
69   , ftwoTrackArray(NULL)
70   , fTotalD0(0)
71   , fuseKF(false)
72   , fSendCandidate(false)
73   , fCandidateTree(NULL)
74   , fCandidateArray(NULL)
75   , fWriteFile(false)
76 {
77   
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
86 const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
87
88 AliHLTD0Trigger::~AliHLTD0Trigger()
89 {
90   // destructor 
91 }
92   
93 const char* AliHLTD0Trigger::GetTriggerName() const
94 {
95   // see header file for class documentation
96   return "D0Trigger";
97 }
98
99 AliHLTComponent* AliHLTD0Trigger::Spawn()
100 {
101   // see header file for class documentation
102   return new AliHLTD0Trigger;
103 }
104
105 int AliHLTD0Trigger::DoTrigger()
106 {
107   // -- Iterator over Data Blocks --
108   //const AliHLTComponentBlockData* iter = NULL;
109
110   if (!IsDataEvent()) return 0;
111
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
114   Int_t nD0=0;
115   TString description;
116
117   const AliESDVertex *Vertex=NULL;
118   if(fCandidateArray){fCandidateArray->Clear();}
119   
120   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {   
121     if(fUseV0){
122       nD0=RecV0(iter);
123     }
124     else{
125        AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
126        event->GetStdContent();
127        Double_t field = event->GetMagneticField();
128        Vertex = event->GetPrimaryVertexTracks();
129        if(!Vertex || Vertex->GetNContributors()<2){
130          HLTDebug("Contributors in ESD vertex to low or not been set");
131          continue;
132        }
133        for(Int_t it=0;it<event->GetNumberOfTracks();it++){
134          SingleTrackSelect(event->GetTrack(it),Vertex,field);
135        }
136        RecD0(nD0,Vertex,field);       
137     }
138   }
139
140   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut); 
141         iter != NULL; iter = GetNextInputObject() ) { 
142     Vertex = dynamic_cast<const AliESDVertex*>(iter);
143     if(!Vertex){
144       HLTError("Vertex object is corrupted");
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     
154     Double_t field = GetBz();
155     if (Vertex){
156       for(UInt_t i=0;i<tracksVector.size();i++){
157         SingleTrackSelect(&tracksVector[i],Vertex,field);
158       }
159       RecD0(nD0,Vertex,field);
160     }    
161   }
162   
163   fTotalD0+=nD0;
164
165   ftwoTrackArray->Clear();
166   fPos.clear();
167   fNeg.clear();
168      
169   HLTDebug("Number of D0 found: %d",nD0);
170   HLTDebug("Total Number of D0 found: %d",fTotalD0);
171  
172   if(fplothisto){
173     PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
174     PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
175   }
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
186   //if (iResult>=0) {
187   if (1) {
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);
214   //return iResult;
215   return 0;
216 }
217
218 int AliHLTD0Trigger::DoInit(int argc, const char** argv)
219 {
220   // component initialization
221
222   fd0calc = new AliHLTD0toKpi();
223   ftwoTrackArray = new TObjArray(2);
224   // see header file for class documentation
225   fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
226   fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
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);
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
245   return iResult;
246 }
247
248 int AliHLTD0Trigger::DoDeinit()
249 {
250   // see header file for class documentation
251   if(fd0calc){delete fd0calc; fd0calc = NULL;}  
252   if(fD0mass){delete fD0mass; fD0mass = NULL;}
253   if(fD0pt){delete fD0pt; fD0pt=NULL;}
254   if(ftwoTrackArray){delete ftwoTrackArray; ftwoTrackArray=NULL;}
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;}
262   return 0;
263 }
264
265 int 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
278 int 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) {
336     fplothisto=true;
337     return 1;
338   }
339   if (argument.CompareTo("-usev0")==0) {
340     fUseV0=true;
341     return 1;
342   }
343   if (argument.CompareTo("-useKF")==0) {
344     fuseKF=true;
345     return 1;
346   }
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  
356   // unknown argument
357   return -EINVAL;
358 }
359
360 void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, const AliESDVertex* pv, Double_t field){
361   // Offline uses || on the cuts of decay products
362   if (!pv) return;
363
364   Double_t pvpos[3];
365   pv->GetXYZ(pvpos);
366
367   if(t->Pt()<fPtMin){return;}
368   if(TMath::Abs(t->GetD(pvpos[0],pvpos[1],field)) > fd0){return;}
369
370   if(t->Charge()>0){
371     fPos.push_back(t);
372   }
373   else{
374     fNeg.push_back(t);
375   }
376 }
377
378 void AliHLTD0Trigger::RecD0(Int_t& nD0,const AliESDVertex* pv, Double_t field){
379   // Default way of reconstructing D0's. Both from ESD and Barreltracks
380   Double_t D0,D0bar,xdummy,ydummy; 
381   Double_t d0[2];
382   Double_t svpos[3];
383   Double_t pvpos[3];
384   
385   if(!pv){
386     HLTError("No Vertex is set");
387     return;
388   }
389   pv->GetXYZ(pvpos);
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           
396       tP->PropagateToDCA(pv,field,kVeryBig);  //do I need this??????
397       tN->PropagateToDCA(pv,field,kVeryBig);
398       
399       Double_t dcatPtN = tP->GetDCA(tN,field,xdummy,ydummy);
400       if(dcatPtN>fdca) { continue; }
401       
402       ftwoTrackArray->AddAt(tP,0);
403       ftwoTrackArray->AddAt(tN,1);
404       AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,field,pv,fuseKF);
405       if(!vertexp1n1) { 
406         ftwoTrackArray->Clear();
407         continue; 
408       }
409       
410       vertexp1n1->GetXYZ(svpos);
411       
412       tP->PropagateToDCA(vertexp1n1,field,kVeryBig); 
413       tN->PropagateToDCA(vertexp1n1,field,kVeryBig);
414       
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);
437       if(TMath::Abs(minv-fD0PDG)>finvMass && TMath::Abs(minv-fD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
438       */
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       
446       fd0calc->CosThetaStar(tN,tP,D0,D0bar);
447       if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){
448         vertexp1n1->RemoveCovMatrix();
449         delete vertexp1n1;
450         ftwoTrackArray->Clear();
451         continue;
452       }
453       
454       d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
455       d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
456       if((d0[0]*d0[1]) > fd0d0){
457         vertexp1n1->RemoveCovMatrix();
458         delete vertexp1n1;
459         ftwoTrackArray->Clear();
460         continue;
461       }
462       
463       if(fd0calc->PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){
464         vertexp1n1->RemoveCovMatrix();
465         delete vertexp1n1;
466         ftwoTrackArray->Clear();            
467         continue;
468       }
469       
470       if(fplothisto){
471         //fD0mass->Fill(minv);
472         //fD0mass->Fill(minvbar);
473         fD0mass->Fill(fd0calc->InvMass(tN,tP));
474         fD0mass->Fill(fd0calc->InvMass(tP,tN));
475         fD0pt->Fill(fd0calc->Pt(tP,tN));
476         /*
477         if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
478           fD0mass->Fill(fd0calc->InvMass(tN,tP));
479         }
480         else{
481           fD0mass->Fill(fd0calc->InvMass(tP,tN));
482           }
483         */
484       }
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());
488       }
489
490       nD0++;
491       vertexp1n1->RemoveCovMatrix();
492       delete vertexp1n1;
493       ftwoTrackArray->Clear();
494     }
495   }
496   ftwoTrackArray->Clear();
497 }
498
499 Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
500   // Reconstruction D0's from the V0 found by the V0 finder
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 ) );
508   if (!event) {
509     HLTError("input object is not of type AliESDEvent");
510     return 0;
511   }
512
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       
541       if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - fD0PDG) > finvMass){continue;}
542       fd0calc->CosThetaStar(tN,tP,D0,D0bar);
543       if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
544       if((d0[0]*d0[1]) > fd0d0){continue;}
545       if(fd0calc->PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
546       
547       nD0++;
548       if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
549         fD0mass->Fill(fd0calc->InvMass(tN,tP));
550       }
551       else{
552         fD0mass->Fill(fd0calc->InvMass(tP,tN));
553       }
554   }
555   return nD0;
556   
557 }
558 void 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 }