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   , fTotalD0Onetrue(0)
  72   , fTotalD0true(0)
  73   , fVertex(NULL)
  74   , fField(0)
  75   , fEvent(NULL)
  76   , fuseKF(false)
  77 {
  78
  79   // see header file for class documentation
  80   // or
  81   // refer to README to build package
  82   // or
  83   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  84
  85 }
  86
  87 const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
  88
  89 AliHLTD0Trigger::~AliHLTD0Trigger()
  90 {
  91   // TODO: code audit 2010-07-23 delete obsolete code
  92   //if(fd0calc){delete fd0calc;}
  93   //if(fD0mass){delete fD0mass;}
  94   //if(ftwoTrackArray){delete ftwoTrackArray;}
  95   // see header file for class documentation
  96 }
  97
  98 const char* AliHLTD0Trigger::GetTriggerName() const
  99 {
 100   // see header file for class documentation
 101   return "D0Trigger";
 102 }
 103
 104 AliHLTComponent* AliHLTD0Trigger::Spawn()
 105 {
 106   // see header file for class documentation
 107   return new AliHLTD0Trigger;
 108 }
 109
 110 int AliHLTD0Trigger::DoTrigger()
 111 {
 112   // -- Iterator over Data Blocks --
 113   //const AliHLTComponentBlockData* iter = NULL;
 114
 115   if (!IsDataEvent()) return 0;
 116
 117   HLTDebug("Used  Cuts: pT: %f , DCA: %f , InvMass: %f , CosThetaStar: %f , d0: %f , d0xd0: %f , CosPointingAngle: %f",fPtMin,fdca,finvMass,fcosThetaStar,fd0,fd0d0,fcosPoint);
 118
 119   Int_t nD0=0;
 120   Int_t nD0Onetrue=0;
 121   Int_t nD0true=0;
 122   TString description;
 123
 124   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeMCObject|kAliHLTDataOriginHLT); iter != NULL; iter = GetNextInputObject() ) {
 125     fEvent = dynamic_cast<AliHLTMCEvent*>(const_cast<TObject*>( iter ) );
 126     if ( ! fEvent ) {
 127       HLTError( "No MC Event present!" );
 128       break;
 129     }
 130   }
 131
 132   // TODO: code audit 2010-07-23 get rid of fVertex as member variable, the object
 133   // is changing every event and already available after the object extraction
 134   // change internal functions to pass the vertex as parameter. Setting the variable
 135   // to NULL here to be coherent with checks further down.
 136   fVertex=NULL;
 137
 138   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
 139     if(fUseV0){
 140       nD0=RecV0(iter);
 141     }
 142     else{
 143        AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 144        event->GetStdContent();
 145        fField = event->GetMagneticField();
 146        const AliESDVertex* pv = event->GetPrimaryVertexTracks();
 147        // TODO: code audit 2010-07-23 this is a memory leak, fVertex is not deleted
 148        //fVertex =  new AliESDVertex(*pv);
 149        // use the vertex object from the ESD, right now we need to cast away the
 150        // const'ness, think about changing the functions to use const AliESDVertex* objects
 151        fVertex =  const_cast<AliESDVertex*>(pv);
 152        if(fVertex->GetNContributors()<2){
 153          HLTDebug("Contributors in ESD vertex to low or not been set");
 154          continue;
 155        }
 156        for(Int_t it=0;it<event->GetNumberOfTracks();it++){
 157          // TODO: code audit 2010-07-23 pass vertex object as parameter
 158          SingleTrackSelect(event->GetTrack(it));
 159        }
 160
 161        // TODO: code audit 2010-07-23 pass vertex object as parameter
 162        if (fVertex) RecD0(nD0,nD0Onetrue,nD0true);
 163     }
 164   }
 165
 166   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut);
 167         iter != NULL; iter = GetNextInputObject() ) {
 168     fVertex = dynamic_cast<AliESDVertex*>(const_cast<TObject*>( iter ));
 169     if(!fVertex){
 170       HLTError("Vertex object is corrupted");
 171       //iResult = -EINVAL;
 172     }
 173   }
 174
 175   for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
 176         iter != NULL; iter = GetNextInputBlock() ) {
 177     vector<AliHLTGlobalBarrelTrack> tracksVector;
 178     AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
 179
 180     fField = GetBz();
 181     for(UInt_t i=0;i<tracksVector.size();i++){
 182       SingleTrackSelect(&tracksVector[i]);
 183     }
 184     // TODO: code audit 2010-07-23 pass vertex object as parameter
 185     if (fVertex) RecD0(nD0,nD0Onetrue,nD0true);
 186   }
 187
 188   fTotalD0+=nD0;
 189   fTotalD0Onetrue += nD0Onetrue;
 190   fTotalD0true += nD0true;
 191
 192   ftwoTrackArray->Clear();
 193   fPos.clear();
 194   fNeg.clear();
 195
 196   HLTDebug("Number of D0 found: %d",nD0);
 197   HLTDebug("Number of D0 candidates with one track as D0 decay: %d",nD0Onetrue);
 198   HLTDebug("Number of True D0 found: %d",nD0true);
 199   HLTDebug("Total Number of D0 found: %d",fTotalD0);
 200   HLTDebug("Number of D0 candidates with one track as D0 decay: %d",fTotalD0Onetrue);
 201   HLTDebug("Total Number of True D0 found: %d",fTotalD0true);
 202
 203   if(fplothisto){
 204     PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
 205     PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
 206   }
 207
 208   //if (iResult>=0) {
 209   if (1) {
 210
 211     if (nD0>=1) {
 212       description.Form("Event contains %d D0(s)", nD0);
 213       SetDescription(description.Data());
 214       // Enable the central detectors for readout.
 215       GetReadoutList().Enable(
 216                               AliHLTReadoutList::kITSSPD |
 217                               AliHLTReadoutList::kITSSDD |
 218                               AliHLTReadoutList::kITSSSD |
 219                               AliHLTReadoutList::kTPC |
 220                               AliHLTReadoutList::kTRD |
 221                               AliHLTReadoutList::kTOF |
 222                               AliHLTReadoutList::kHMPID |
 223                               AliHLTReadoutList::kPHOS
 224                               );
 225       // Add the available HLT information for readout too.
 226       GetTriggerDomain().Add("CLUSTERS", "TPC ");
 227       TriggerEvent(true);
 228       return 0;
 229     }
 230     description.Form("No D0");
 231   } else {
 232     description.Form("No input blocks found");
 233   }
 234   SetDescription(description.Data());
 235   TriggerEvent(false);
 236   //return iResult;
 237   return 0;
 238 }
 239
 240 int AliHLTD0Trigger::DoInit(int argc, const char** argv)
 241 {
 242
 243   fd0calc = new AliHLTD0toKpi();
 244   ftwoTrackArray = new TObjArray(2);
 245
 246   fplothisto=false;
 247   // see header file for class documentation
 248   fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 249   fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
 250   // first configure the default
 251   int iResult=0;
 252   if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
 253
 254   // configure from the command line parameters if specified
 255   if (iResult>=0 && argc>0)
 256     iResult=ConfigureFromArgumentString(argc, argv);
 257   return iResult;
 258 }
 259
 260 int AliHLTD0Trigger::DoDeinit()
 261 {
 262   // see header file for class documentation
 263
 264   // TODO: code audit 2010-07-23 deleted variables to be set to NULL
 265   if(fd0calc){delete fd0calc;}
 266   if(fD0mass){delete fD0mass;}
 267   // TODO: code audit 2010-07-23 delete obsolete code
 268   //if(ftwoTrackArray){delete ftwoTrackArray;}
 269   //if(fVertex){delete fVertex;}
 270   return 0;
 271 }
 272
 273 int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
 274 {
 275   // see header file for class documentation
 276
 277   // configure from the specified antry or the default one
 278   const char* entry=cdbEntry;
 279   if (!entry || entry[0]==0) {
 280     entry=fgkOCDBEntry;
 281   }
 282
 283   return ConfigureFromCDBTObjString(entry);
 284 }
 285
 286 int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
 287 {
 288   // see header file for class documentation
 289   if (argc<=0) return 0;
 290   int i=0;
 291   TString argument=argv[i];
 292   // -minpt for decay
 293   if (argument.CompareTo("-pt")==0) {
 294     if (++i>=argc) return -EPROTO;
 295     argument=argv[i];
 296     fPtMin=argument.Atof();
 297     return 2;
 298   }
 299   // minimum dca for decay tracks
 300   if (argument.CompareTo("-dca")==0) {
 301     if (++i>=argc) return -EPROTO;
 302     argument=argv[i];
 303     fdca=argument.Atof();
 304     return 2;
 305   }
 306   // inv. mass half width.
 307   if (argument.CompareTo("-invmass")==0) {
 308     if (++i>=argc) return -EPROTO;
 309     argument=argv[i];
 310     finvMass=argument.Atof();
 311     return 2;
 312   }
 313
 314 // cos theta for decay angle
 315   if (argument.CompareTo("-costhetastar")==0) {
 316     if (++i>=argc) return -EPROTO;
 317     argument=argv[i];
 318     fcosThetaStar=argument.Atof();
 319     return 2;
 320   }
 321
 322   // impact parameter for decay
 323   if (argument.CompareTo("-d0")==0) {
 324     if (++i>=argc) return -EPROTO;
 325     argument=argv[i];
 326     fd0=argument.Atof();
 327     return 2;
 328   }
 329   // product of impact parameter
 330   if (argument.CompareTo("-d0d0")==0) {
 331     if (++i>=argc) return -EPROTO;
 332     argument=argv[i];
 333     fd0d0=argument.Atof();
 334     return 2;
 335   }
 336   // product of impact parameter
 337   if (argument.CompareTo("-cospoint")==0) {
 338     if (++i>=argc) return -EPROTO;
 339     argument=argv[i];
 340     fcosPoint=argument.Atof();
 341     return 2;
 342   }
 343   if (argument.CompareTo("-plothistogram")==0) {
 344     fplothisto=true;
 345     return 1;
 346   }
 347   if (argument.CompareTo("-usev0")==0) {
 348     fUseV0=true;
 349     return 1;
 350   }
 351   if (argument.CompareTo("-useKF")==0) {
 352     fuseKF=true;
 353     return 1;
 354   }
 355   // unknown argument
 356   return -EINVAL;
 357 }
 358
 359 void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t){
 360   // Offline har || på disse kuttene på de to henfallsproduktene
 361   // TODO: code audit 2010-07-23 translate all comments
 362   if (!fVertex) return;
 363
 364   Double_t pv[3];
 365   fVertex->GetXYZ(pv);
 366
 367   if(t->Pt()<fPtMin){return;}
 368   if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > 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,Int_t& nD0Onetrue ,Int_t& nD0true){
 379
 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(!fVertex){
 386     HLTError("No Vertex is set");
 387     return;
 388   }
 389   fVertex->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(fVertex,fField,kVeryBig);  //do I need this??????
 397       tN->PropagateToDCA(fVertex,fField,kVeryBig);
 398
 399       Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
 400       if(dcatPtN>fdca) { continue; }
 401
 402       ftwoTrackArray->AddAt(tP,0);
 403       ftwoTrackArray->AddAt(tN,1);
 404       AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,fField,fVertex,fuseKF);
 405       if(!vertexp1n1) {
 406         ftwoTrackArray->Clear();
 407         continue;
 408       }
 409
 410       vertexp1n1->GetXYZ(svpos);
 411
 412       tP->PropagateToDCA(vertexp1n1,fField,kVeryBig);
 413       tN->PropagateToDCA(vertexp1n1,fField,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-mD0PDG)>finvMass && TMath::Abs(minv-mD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
 438       */
 439       if((TMath::Abs(fd0calc->InvMass(tN,tP)-mD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-mD0PDG) > finvMass){continue;}
 440       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
 441       if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){continue;}
 442       d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
 443       d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
 444       if((d0[0]*d0[1]) > fd0d0){continue;}
 445       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
 446
 447       if(fplothisto){
 448         //fD0mass->Fill(minv);
 449         //fD0mass->Fill(minvbar);
 450         fD0mass->Fill(fd0calc->InvMass(tN,tP));
 451         fD0mass->Fill(fd0calc->InvMass(tP,tN));
 452         fD0pt->Fill(fd0calc->Pt(tP,tN));
 453         /*
 454         if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
 455           fD0mass->Fill(fd0calc->InvMass(tN,tP));
 456         }
 457         else{
 458           fD0mass->Fill(fd0calc->InvMass(tP,tN));
 459           }
 460         */
 461       }
 462
 463       if(CheckTrackMC(tP,tN)==2){
 464         nD0true++;
 465       }
 466       if(CheckTrackMC(tP,tN)==1){
 467         nD0Onetrue++;
 468       }
 469
 470       nD0++;
 471       delete vertexp1n1;
 472     }
 473   }
 474 }
 475
 476 Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
 477   int nD0=0;
 478   Double_t d0[2];
 479   Double_t D0,D0bar;
 480   Double_t svpos[3];
 481   Double_t pvpos[3];
 482
 483   AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 484   event->GetStdContent();
 485   Int_t nV0 = event->GetNumberOfV0s();
 486   Double_t field = event->GetMagneticField();
 487   const AliESDVertex* pv = event->GetPrimaryVertexTracks();
 488   pv->GetXYZ(pvpos);
 489
 490   for (Int_t iv=0; iv<nV0; iv++) {
 491
 492       AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
 493       AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());
 494
 495       if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
 496
 497       d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
 498       d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
 499
 500       if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
 501
 502       event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
 503
 504       if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
 505         HLTInfo("Tracks could not be propagated to secondary vertex");
 506         continue;
 507       }
 508
 509       Double_t tmp1, tmp2;
 510       if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
 511
 512       if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - mD0PDG) > finvMass){continue;}
 513       fd0calc->cosThetaStar(tN,tP,D0,D0bar);
 514       if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
 515       if((d0[0]*d0[1]) > fd0d0){continue;}
 516       if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
 517
 518       nD0++;
 519       if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
 520         fD0mass->Fill(fd0calc->InvMass(tN,tP));
 521       }
 522       else{
 523         fD0mass->Fill(fd0calc->InvMass(tP,tN));
 524       }
 525   }
 526   return nD0;
 527
 528 }
 529 int AliHLTD0Trigger::CheckTrackMC(AliExternalTrackParam* pt, AliExternalTrackParam* pn){
 530
 531   if(!fEvent){return 0;}
 532
 533   int lP = pt->GetLabel();
 534   int lN = pn->GetLabel();
 535
 536   if(lN>=0 && lP>=0){
 537
 538     int imP = (fEvent->Particle(lP))->GetFirstMother();
 539     int imN = (fEvent->Particle(lN))->GetFirstMother();
 540
 541     if(imP>=0 && imN>=0){
 542       TParticle * mP = fEvent->Particle(imP);
 543       TParticle * mN = fEvent->Particle(imN);
 544       if(fabs(mP->GetPdgCode())==421 && fabs(mN->GetPdgCode())==421 && imP == imN){
 545         return 2;
 546       }
 547       if(fabs(mP->GetPdgCode())==421 || fabs(mN->GetPdgCode())==421 && imP == imN){
 548         return 1;
 549       }
 550     }
 551   }
 552   return 0;
 553 }
 554 void AliHLTD0Trigger::GetOCDBObjectDescription( TMap* const targetMap)
 555 {
 556   // Get a list of OCDB object description.
 557   if (!targetMap) return;
 558   targetMap->Add(new TObjString("HLT/ConfigHLT/D0Trigger"),new TObjString("Object with default cuts for D0 reconstruction" ));
 559 }