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