correcting compilation error
[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 "AliHLTMCEvent.h"
45 #include "AliStack.h"
46 #include "TParticle.h"
47 #include "TMap.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   , fTotalD0Onetrue(0)
72   , fTotalD0true(0)
73   , fField(0)
74   , fEvent(NULL)
75   , fuseKF(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  
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   Int_t nD0Onetrue=0;
116   Int_t nD0true=0;
117   TString description;
118
119   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeMCObject|kAliHLTDataOriginHLT); iter != NULL; iter = GetNextInputObject() ) {   
120     fEvent = dynamic_cast<AliHLTMCEvent*>(const_cast<TObject*>( iter ) );
121     if ( ! fEvent ) {
122       HLTError( "No MC Event present!" );
123       break;
124     }
125   }
126
127   const AliESDVertex *Vertex=NULL;
128
129   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {   
130     if(fUseV0){
131       nD0=RecV0(iter);
132     }
133     else{
134        AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
135        event->GetStdContent();
136        fField = event->GetMagneticField();
137        Vertex = event->GetPrimaryVertexTracks();
138        if(!Vertex || Vertex->GetNContributors()<2){
139          HLTDebug("Contributors in ESD vertex to low or not been set");
140          continue;
141        }
142        for(Int_t it=0;it<event->GetNumberOfTracks();it++){
143          SingleTrackSelect(event->GetTrack(it),Vertex);
144        }
145        RecD0(nD0,nD0Onetrue,nD0true,Vertex);       
146     }
147   }
148
149   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut); 
150         iter != NULL; iter = GetNextInputObject() ) { 
151     Vertex = dynamic_cast<const AliESDVertex*>(iter);
152     if(!Vertex){
153       HLTError("Vertex object is corrupted");
154       //iResult = -EINVAL;    
155     }
156   }  
157    
158   for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS); 
159         iter != NULL; iter = GetNextInputBlock() ) { 
160     vector<AliHLTGlobalBarrelTrack> tracksVector;
161     AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
162     
163     fField = GetBz();
164     if (Vertex){
165       for(UInt_t i=0;i<tracksVector.size();i++){
166         SingleTrackSelect(&tracksVector[i],Vertex);
167       }
168       RecD0(nD0,nD0Onetrue,nD0true,Vertex);
169     }    
170   }
171   
172   fTotalD0+=nD0;
173   fTotalD0Onetrue += nD0Onetrue;
174   fTotalD0true += nD0true;
175
176   ftwoTrackArray->Clear();
177   fPos.clear();
178   fNeg.clear();
179      
180   HLTDebug("Number of D0 found: %d",nD0);
181   HLTDebug("Number of D0 candidates with one track as D0 decay: %d",nD0Onetrue);
182   HLTDebug("Number of True D0 found: %d",nD0true);
183   HLTDebug("Total Number of D0 found: %d",fTotalD0);
184   HLTDebug("Number of D0 candidates with one track as D0 decay: %d",fTotalD0Onetrue);
185   HLTDebug("Total Number of True D0 found: %d",fTotalD0true);
186
187   if(fplothisto){
188     PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
189     PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
190   }
191   
192   //if (iResult>=0) {
193   if (1) {
194    
195     if (nD0>=1) {
196       description.Form("Event contains %d D0(s)", nD0);
197       SetDescription(description.Data());
198       // Enable the central detectors for readout.
199       GetReadoutList().Enable(
200                               AliHLTReadoutList::kITSSPD |
201                               AliHLTReadoutList::kITSSDD |
202                               AliHLTReadoutList::kITSSSD |
203                               AliHLTReadoutList::kTPC |
204                               AliHLTReadoutList::kTRD |
205                               AliHLTReadoutList::kTOF |
206                               AliHLTReadoutList::kHMPID |
207                               AliHLTReadoutList::kPHOS
208                               );
209       // Add the available HLT information for readout too.
210       GetTriggerDomain().Add("CLUSTERS", "TPC ");
211       TriggerEvent(true);
212       return 0;
213     }
214     description.Form("No D0");
215   } else {
216     description.Form("No input blocks found");
217   }
218   SetDescription(description.Data());
219   TriggerEvent(false);
220   //return iResult;
221   return 0;
222 }
223
224 int AliHLTD0Trigger::DoInit(int argc, const char** argv)
225 {
226   
227   fd0calc = new AliHLTD0toKpi();
228   ftwoTrackArray = new TObjArray(2);
229
230   fplothisto=false;
231   // see header file for class documentation
232   fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
233   fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
234   // first configure the default
235   int iResult=0;
236   if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
237
238   // configure from the command line parameters if specified
239   if (iResult>=0 && argc>0)
240     iResult=ConfigureFromArgumentString(argc, argv);
241   return iResult;
242 }
243
244 int AliHLTD0Trigger::DoDeinit()
245 {
246   // see header file for class documentation
247   if(fd0calc){delete fd0calc; fd0calc = NULL;}  
248   if(fD0mass){delete fD0mass; fD0mass = NULL;}
249   if(fD0pt){delete fD0pt; fD0pt=NULL;}
250   if(ftwoTrackArray){delete ftwoTrackArray; ftwoTrackArray=NULL;}
251   
252   return 0;
253 }
254
255 int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
256 {
257   // see header file for class documentation
258
259   // configure from the specified antry or the default one
260   const char* entry=cdbEntry;
261   if (!entry || entry[0]==0) {
262     entry=fgkOCDBEntry;
263   }
264
265   return ConfigureFromCDBTObjString(entry);
266 }
267
268 int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
269 {
270   // see header file for class documentation
271   if (argc<=0) return 0;
272   int i=0;
273   TString argument=argv[i];
274   // -minpt for decay
275   if (argument.CompareTo("-pt")==0) {
276     if (++i>=argc) return -EPROTO;
277     argument=argv[i];
278     fPtMin=argument.Atof();
279     return 2;
280   }    
281   // minimum dca for decay tracks
282   if (argument.CompareTo("-dca")==0) {
283     if (++i>=argc) return -EPROTO;
284     argument=argv[i];
285     fdca=argument.Atof();
286     return 2;
287   }
288   // inv. mass half width.
289   if (argument.CompareTo("-invmass")==0) {
290     if (++i>=argc) return -EPROTO;
291     argument=argv[i];
292     finvMass=argument.Atof();
293     return 2;
294   }
295
296 // cos theta for decay angle
297   if (argument.CompareTo("-costhetastar")==0) {
298     if (++i>=argc) return -EPROTO;
299     argument=argv[i];
300     fcosThetaStar=argument.Atof();
301     return 2;
302   }
303
304   // impact parameter for decay
305   if (argument.CompareTo("-d0")==0) {
306     if (++i>=argc) return -EPROTO;
307     argument=argv[i];
308     fd0=argument.Atof();
309     return 2;
310   }
311   // product of impact parameter
312   if (argument.CompareTo("-d0d0")==0) {
313     if (++i>=argc) return -EPROTO;
314     argument=argv[i];
315     fd0d0=argument.Atof();
316     return 2;
317   }
318   // product of impact parameter
319   if (argument.CompareTo("-cospoint")==0) {
320     if (++i>=argc) return -EPROTO;
321     argument=argv[i];
322     fcosPoint=argument.Atof();
323     return 2;
324   }
325   if (argument.CompareTo("-plothistogram")==0) {
326     fplothisto=true;
327     return 1;
328   }
329   if (argument.CompareTo("-usev0")==0) {
330     fUseV0=true;
331     return 1;
332   }
333   if (argument.CompareTo("-useKF")==0) {
334     fuseKF=true;
335     return 1;
336   }
337   // unknown argument
338   return -EINVAL;
339 }
340
341 void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, const AliESDVertex* pv){
342   // Offline uses || on the cuts of decay products
343   if (!pv) return;
344
345   Double_t pvpos[3];
346   pv->GetXYZ(pvpos);
347
348   if(t->Pt()<fPtMin){return;}
349   if(TMath::Abs(t->GetD(pvpos[0],pvpos[1],fField)) > fd0){return;}
350
351   if(t->Charge()>0){
352     fPos.push_back(t);
353   }
354   else{
355     fNeg.push_back(t);
356   }
357 }
358
359 void AliHLTD0Trigger::RecD0(Int_t& nD0,Int_t& nD0Onetrue ,Int_t& nD0true,const AliESDVertex* pv){
360   // Default way of reconstructing D0's. Both from ESD and Barreltracks
361   Double_t D0,D0bar,xdummy,ydummy; 
362   Double_t d0[2];
363   Double_t svpos[3];
364   Double_t pvpos[3];
365   
366   if(!pv){
367     HLTError("No Vertex is set");
368     return;
369   }
370   pv->GetXYZ(pvpos);
371     
372   for(UInt_t ip=0;ip<fPos.size();ip++){
373     AliExternalTrackParam *tP=fPos[ip];
374     for(UInt_t in=0;in<fNeg.size();in++){
375       AliExternalTrackParam *tN=fNeg[in];
376           
377       tP->PropagateToDCA(pv,fField,kVeryBig);  //do I need this??????
378       tN->PropagateToDCA(pv,fField,kVeryBig);
379       
380       Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
381       if(dcatPtN>fdca) { continue; }
382       
383       ftwoTrackArray->AddAt(tP,0);
384       ftwoTrackArray->AddAt(tN,1);
385       AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,fField,pv,fuseKF);
386       if(!vertexp1n1) { 
387         ftwoTrackArray->Clear();
388         continue; 
389       }
390       
391       vertexp1n1->GetXYZ(svpos);
392       
393       tP->PropagateToDCA(vertexp1n1,fField,kVeryBig); 
394       tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
395       
396       /*
397       Double_t px[2],py[2],pz[2];
398       Double_t momentum[3];
399       tP->GetPxPyPz(momentum);
400       px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
401       tN->GetPxPyPz(momentum);
402       px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
403
404       Short_t dummycharge=0;
405       Double_t *dummyd0 = new Double_t[2];
406       Int_t nprongs = 2;
407
408       for(Int_t ipr=0;ipr<nprongs;ipr++) dummyd0[ipr]=0.;
409       AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
410       delete [] dummyd0; dummyd0=NULL;
411
412       UInt_t pdg2[2],pdg2bar[2];
413       Double_t mPDG,minv,minvbar;
414
415       pdg2[0]=211; pdg2[1]=321; pdg2bar[0]=321; pdg2bar[1]=211;
416       minv = rd->InvMass(nprongs,pdg2);
417       minvbar = rd->InvMass(nprongs,pdg2bar);
418       if(TMath::Abs(minv-fD0PDG)>finvMass && TMath::Abs(minv-fD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
419       */
420       if((TMath::Abs(fd0calc->InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-fD0PDG) > finvMass){continue;}
421       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
422       if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){continue;}
423       d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
424       d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
425       if((d0[0]*d0[1]) > fd0d0){continue;}
426       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
427       
428       if(fplothisto){
429         //fD0mass->Fill(minv);
430         //fD0mass->Fill(minvbar);
431         fD0mass->Fill(fd0calc->InvMass(tN,tP));
432         fD0mass->Fill(fd0calc->InvMass(tP,tN));
433         fD0pt->Fill(fd0calc->Pt(tP,tN));
434         /*
435         if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
436           fD0mass->Fill(fd0calc->InvMass(tN,tP));
437         }
438         else{
439           fD0mass->Fill(fd0calc->InvMass(tP,tN));
440           }
441         */
442       }
443       
444       if(CheckTrackMC(tP,tN)==2){
445         nD0true++;
446       }
447       if(CheckTrackMC(tP,tN)==1){
448         nD0Onetrue++;
449       }
450
451       nD0++;
452       vertexp1n1->RemoveCovMatrix();
453       delete vertexp1n1;
454     }
455   }
456 }
457
458 Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
459   // Reconstruction D0's from the V0 found by the V0 finder
460   int nD0=0;
461   Double_t d0[2];
462   Double_t D0,D0bar; 
463   Double_t svpos[3];
464   Double_t pvpos[3];
465   
466   AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
467   event->GetStdContent();
468   Int_t nV0 = event->GetNumberOfV0s();
469   Double_t field = event->GetMagneticField();
470   const AliESDVertex* pv = event->GetPrimaryVertexTracks();
471   pv->GetXYZ(pvpos);
472
473   for (Int_t iv=0; iv<nV0; iv++) {
474       
475       AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
476       AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());      
477       
478       if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
479       
480       d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
481       d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
482       
483       if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
484       
485       event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
486       
487       if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
488         HLTInfo("Tracks could not be propagated to secondary vertex");
489         continue;
490       }
491       
492       Double_t tmp1, tmp2;
493       if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
494       
495       if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - fD0PDG) > finvMass){continue;}
496       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
497       if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
498       if((d0[0]*d0[1]) > fd0d0){continue;}
499       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
500       
501       nD0++;
502       if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
503         fD0mass->Fill(fd0calc->InvMass(tN,tP));
504       }
505       else{
506         fD0mass->Fill(fd0calc->InvMass(tP,tN));
507       }
508   }
509   return nD0;
510   
511 }
512 int AliHLTD0Trigger::CheckTrackMC(AliExternalTrackParam* pt, AliExternalTrackParam* pn){
513   // Checks if decay products both came from a D0 or one.
514   if(!fEvent){return 0;}
515
516   int lP = pt->GetLabel();
517   int lN = pn->GetLabel();
518
519   if(lN>=0 && lP>=0){
520     
521     int imP = (fEvent->Particle(lP))->GetFirstMother();
522     int imN = (fEvent->Particle(lN))->GetFirstMother();
523     
524     if(imP>=0 && imN>=0){
525       TParticle * mP = fEvent->Particle(imP);
526       TParticle * mN = fEvent->Particle(imN);
527       if((fabs(mP->GetPdgCode())==421 && fabs(mN->GetPdgCode())==421) && imP == imN){
528         return 2;
529       }
530       if((fabs(mP->GetPdgCode())==421 || fabs(mN->GetPdgCode())==421) && imP == imN){
531         return 1;
532       }
533     }
534   }
535   return 0;
536 }
537 void AliHLTD0Trigger::GetOCDBObjectDescription( TMap* const targetMap)
538 {
539   // Get a list of OCDB object description.
540   if (!targetMap) return;
541   targetMap->Add(new TObjString("HLT/ConfigHLT/D0Trigger"),new TObjString("Object with default cuts for D0 reconstruction" ));
542 }