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