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 }