HLT/QA/tasks/AliAnalysisTaskD0Trigger.cxx

   1 // $Id$
   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 <st05886@alf.uib.no>,                *
   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
  19 /** @file  AliAnalysisTaskD0Trigger.cxx
  20     @author Gaute Ovrebekk
  21     @date
  22     @brief An analysis task for the D0 Trigger.
  23 */
  24
  25 class AliAnalysisTask;
  26 class AliAnalysisManager;
  27
  28 #include "TH1F.h"
  29 #include "TCanvas.h"
  30 #include "AliESDEvent.h"
  31 #include "AliESDInputHandler.h"
  32 #include "TObjArray.h"
  33 #include "TDatabasePDG.h"
  34 #include "AliESDVertex.h"
  35 #include "TMath.h"
  36 #include "AliExternalTrackParam.h"
  37 #include "AliAODVertex.h"
  38 #include "TDatabasePDG.h"
  39 #include "AliESDtrack.h"
  40 #include "TVector3.h"
  41 #include "AliVertexerTracks.h"
  42 #include "AliKFVertex.h"
  43 #include "TDatabasePDG.h"
  44 #include "TVector3.h"
  45
  46 #include "AliAnalysisTaskD0Trigger.h"
  47
  48 ClassImp(AliAnalysisTaskD0Trigger)
  49
  50 AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger()
  51 :
  52 AliAnalysisTaskSE()
  53   , fOutputList(0)
  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   , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  62   , fD0massHLT(NULL)
  63   , fD0ptHLT(NULL)
  64   , fD0massOFF(NULL)
  65   , fD0ptOFF(NULL)
  66   , fPos()
  67   , fNeg()
  68   , ftwoTrackArray(NULL)
  69   , fTotalD0HLT(0)
  70   , fTotalD0OFF(0)
  71   , fField(0)
  72   , fNevents(0)
  73   , fuseKF(false)
  74 {
  75   // Constructor
  76   // Define input and output slots here
  77   // Input slot #0 works with a TChain
  78   // DefineInput(0, TChain::Class());
  79   // Output slot #0 writes into a TH1 container
  80
  81   //DefineOutput(1, TList::Class());
  82 }
  83
  84 AliAnalysisTaskD0Trigger::AliAnalysisTaskD0Trigger(const char *name,float cuts[7])
  85 :
  86 AliAnalysisTaskSE(name)
  87   , fOutputList(0)
  88   , fPtMin(cuts[0])
  89   , fdca(cuts[1])
  90   , finvMass(cuts[2])
  91   , fcosThetaStar(cuts[3])
  92   , fd0(cuts[4])
  93   , fd0d0(cuts[5])
  94   , fcosPoint(cuts[6])
  95   , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  96   , fD0massHLT(NULL)
  97   , fD0ptHLT(NULL)
  98   , fD0massOFF(NULL)
  99   , fD0ptOFF(NULL)
 100   , fPos()
 101   , fNeg()
 102   , ftwoTrackArray(NULL)
 103   , fTotalD0HLT(0)
 104   , fTotalD0OFF(0)
 105   , fField(0)
 106   , fNevents(0)
 107   , fuseKF(false)
 108 {
 109   // Constructor
 110   // Define input and output slots here
 111   // Input slot #0 works with a TChain
 112   // DefineInput(0, TChain::Class());
 113   // Output slot #0 writes into a TH1 container
 114
 115   DefineOutput(1, TList::Class());
 116 }
 117 void AliAnalysisTaskD0Trigger::UserCreateOutputObjects(){
 118   // Create histograms
 119
 120   OpenFile(1);
 121
 122   fOutputList = new TList();
 123   fOutputList->SetName(GetName());
 124
 125   fD0massHLT = new TH1F("hMassHLT","D^{0} mass plot from HLT reconstruction",100,1.7,2);
 126   fD0ptHLT = new TH1F("hPtHLT","D^{0} Pt plot HLT reconstruction",20,0,20);
 127
 128   fD0massOFF = new TH1F("hMassOFF","D^{0} mass plot Offline reconstruction",100,1.7,2);
 129   fD0ptOFF = new TH1F("hPtOFF","D^{0} Pt plot Offline reconstruction",20,0,20);
 130
 131   fOutputList->Add(fD0massHLT);
 132   fOutputList->Add(fD0ptHLT);
 133   fOutputList->Add(fD0massOFF);
 134   fOutputList->Add(fD0ptOFF);
 135
 136 }
 137
 138 void AliAnalysisTaskD0Trigger::NotifyRun(){
 139   // see header file of AliAnalysisTask for documentation
 140 }
 141
 142 void AliAnalysisTaskD0Trigger::UserExec(Option_t *){
 143   // see header file of AliAnalysisTask for documentation
 144
 145   AliESDEvent *esdOFF = dynamic_cast<AliESDEvent*>(InputEvent());
 146
 147   if(!esdOFF){
 148     Printf("ERROR: fESD not available");
 149     return;
 150   }
 151
 152   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(fInputHandler);
 153   AliESDEvent *esdHLT = NULL;
 154   if(esdH) esdHLT = esdH->GetHLTEvent();
 155
 156   if(!esdHLT){
 157     Printf("ERROR: HLTesd not available");
 158     return;
 159   }
 160
 161   fNevents++;
 162
 163   ftwoTrackArray = new TObjArray(2);
 164   Int_t nD0=0;
 165
 166   fField = esdHLT->GetMagneticField();
 167   const AliESDVertex* pvHLT = esdHLT->GetPrimaryVertexTracks();
 168   AliESDVertex *pVertexHLT =  new AliESDVertex(*pvHLT);
 169   if(pVertexHLT->GetNContributors()<2){
 170     //Printf("ERROR: Contributiors to HLT Primary vertex is to low or not set");
 171     return;
 172   }
 173
 174   for(Int_t it=0;it<esdHLT->GetNumberOfTracks();it++){
 175     SingleTrackSelect(esdHLT->GetTrack(it),pVertexHLT);
 176   }
 177
 178   RecD0(nD0,pVertexHLT,true);
 179   fTotalD0HLT+=nD0;
 180
 181   fPos.clear();
 182   fNeg.clear();
 183
 184   nD0=0;
 185
 186   fField = esdOFF->GetMagneticField();
 187   const AliESDVertex* pvOFF = esdOFF->GetPrimaryVertexTracks();
 188   AliESDVertex *pVertexOFF =  new AliESDVertex(*pvOFF);
 189   if(pVertexOFF->GetNContributors()<2){
 190     //Printf("ERROR: Contributiors to OFFLINE Primary vertex is to low or not set");
 191     return;
 192   }
 193   for(Int_t it=0;it<esdOFF->GetNumberOfTracks();it++){
 194     SingleTrackSelect(esdOFF->GetTrack(it),pVertexOFF);
 195   }
 196
 197   RecD0(nD0,pVertexOFF,false);
 198   fTotalD0OFF+=nD0;
 199
 200   fPos.clear();
 201   fNeg.clear();
 202
 203   // Post output data.
 204   PostData(1, fOutputList);
 205   ftwoTrackArray->Clear();
 206   delete pVertexHLT;
 207   delete pVertexOFF;
 208 }
 209
 210 void AliAnalysisTaskD0Trigger::Terminate(Option_t *){
 211   Printf("Event Number: %d",fNevents);
 212   Printf("Total Number of D0 foundfor HLT: %d",fTotalD0HLT);
 213   Printf("Total Number of D0 found for OFFLINE: %d",fTotalD0OFF);
 214 }
 215
 216 void AliAnalysisTaskD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, AliESDVertex *pV){
 217   // Offline har || på disse kuttene på de to henfallsproduktene
 218   Double_t pv[3];
 219   pV->GetXYZ(pv);
 220
 221   if(t->Pt()<fPtMin){return;}
 222   if(TMath::Abs(t->GetD(pv[0],pv[1],fField)) > fd0){return;}
 223
 224   if(t->Charge()>0){
 225     fPos.push_back(t);
 226   }
 227   else{
 228     fNeg.push_back(t);
 229   }
 230 }
 231
 232 void AliAnalysisTaskD0Trigger::RecD0(Int_t& nD0, AliESDVertex *pV,bool isHLT){
 233   // Reconstructing D0
 234   Double_t starD0,starD0bar,xdummy,ydummy;
 235   Double_t d0[2];
 236   Double_t svpos[3];
 237   Double_t pvpos[3];
 238   ftwoTrackArray->Clear();
 239
 240   if(!pV){
 241     Printf("No Primary Vertex");
 242     return;
 243   }
 244   pV->GetXYZ(pvpos);
 245
 246   for(UInt_t ip=0;ip<fPos.size();ip++){
 247     AliExternalTrackParam *tP=fPos[ip];
 248     for(UInt_t in=0;in<fNeg.size();in++){
 249       AliExternalTrackParam *tN=fNeg[in];
 250
 251       tP->PropagateToDCA(pV,fField,kVeryBig);  //do I need this??????
 252       tN->PropagateToDCA(pV,fField,kVeryBig);
 253
 254       Double_t dcatPtN = tP->GetDCA(tN,fField,xdummy,ydummy);
 255       if(dcatPtN>fdca) { continue; }
 256
 257       ftwoTrackArray->AddAt(tP,0);
 258       ftwoTrackArray->AddAt(tN,1);
 259       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(ftwoTrackArray,fField,pV,fuseKF);
 260       if(!vertexp1n1) {
 261         ftwoTrackArray->Clear();
 262         continue;
 263       }
 264
 265       vertexp1n1->GetXYZ(svpos);
 266
 267       tP->PropagateToDCA(vertexp1n1,fField,kVeryBig);
 268       tN->PropagateToDCA(vertexp1n1,fField,kVeryBig);
 269
 270       if((TMath::Abs(InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((InvMass(tP,tN))-fD0PDG) > finvMass){continue;}
 271       CosThetaStar(tN,tP,starD0,starD0bar);
 272       if(TMath::Abs(starD0) > fcosThetaStar && TMath::Abs(starD0bar) > fcosThetaStar){continue;}
 273       d0[0] = tP->GetD(pvpos[0],pvpos[1],fField);
 274       d0[1] = tN->GetD(pvpos[0],pvpos[1],fField);
 275       if((d0[0]*d0[1]) > fd0d0){continue;}
 276       if(PointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
 277
 278       if(isHLT){
 279         fD0massHLT->Fill(InvMass(tN,tP));
 280         fD0massHLT->Fill(InvMass(tP,tN));
 281         fD0ptHLT->Fill(Pt(tP,tN));
 282       }
 283       else{
 284         fD0massOFF->Fill(InvMass(tN,tP));
 285         fD0massOFF->Fill(InvMass(tP,tN));
 286         fD0ptOFF->Fill(Pt(tP,tN));
 287       }
 288       nD0++;
 289       delete vertexp1n1;
 290     }
 291   }
 292 }
 293
 294 Double_t AliAnalysisTaskD0Trigger::InvMass(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
 295 {
 296   // Calculating the invariant mass
 297   Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 298   Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
 299
 300   Double_t energy[2];
 301   energy[1] = TMath::Sqrt(mK*mK+d1->GetP()*d1->GetP());
 302   energy[0] = TMath::Sqrt(mpi*mpi+d2->GetP()*d2->GetP());
 303
 304   Double_t p1[3],p2[3];
 305   d1->GetPxPyPz(p1);
 306   d2->GetPxPyPz(p2);
 307
 308   Double_t momTot2 = (p1[0]+p2[0])*(p1[0]+p2[0])+
 309                      (p1[1]+p2[1])*(p1[1]+p2[1])+
 310                      (p1[2]+p2[2])*(p1[2]+p2[2]);
 311
 312   return TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
 313
 314 }
 315
 316 void AliAnalysisTaskD0Trigger::CosThetaStar(AliExternalTrackParam* d1, AliExternalTrackParam* d2,Double_t &D0,Double_t &D0bar)
 317 {
 318   //Calculating the decay angle
 319   Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 320   Double_t mpi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 321   Double_t mK=TDatabasePDG::Instance()->GetParticle(321)->Mass();
 322
 323   Double_t pStar = TMath::Sqrt(TMath::Power(mD0*mD0-mK*mK-mpi*mpi,2.)-4.*mK*mK*mpi*mpi)/(2.*mD0);
 324
 325   Double_t px = d1->Px()+d2->Px();
 326   Double_t py = d1->Py()+d2->Py();
 327   Double_t pz = d1->Pz()+d2->Pz();
 328   Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
 329   Double_t energy = TMath::Sqrt(p*p+mD0*mD0);
 330
 331   Double_t beta = p/energy;
 332   Double_t gamma = energy/mD0;
 333
 334   Double_t qL;
 335   TVector3 mom(d1->Px(),d1->Py(),d1->Pz());
 336   TVector3 momD(px,py,pz);
 337   qL = mom.Dot(momD)/momD.Mag();
 338
 339   D0 = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+mK*mK))/pStar;
 340
 341   TVector3 mom2(d2->Px(),d2->Py(),d2->Pz());
 342   TVector3 momD2(px,py,pz);
 343   qL = mom2.Dot(momD2)/momD2.Mag();
 344
 345   D0bar = (qL/gamma-beta*TMath::Sqrt(pStar*pStar+mK*mK))/pStar;
 346
 347 }
 348
 349 Double_t AliAnalysisTaskD0Trigger::PointingAngle(AliExternalTrackParam* n, AliExternalTrackParam* p, Double_t *pv, Double_t *sv)
 350 {
 351   // Calcutating the pointing angle
 352   TVector3 mom(n->Px()+p->Px(),n->Py()+p->Py(),n->Pz()+p->Pz());
 353   TVector3 flight(sv[0]-pv[0],sv[1]-pv[1],sv[2]-pv[2]);
 354
 355   double pta = mom.Angle(flight);
 356
 357   return TMath::Cos(pta);
 358 }
 359
 360 AliAODVertex* AliAnalysisTaskD0Trigger::ReconstructSecondaryVertex(TObjArray *trkArray, Double_t b, const AliESDVertex *v, bool useKF)
 361 {
 362   // Finding the vertex of the two tracks in trkArray
 363   AliESDVertex *vertexESD = 0;
 364   AliAODVertex *vertexAOD = 0;
 365
 366   if(!useKF){
 367     AliVertexerTracks *vertexer = new AliVertexerTracks(b);
 368     AliESDVertex* vertex =  const_cast<AliESDVertex*>(v);
 369     vertexer->SetVtxStart(vertex);
 370     //if(isESD){vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);}
 371     UShort_t *id = new UShort_t[2];
 372     AliExternalTrackParam *t1 = (AliExternalTrackParam*) trkArray->At(0);
 373     AliExternalTrackParam *t2 = (AliExternalTrackParam*) trkArray->At(1);
 374     id[0]=(UShort_t) t1->GetID();
 375     id[1]=(UShort_t) t2->GetID();
 376     vertexESD = (AliESDVertex*)vertexer->VertexForSelectedTracks(trkArray,id);
 377     delete [] id;
 378     delete vertexer; vertexer=NULL;
 379
 380     if(!vertexESD) return vertexAOD;
 381
 382     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) {
 383       //AliDebug(2,"vertexing failed");
 384       delete vertexESD; vertexESD=NULL;
 385       delete vertex;
 386       return vertexAOD;
 387     }
 388   }
 389   else{
 390     AliKFParticle::SetField(b);
 391
 392     AliKFVertex vertexKF;
 393
 394     Int_t nTrks = trkArray->GetEntriesFast();
 395     for(Int_t i=0; i<nTrks; i++) {
 396       AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
 397       AliKFParticle daughterKF(*esdTrack,211);
 398       vertexKF.AddDaughter(daughterKF);
 399     }
 400     vertexESD = new AliESDVertex(vertexKF.Parameters(),
 401                                  vertexKF.CovarianceMatrix(),
 402                                  vertexKF.GetChi2(),
 403                                  vertexKF.GetNContributors());
 404   }
 405   // convert to AliAODVertex
 406   Double_t pos[3],cov[6],chi2perNDF;
 407   vertexESD->GetXYZ(pos); // position
 408   vertexESD->GetCovMatrix(cov); //covariance matrix
 409   chi2perNDF = vertexESD->GetChi2toNDF();
 410   //dispersion = vertexESD->GetDispersion();
 411   delete vertexESD; vertexESD=NULL;
 412
 413   Int_t nprongs= trkArray->GetEntriesFast();
 414   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 415
 416   return vertexAOD;
 417
 418 }
 419
 420 Double_t AliAnalysisTaskD0Trigger::Pt(AliExternalTrackParam* d1, AliExternalTrackParam* d2)
 421 {
 422   //Calculating pT of the two tracks
 423   Double_t p1[3],p2[3];
 424   d1->GetPxPyPz(p1);
 425   d2->GetPxPyPz(p2);
 426
 427   Double_t pt2 = (p1[0]+p2[0])*(p1[0]+p2[0]) + (p1[1]+p2[1])*(p1[1]+p2[1]);
 428
 429   return TMath::Sqrt(pt2);
 430 }