Update and addition of LS analysis (Renu, Giacomo, Francesco)

diff --git a/PWG3/vertexingHF/AddTaskDplus.C b/PWG3/vertexingHF/AddTaskDplus.C
index fe237feae4196701affc5c329efe1c9415fddbd9..a389d7a8e2f5b1dc631180aa8bca9f0496f5c7f2 100644 (file)
--- a/PWG3/vertexingHF/AddTaskDplus.C
+++ b/PWG3/vertexingHF/AddTaskDplus.C
@@ -1,12 +1,11 @@
 AliAnalysisTaskSEDplus *AddTaskDplus(Bool_t storeNtuple=kFALSE)
 {
   //                                                                                                                                    
 AliAnalysisTaskSEDplus *AddTaskDplus(Bool_t storeNtuple=kFALSE)
 {
   //                                                                                                                                    

-  // Test macro for the AliAnalysisTaskSE for heavy-flavour candidates                                                                  
-  // association with MC truth (using MC info in AOD)                                                                                   
-                                                                                               
-  //                                                                                                                                    
-
+  // Test macro for the AliAnalysisTaskSE for D+ candidates 

 
 

+  //Invariant mass histogram and                                                 
+  // association with MC truth (using MC info in AOD)                                                                                   
+  //  R. Bala, bala@to.infn.it                                                                                                                                  

   // Get the pointer to the existing analysis manager via the static access method.                                                     
   //==============================================================================                                                      
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   // Get the pointer to the existing analysis manager via the static access method.                                                     
   //==============================================================================                                                      
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();


@@ -18,6 +17,8 @@ AliAnalysisTaskSEDplus *AddTaskDplus(Bool_t storeNtuple=kFALSE)
 
   // Aanalysis task                                                                                                                     
   AliAnalysisTaskSEDplus *dplusTask = new AliAnalysisTaskSEDplus("DplusAnalysis",storeNtuple);
 
   // Aanalysis task                                                                                                                     
   AliAnalysisTaskSEDplus *dplusTask = new AliAnalysisTaskSEDplus("DplusAnalysis",storeNtuple);

+  dplusTask->SetReadMC(kTRUE);
+  dplusTask->SetDoLikeSign(kTRUE);

   dplusTask->SetDebugLevel(0);
   mgr->AddTask(dplusTask);
 
   dplusTask->SetDebugLevel(0);
   mgr->AddTask(dplusTask);
 


@@ -32,11 +33,8 @@ AliAnalysisTaskSEDplus *AddTaskDplus(Bool_t storeNtuple=kFALSE)
     AliAnalysisDataContainer *coutputDplus2 = mgr->CreateContainer("coutputDplus2",TNtuple::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                                 "InvMassDplus_nt1.root");
     AliAnalysisDataContainer *coutputDplus2 = mgr->CreateContainer("coutputDplus2",TNtuple::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                                 "InvMassDplus_nt1.root");

-    AliAnalysisDataContainer *coutputDplus3 = mgr->CreateContainer("coutputDplus3",TNtuple::Class(),
-                                                           AliAnalysisManager::kOutputContainer,
-                                                                "InvMassDplus_nt2.root");
+

     coutputDplus2->SetSpecialOutput();
     coutputDplus2->SetSpecialOutput();

-    coutputDplus3->SetSpecialOutput();

   }
   mgr->ConnectInput(dplusTask,0,mgr->GetCommonInputContainer());
 
   }
   mgr->ConnectInput(dplusTask,0,mgr->GetCommonInputContainer());
 


@@ -44,7 +42,6 @@ AliAnalysisTaskSEDplus *AddTaskDplus(Bool_t storeNtuple=kFALSE)
   
   if(storeNtuple){
     mgr->ConnectOutput(dplusTask,2,coutputDplus2);
   
   if(storeNtuple){
     mgr->ConnectOutput(dplusTask,2,coutputDplus2);

-    mgr->ConnectOutput(dplusTask,3,coutputDplus3);

   }
   return dplusTask;
 }
   }
   return dplusTask;
 }

diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
index 6de97168cfa105aef7585762dd44b22976fea088..8acfa3a57af1a9679909d72a8e273cd510b56dbf 100644 (file)
--- a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
+++ b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
@@ -17,7 +17,9 @@
 //
 // AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
 // decay candidates with the MC truth.
 //
 // AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
 // decay candidates with the MC truth.

-// Author: Renu Bala, bala@to.infn.it
+// Authors: Renu Bala, bala@to.infn.it
+// F. Prino, prino@to.infn.it
+// G. Ortona, ortona@to.infn.it

 /////////////////////////////////////////////////////////////
 
 #include <TClonesArray.h>
 /////////////////////////////////////////////////////////////
 
 #include <TClonesArray.h>


@@ -25,6 +27,7 @@
 #include <TList.h>
 #include <TString.h>
 #include <TH1F.h>
 #include <TList.h>
 #include <TString.h>
 #include <TH1F.h>

+#include <TDatabasePDG.h>

 
 #include "AliAnalysisManager.h"
 #include "AliAODHandler.h"
 
 #include "AliAnalysisManager.h"
 #include "AliAODHandler.h"


@@ -45,9 +48,18 @@ ClassImp(AliAnalysisTaskSEDplus)
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
 AliAnalysisTaskSE(),
 fOutput(0), 
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
 AliAnalysisTaskSE(),
 fOutput(0), 

+fHistNEvents(0),

 fNtupleDplus(0),
 fNtupleDplus(0),

-fNtupleDplusbackg(0),
+fHistLS(0),
+fHistLStrip(0),
+fHistLStripcut(0),
+fHistOS(0),
+fHistOSbkg(0),
+fHistLSnoweight(0),
+fHistLSsinglecut(0),

 fFillNtuple(kFALSE),
 fFillNtuple(kFALSE),

+fReadMC(kFALSE),
+fDoLS(kFALSE),

 fVHF(0)
 {
   // Default constructor
 fVHF(0)
 {
   // Default constructor


@@ -56,10 +68,19 @@ fVHF(0)
 //________________________________________________________________________
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
 AliAnalysisTaskSE(name),
 //________________________________________________________________________
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
 AliAnalysisTaskSE(name),

-fOutput(0), 
+fOutput(0),
+fHistNEvents(0),

 fNtupleDplus(0),
 fNtupleDplus(0),

-fNtupleDplusbackg(0),
+fHistLS(0),
+fHistLStrip(0),
+fHistLStripcut(0),
+fHistOS(0),
+fHistOSbkg(0),
+fHistLSnoweight(0),
+fHistLSsinglecut(0),

 fFillNtuple(fillNtuple),
 fFillNtuple(fillNtuple),

+fReadMC(kFALSE),
+fDoLS(kFALSE),

 fVHF(0)
 {
   // Default constructor
 fVHF(0)
 {
   // Default constructor


@@ -70,8 +91,6 @@ fVHF(0)
   if(fFillNtuple){
     // Output slot #2 writes into a TNtuple container
     DefineOutput(2,TNtuple::Class());  //My private output
   if(fFillNtuple){
     // Output slot #2 writes into a TNtuple container
     DefineOutput(2,TNtuple::Class());  //My private output

-    // Output slot #3 writes into a TNtuple container
-    DefineOutput(3,TNtuple::Class());  //My private output

   }
 }
 
   }
 }
 


@@ -83,10 +102,15 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     delete fOutput;
     fOutput = 0;
   }
     delete fOutput;
     fOutput = 0;
   }

+  if(fNtupleDplus){
+    delete fNtupleDplus;
+    fNtupleDplus=0;
+  }

   if (fVHF) {
     delete fVHF;
     fVHF = 0;
   }
   if (fVHF) {
     delete fVHF;
     fVHF = 0;
   }

+

 }  
 
 //________________________________________________________________________
 }  
 
 //________________________________________________________________________


@@ -147,13 +171,45 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     fOutput->Add(hstc);
     fOutput->Add(hbtc);
   }
     fOutput->Add(hstc);
     fOutput->Add(hbtc);
   }

-
+  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
+  fHistNEvents->Sumw2();
+  fHistNEvents->SetMinimum(0);
+  fOutput->Add(fHistNEvents);
+  if(fDoLS){
+    Double_t massDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();

   
   

+    fHistLS = new TH1F("fHistLS","fHistLS",100,massDplus-0.2,massDplus+0.2);
+    fHistLS->Sumw2();
+    fOutput->Add(fHistLS);
+
+    fHistLStrip = new TH1F("fHistLStrip","fHistLStrip",100,massDplus-0.2,massDplus+0.2);
+    fHistLStrip->Sumw2();
+    fOutput->Add(fHistLStrip);
+
+    fHistLStripcut = new TH1F("fHistLStripcut","fHistLStripcut",100,massDplus-0.2,massDplus+0.2);
+    fHistLStripcut->Sumw2();
+    fOutput->Add(fHistLStripcut);
+
+    fHistOS = new TH1F("fHistOS","fHistOS",100,massDplus-0.2,massDplus+0.2);
+    fHistOS->Sumw2();
+    fOutput->Add(fHistOS);
+
+    fHistOSbkg = new TH1F("fHistOSbkg","fHistOSbkg",100,massDplus-0.2,massDplus+0.2);
+    fHistOSbkg->Sumw2();
+    fOutput->Add(fHistOSbkg);
+
+    fHistLSnoweight = new TH1F("fHistLSnoweight","fHistLSnoweight",100,massDplus-0.2,massDplus+0.2);
+    fHistLSnoweight->Sumw2();
+    fOutput->Add(fHistLSnoweight);
+
+    fHistLSsinglecut = new TH1F("fHistLSsinglecut","fHistLSsinglecut",100,massDplus-0.2,massDplus+0.2);
+    fHistLSsinglecut->Sumw2();
+    fOutput->Add(fHistLSsinglecut);
+  }
+

   if(fFillNtuple){
     OpenFile(2); // 2 is the slot number of the ntuple
   if(fFillNtuple){
     OpenFile(2); // 2 is the slot number of the ntuple

-    fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:PtK:Ptpi2:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2");
-    OpenFile(3); // 3 is the slot number of the ntuple
-    fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
+    fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2");

   }
 
   return;
   }
 
   return;


@@ -168,9 +224,13 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   
   
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
   
   
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());

+  fHistNEvents->Fill(0); // count event
+  // Post the data already here
+  PostData(1,fOutput);
+  

 
   TClonesArray *array3Prong = 0;
 
   TClonesArray *array3Prong = 0;

-
+  TClonesArray *arrayLikeSign =0;

   if(!aod && AODEvent() && IsStandardAOD()) {
     // In case there is an AOD handler writing a standard AOD, use the AOD 
     // event in memory rather than the input (ESD) event.    
   if(!aod && AODEvent() && IsStandardAOD()) {
     // In case there is an AOD handler writing a standard AOD, use the AOD 
     // event in memory rather than the input (ESD) event.    


@@ -183,39 +243,52 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
       AliAODEvent *aodFromExt = ext->GetAOD();
       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
       AliAODEvent *aodFromExt = ext->GetAOD();
       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");

+      arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");

     }
   } else {
     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
     }
   } else {
     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");

+    arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");

   }
 
   if(!array3Prong) {
     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
     return;
   }
   }
 
   if(!array3Prong) {
     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
     return;
   }

+  if(!arrayLikeSign) {
+    printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
+    // do in any case the OS analysis.
+  }
+
+ 
+  TClonesArray *arrayMC=0;
+  AliAODMCHeader *mcHeader=0;

 
   // AOD primary vertex
   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
   //    vtx1->Print();
   
   // load MC particles
 
   // AOD primary vertex
   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
   //    vtx1->Print();
   
   // load MC particles

-  TClonesArray *arrayMC = 
-    (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
-  if(!arrayMC) {
-    printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
-    return;
-  }
+  if(fReadMC){
+    
+    arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+    if(!arrayMC) {
+      printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
+      return;
+    }

     
   // load MC header
     
   // load MC header

-  AliAODMCHeader *mcHeader = 
-    (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
-  if(!mcHeader) {
-    printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
-    return;
+    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+    if(!mcHeader) {
+      printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
+      return;
+    }

   }
   }

+  

   Int_t n3Prong = array3Prong->GetEntriesFast();
   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
   Int_t n3Prong = array3Prong->GetEntriesFast();
   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);

-
-
+  
+  
+  Int_t nOS=0;

   Int_t pdgDgDplustoKpipi[3]={321,211,211};
   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
   Int_t pdgDgDplustoKpipi[3]={321,211,211};
   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {


@@ -227,6 +300,7 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       d->SetOwnPrimaryVtx(vtx1);
       unsetvtx=kTRUE;
     }
       d->SetOwnPrimaryVtx(vtx1);
       unsetvtx=kTRUE;
     }

+

     if(d->SelectDplus(fVHF->GetDplusCuts())) {
       Int_t iPtBin=0;
       Double_t ptCand = d->Pt();
     if(d->SelectDplus(fVHF->GetDplusCuts())) {
       Int_t iPtBin=0;
       Double_t ptCand = d->Pt();


@@ -234,26 +308,55 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        iPtBin=0;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;
        iPtBin=0;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;

-       cutsDplus[10]=0.979;
+       cutsDplus[9]=0.979;
+       cutsDplus[10]=0.0000055;

       }
       else if(ptCand>2. && ptCand<3){ 
        iPtBin=1;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.991;
       }
       else if(ptCand>2. && ptCand<3){ 
        iPtBin=1;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.991;

+       cutsDplus[10]=0.000005;

       }else if(ptCand>3. && ptCand<5){ 
        iPtBin=2;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.9955;
       }else if(ptCand>3. && ptCand<5){ 
        iPtBin=2;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.9955;

+       cutsDplus[10]=0.0000035;

       }else{
        iPtBin=3;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.997;
       }else{
        iPtBin=3;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
        cutsDplus[9]=0.997;

+       cutsDplus[10]=0.000001;

       }
       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
       }
       Bool_t passTightCuts=d->SelectDplus(cutsDplus);

-      Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
+      Int_t labDp=-1;
+      Float_t deltaPx=0.;
+      Float_t deltaPy=0.;
+      Float_t deltaPz=0.;
+      Float_t truePt=0.;
+      Float_t xDecay=0.;
+      Float_t yDecay=0.;
+      Float_t zDecay=0.;
+      Float_t pdgCode=-2;
+      if(fReadMC){
+       labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
+       if(labDp>=0){
+         AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
+         AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
+         deltaPx=partDp->Px()-d->Px();
+         deltaPy=partDp->Py()-d->Py();
+         deltaPz=partDp->Pz()-d->Pz();
+         truePt=partDp->Pt();
+         xDecay=dg0->Xv();       
+         yDecay=dg0->Yv();       
+         zDecay=dg0->Zv();
+         pdgCode=TMath::Abs(partDp->GetPdgCode());
+       }else{
+         pdgCode=-1;
+       }
+      }

       Double_t invMass=d->InvMassDplus();
 
       TString hisNameA(Form("hMassPt%d",iPtBin));
       Double_t invMass=d->InvMassDplus();
 
       TString hisNameA(Form("hMassPt%d",iPtBin));


@@ -268,90 +371,182 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
       }
 
        ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
       }
 

+      Float_t tmp[22];
+      if(fFillNtuple){           
+       tmp[0]=pdgCode;
+       tmp[1]=deltaPx;
+       tmp[2]=deltaPy;
+       tmp[3]=deltaPz;
+       tmp[4]=truePt;
+       tmp[5]=xDecay;    
+       tmp[6]=yDecay;    
+       tmp[7]=zDecay;    
+       tmp[8]=d->PtProng(0);
+       tmp[9]=d->PtProng(1);
+       tmp[10]=d->PtProng(2);
+       tmp[11]=d->Pt();
+       tmp[12]=d->CosPointingAngle();
+       tmp[13]=d->DecayLength();
+       tmp[14]=d->Xv();
+       tmp[15]=d->Yv();
+       tmp[16]=d->Zv();
+       tmp[17]=d->InvMassDplus();
+       tmp[18]=d->GetSigmaVert();
+       tmp[19]=d->Getd0Prong(0);
+       tmp[20]=d->Getd0Prong(1);
+       tmp[21]=d->Getd0Prong(2);         
+       fNtupleDplus->Fill(tmp);
+       PostData(2,fNtupleDplus);
+      }

 
 

-      if(labDp>=0) {
-       ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
-       if(passTightCuts){
-         ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
-       }
-       if(fFillNtuple){
-         AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
-         AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));         
-         Float_t tmp[20];
-         tmp[0]=TMath::Abs(partDp->GetPdgCode());
-         tmp[1]=partDp->Px()-d->Px();
-         tmp[2]=partDp->Py()-d->Py();
-         tmp[3]=partDp->Pz()-d->Pz();
-         tmp[4]=d->PtProng(0);
-         tmp[5]=d->PtProng(1);
-         tmp[6]=d->PtProng(2);
-         tmp[7]=d->Pt();
-         tmp[8]=partDp->Pt();
-         tmp[9]=d->CosPointingAngle();
-         tmp[10]=d->DecayLength();
-         tmp[11]=dg0->Xv();
-         tmp[12]=d->Xv();
-         tmp[13]=d->Yv();
-         tmp[14]=d->Zv();
-         tmp[15]=d->InvMassDplus();
-         tmp[16]=d->GetSigmaVert();
-         tmp[17]=d->Getd0Prong(0);
-         tmp[18]=d->Getd0Prong(1);
-         tmp[19]=d->Getd0Prong(2);       
-         fNtupleDplus->Fill(tmp);
-         PostData(2,fNtupleDplus);
-       }
-      }else{
-       ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
-       if(passTightCuts){
-         ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
-       }
-       if(fFillNtuple){
-         Float_t tmpbkg[14];
-         tmpbkg[0]=d->PtProng(0);
-         tmpbkg[1]=d->PtProng(1);
-         tmpbkg[2]=d->PtProng(2);
-         tmpbkg[3]=d->Pt();
-         tmpbkg[4]=d->CosPointingAngle();
-         tmpbkg[5]=d->DecayLength();
-         tmpbkg[6]=d->Xv();
-         tmpbkg[7]=d->Yv();
-         tmpbkg[8]=d->Zv();
-         tmpbkg[9]=d->InvMassDplus();
-         tmpbkg[10]=d->GetSigmaVert();
-         tmpbkg[11]=d->Getd0Prong(0);
-         tmpbkg[12]=d->Getd0Prong(1);
-         tmpbkg[13]=d->Getd0Prong(2);
-         fNtupleDplusbackg->Fill(tmpbkg);          
-         PostData(3,fNtupleDplusbackg);
+      if(fReadMC){
+       if(labDp>=0) {
+         ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
+         if(passTightCuts){
+           ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
+         }
+         
+       }else{
+         ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
+         if(passTightCuts){
+           ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
+         }
+         fHistOSbkg->Fill(invMass);

        }
       }
        }
       }

-      PostData(1,fOutput);    
+
+      fHistOS->Fill(invMass);
+      nOS++;

     }
     if(unsetvtx) d->UnsetOwnPrimaryVtx();
   }
     }
     if(unsetvtx) d->UnsetOwnPrimaryVtx();
   }

-  
+ 
+  //start LS analysis
+  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
+   
+  PostData(1,fOutput);    

   return;
 }
 
   return;
 }
 

-//________________________________________________________________________
+
+
+//_________________________________________________________________
+void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
+
+/*
+ * Fill the Like Sign histograms
+ */
+
+
+
+
+  //count pos/neg tracks
+  Int_t nPosTrks=0,nNegTrks=0;
+  //counter for particles passing single particle cuts
+  Int_t nspcplus=0;
+  Int_t nspcminus=0;
+
+  for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
+    AliAODTrack *track = aod->GetTrack(it);
+    if(track->Charge()>0){
+      nPosTrks++;
+      if(track->Pt()>=0.4){
+       nspcplus++;
+      }
+    }else{      
+      nNegTrks++;
+      if(track->Pt()>=0.4){
+       nspcminus++;
+      }
+    }
+  }
+
+  Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
+
+  Int_t nDplusLS=0;
+  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
+ 
+  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
+    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
+    Bool_t unsetvtx=kFALSE;
+    if(!d->GetOwnPrimaryVtx()) {
+      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+      unsetvtx=kTRUE;
+    }
+    if(d->SelectDplus(fVHF->GetDplusCuts()))nDplusLS++;
+    if(unsetvtx) d->UnsetOwnPrimaryVtx();
+  }
+
+  Float_t wei2=0;
+  if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
+  Float_t wei3=0;
+  if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
+
+ // loop over  sign candidates
+  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
+    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
+    Bool_t unsetvtx=kFALSE;
+    if(!d->GetOwnPrimaryVtx()) {
+      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+      unsetvtx=kTRUE;
+    }
+ 
+    if(d->SelectDplus(fVHF->GetDplusCuts())){
+      Int_t sign= d->GetCharge();
+      Float_t wei1=1;
+      Float_t wei4=1;
+      if(sign>0&&nPosTrks>2&&nspcplus>2){//wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
+       wei1=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
+       wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
+      }
+      if(sign<0&&nNegTrks>2&&nspcminus>2){
+       wei1=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
+       wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
+      }
+      Double_t invMass=d->InvMassDplus();
+      fHistLS->Fill(invMass,wei1);
+      fHistLSnoweight->Fill(invMass);
+      fHistLStrip->Fill(invMass,wei2);
+      fHistLStripcut->Fill(invMass,wei3);
+      fHistLSsinglecut->Fill(invMass,wei4);
+     }
+    if(unsetvtx) d->UnsetOwnPrimaryVtx();
+  }
+  
+  //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
+  //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
+
+  //  printf("LS analysis...done\n");
+
+}
+
+//_________________________________________________________________
+
+

 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
 {
   // Terminate analysis
   //
   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
 {
   // Terminate analysis
   //
   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");

-

   fOutput = dynamic_cast<TList*> (GetOutputData(1));
   if (!fOutput) {     
     printf("ERROR: fOutput not available\n");
     return;
   }
   fOutput = dynamic_cast<TList*> (GetOutputData(1));
   if (!fOutput) {     
     printf("ERROR: fOutput not available\n");
     return;
   }

-
+  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));

   if(fFillNtuple){
     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
   if(fFillNtuple){
     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));

-    fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
+  }
+  fHistOS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOS"));
+  fHistOSbkg =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOSbkg"));
+  if(fDoLS){
+    fHistLS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLS"));
+    fHistLStrip =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStrip"));
+    fHistLStripcut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStripcut"));
+    fHistLSnoweight =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSnoweight"));
+    fHistLSsinglecut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSsinglecut"));

   }
 
   }
 

+

   return;
 }
   return;
 }

-

diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.h b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.h
index fc09a32ada9c6a84a419489f9e5a190d1d798df0..66b65a3614bb4a36416c3bfb27c36f93ff4cbc47 100644 (file)
--- a/PWG3/vertexingHF/AliAnalysisTaskSEDplus.h
+++ b/PWG3/vertexingHF/AliAnalysisTaskSEDplus.h
@@ -6,9 +6,10 @@
 
 //*************************************************************************
 // Class AliAnalysisTaskSEDplus
 
 //*************************************************************************
 // Class AliAnalysisTaskSEDplus

-// AliAnalysisTaskSE for the comparison of heavy-flavour decay candidates
+// AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and 
+//comparison of heavy-flavour decay candidates

 // to MC truth (kinematics stored in the AOD)
 // to MC truth (kinematics stored in the AOD)

-
+// Renu Bala, bala@to.infn.it

 //*************************************************************************
 
 #include <TROOT.h>
 //*************************************************************************
 
 #include <TROOT.h>


@@ -26,24 +27,34 @@ class AliAnalysisTaskSEDplus : public AliAnalysisTaskSE
   AliAnalysisTaskSEDplus();
   AliAnalysisTaskSEDplus(const char *name, Bool_t fillNtuple=kFALSE);
   virtual ~AliAnalysisTaskSEDplus();
   AliAnalysisTaskSEDplus();
   AliAnalysisTaskSEDplus(const char *name, Bool_t fillNtuple=kFALSE);
   virtual ~AliAnalysisTaskSEDplus();

-
-
+  void SetReadMC(Bool_t readMC=kTRUE){fReadMC=readMC;}
+  void SetDoLikeSign(Bool_t dols=kTRUE){fDoLS=dols;}

   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void Init();
   virtual void LocalInit() {Init();}
   virtual void UserExec(Option_t *option);
   virtual void Terminate(Option_t *option);
   // Implementation of interface methods
   virtual void UserCreateOutputObjects();
   virtual void Init();
   virtual void LocalInit() {Init();}
   virtual void UserExec(Option_t *option);
   virtual void Terminate(Option_t *option);

-
-  
+  void LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS);
+    

  private:
 
   AliAnalysisTaskSEDplus(const AliAnalysisTaskSEDplus &source);
   AliAnalysisTaskSEDplus& operator=(const AliAnalysisTaskSEDplus& source); 
   TList   *fOutput; //! list send on output slot 0
  private:
 
   AliAnalysisTaskSEDplus(const AliAnalysisTaskSEDplus &source);
   AliAnalysisTaskSEDplus& operator=(const AliAnalysisTaskSEDplus& source); 
   TList   *fOutput; //! list send on output slot 0

+  TH1F    *fHistNEvents; //!hist. for No. of events

   TNtuple *fNtupleDplus; //! output ntuple
   TNtuple *fNtupleDplus; //! output ntuple

-  TNtuple *fNtupleDplusbackg; //! output ntuple
+  TH1F *fHistLS;     //! LS tripl normalized using all charged tracks in the event
+  TH1F *fHistLStrip; //! LS tripl normalized using OS vs LS number of triplets ratio
+  TH1F *fHistLStripcut;//! LS tripl normalized using OS vs LS number of triplets passing dplus cuts ratio
+  TH1F *fHistOS; //! OS triplets
+  TH1F *fHistOSbkg;//! only Background (need MC truth)
+  TH1F *fHistLSnoweight;//! LS triplets wo normalization
+  TH1F *fHistLSsinglecut;//! LS tripl normalized using tracks with Pt>0.4GeV
+  

   Bool_t fFillNtuple;   // flag for filling ntuple
   Bool_t fFillNtuple;   // flag for filling ntuple

+  Bool_t fReadMC;       //flag for access to MC
+  Bool_t fDoLS;        //flag for switching on/off Like Sign analysis

   AliAnalysisVertexingHF *fVHF;  // Vertexer heavy flavour (used to pass the cuts)
   
   ClassDef(AliAnalysisTaskSEDplus,4); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates
   AliAnalysisVertexingHF *fVHF;  // Vertexer heavy flavour (used to pass the cuts)
   
   ClassDef(AliAnalysisTaskSEDplus,4); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates