Add histograms imppar and massB
authorcbianchi <cbianchi@cern.ch>
Wed, 12 Nov 2014 18:03:57 +0000 (19:03 +0100)
committercbianchi <cbianchi@cern.ch>
Wed, 12 Nov 2014 18:06:53 +0000 (19:06 +0100)
PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx
PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.h
PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx
PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.h
PWGJE/FlavourJetTasks/macros/AddTaskDFilterAndCorrelations.C
PWGJE/FlavourJetTasks/macros/AddTasksFlavourJet.C

index e61d881..ac6abaa 100644 (file)
@@ -102,10 +102,14 @@ fhNtrArr(),
 fhNJetPerEv(),
 fhdeltaRJetTracks(),
 fhsJet(),
+fhImpPar(),
 fhNDPerEvNoJet(),
 fhptDPerEvNoJet(),
 fhNJetPerEvNoD(),
 fhPtJetPerEvNoD(),
+fhVtxX(),
+fhVtxY(),
+fhVtxZ(),
 fhsDstandalone(),
 fhInvMassptD(),
 fhDiffSideBand(),
@@ -133,6 +137,7 @@ fhDeltaRD(),
 fhDeltaRptDptj(),
 fhDeltaRptDptjB(),
 fhsDphiz()
+
 {
    //
    // Default ctor
@@ -189,10 +194,14 @@ fhNtrArr(),
 fhNJetPerEv(),
 fhdeltaRJetTracks(),
 fhsJet(),
+fhImpPar(),
 fhNDPerEvNoJet(),
 fhptDPerEvNoJet(),
 fhNJetPerEvNoD(),
 fhPtJetPerEvNoD(),
+fhVtxX(),
+fhVtxY(),
+fhVtxZ(),
 fhsDstandalone(),
 fhInvMassptD(),
 fhDiffSideBand(),
@@ -429,7 +438,14 @@ Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
    if(!iseventselected) return kFALSE;
    
    fhstat->Fill(1);
-
+   
+   const AliVVertex *vertex = aodEvent->GetPrimaryVertex();
+   Double_t vtx[3];
+   vertex->GetXYZ(vtx);
+   fhVtxX->Fill(vtx[0]);
+   fhVtxY->Fill(vtx[1]);
+   fhVtxZ->Fill(vtx[2]);
+   //Printf(">>>>>>>>.>>>Primary vertex %.3f,%.3f,%.3f",vtx[0], vtx[1],vtx[2]);
    //retrieve charm candidates selected
    Int_t candidates = 0;
    Int_t njets=fJetCont->GetNJets();
@@ -531,13 +547,30 @@ Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
    fhNtrArr->Fill(ntrarr);
    
    for(Int_t i=0;i<ntrarr;i++){
-      AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackCont->GetParticle(i));
+      AliAODTrack *jtrack=static_cast<AliAODTrack*>(fTrackCont->GetParticle(i));
       //reusing histograms
       if(!jtrack) continue;
       fhPtJetTrks->Fill(jtrack->Pt());
       fhPhiJetTrks->Fill(jtrack->Phi());
       fhEtaJetTrks->Fill(jtrack->Eta());
       fhEjetTrks->Fill(jtrack->E());
+      Double_t vDCAglobalxy;
+      //Double_t vDCAglobalz;
+      if(jtrack->IsGlobalConstrained()){
+      // retrieve impact parameter
+      Double_t pos[3];
+
+      jtrack->GetXYZ(pos);
+      
+      vDCAglobalxy  = pos[0] - vtx[0]; //should be impact parameter in transverse plane
+      //vDCAglobalz  = pos[1] - vtx[1]; //should be impact parameter in z direction
+      }
+      else{
+      vDCAglobalxy=jtrack->DCA(); 
+      //vDCAglobalz=jtrack->ZAtDCA();
+      }
+      fhImpPar->Fill(vDCAglobalxy,jtrack->Pt());
+      //Printf("Track position  %.3f,%.3f,%.3f",pos[0],pos[1],pos[2]);
    }
    
    
@@ -883,7 +916,7 @@ Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
    fhstat->SetNdivisions(1);
    fOutput->Add(fhstat);
    
-   const Int_t nbinsmass=200;
+   const Int_t nbinsmass=300;
    const Int_t nbinsptjet=500;
    const Int_t nbinsptD=100;
    const Int_t nbinsz=100;
@@ -937,6 +970,19 @@ Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
    fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
    fhNJetPerEv->Sumw2();
    
+   fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used 
+   //NB at the moment fillet with jet track imp par!!!
+   fhVtxX = new TH1F("hVtxX", "Vertex X;vtx_x",100,-0.5,0.5);
+   fhVtxY = new TH1F("hVtxY", "Vertex Y;vtx_y",100,-0.5,0.5);
+   fhVtxZ = new TH1F("hVtxZ", "Vertex Z;vtx_z",100,-10,10);
+   /*
+   if(fUseMCInfo){
+   //understand how to do this. At the moment fhImpPar is filled with the impact parameter of jet tracks even if it is written differently in the fhImpPar definition
+      fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,nbinsptD, ptDlims[0], ptDlims[1]); //same range of pt of the D, but pt daughter used
+      fOutput->Add(fhImpParB);
+      
+   }
+   */
    
    fOutput->Add(fhEjetTrks);
    fOutput->Add(fhPhiJetTrks);
@@ -949,7 +995,10 @@ Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
    fOutput->Add(fhdeltaRJetTracks);
    fOutput->Add(fhNtrArr);
    fOutput->Add(fhNJetPerEv);
-   
+   fOutput->Add(fhImpPar);
+   fOutput->Add(fhVtxX);
+   fOutput->Add(fhVtxY); 
+   fOutput->Add(fhVtxZ);
    if(fJetOnlyMode && fSwitchOnSparses){
       //thnsparse for jets
       const Int_t nAxis=6;   
@@ -1282,7 +1331,7 @@ void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDe
    }
    
    
-   Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());   
+   //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());   
    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
    if(isselected==1 || isselected==3) {
       
@@ -1690,6 +1739,7 @@ Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *the
    Int_t nchtrk=thejet->GetNumberOfTracks();
    for(Int_t itk=0;itk<nchtrk;itk++){
       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
+      if(!tkjet) continue;
       Int_t idtkjet=tkjet->GetID();
       if(fillH) fhIDjetTracks->Fill(idtkjet);
       for(Int_t id=0;id<ndaugh;id++){
index 7ef4c1f..9f58568 100644 (file)
@@ -195,11 +195,16 @@ private:
    TH1F* fhNJetPerEv;               //!
    TH1F* fhdeltaRJetTracks;         //!
    THnSparse* fhsJet; //available in jet only mode
+   TH2F* fhImpPar;                  //!
+   //TH2F* fhImpParB;                 //!
    // event characteristics;        
    TH1F* fhNDPerEvNoJet;            //!
    TH1F* fhptDPerEvNoJet;           //!
    TH1F* fhNJetPerEvNoD;            //!
    TH1F* fhPtJetPerEvNoD;           //!
+   TH1F* fhVtxX;                    //!
+   TH1F* fhVtxY;                    //!
+   TH1F* fhVtxZ;                    //!
    //D mesons                       
    THnSparse* fhsDstandalone;       //!
    TH2F* fhInvMassptD;              //!
index 9cfe74b..92adbf9 100644 (file)
@@ -61,8 +61,11 @@ fCuts(0),
 fMinMass(0.),
 fMaxMass(0.),
 fCandidateArray(0),
-fSideBandArray(0)
-
+fSideBandArray(0),
+fhImpPar(0),
+fhImpParB(0),
+fhInvMassS(0),
+fhInvMassB(0)
 {
    //
    // Default constructor
@@ -88,7 +91,11 @@ fCuts(cuts),
 fMinMass(0.),
 fMaxMass(0.),
 fCandidateArray(0),
-fSideBandArray(0)
+fSideBandArray(0),
+fhImpPar(0),
+fhImpParB(0),
+fhInvMassS(0),
+fhInvMassB(0)
 {
    //
    // Constructor. Initialization of Inputs and Outputs
@@ -291,7 +298,7 @@ void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
    hstat->Fill(1);
    
    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
-   if(fUseReco) hstat->Fill(2, nD);
+   if(!fUseMCInfo) hstat->Fill(2, nD);
    
    //D* and D0 prongs needed to MatchToMC method
    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
@@ -335,8 +342,14 @@ void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
         if (fCandidateType==kD0toKpi) 
            mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
         
-        if (mcLabel<=0) isMCBkg=kTRUE;
-        else hstat->Fill(2);
+        if (mcLabel<=0) {
+           isMCBkg=kTRUE;
+           hstat->Fill(5);
+        }
+        else {
+           hstat->Fill(2);
+           
+        }
         if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
                 
         if (isMCBkg) smcTruth="B";
@@ -369,19 +382,34 @@ void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
       if (!isSelected) continue;
       
+      Int_t nprongs=charmCand->GetNProngs();
+      
       //for data and MC signal fill fCandidateArray
       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
         // for data or MC with the requirement fUseReco fill with candidates
+              
         if(fUseReco) {
            if (fCandidateType==kDstartoKpipi){
               new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
-              AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
+              AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass())); 
+              fhInvMassS->Fill(dstar->DeltaInvMass());
+
            } else{
               new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
               //Printf("Filling reco");
+              for(Int_t kd=0;kd<nprongs;kd++){
+                 fhImpPar->Fill(charmCand->Getd0Prong(kd),charmCand->PtProng(kd));
+              }
+              Double_t masses[2];
+              masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
+              masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
+              fhInvMassS->Fill(masses[0]);
+              fhInvMassS->Fill(masses[1]);
+              
            }               
            
            hstat->Fill(3);
+           
         }
         // for MC with requirement particle level fill with AliAODMCParticle
         else if (fUseMCInfo) {
@@ -396,9 +424,16 @@ void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
       else if(fUseReco){
         if (fCandidateType==kDstartoKpipi){
            new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
+           fhInvMassB->Fill(dstar->DeltaInvMass());
         }
         if (fCandidateType==kD0toKpi){
            new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
+           Double_t masses[2];
+           masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
+           masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
+           fhInvMassB->Fill(masses[0]);
+           fhInvMassB->Fill(masses[1]);
+
         }
         iSBCand++;
       }
@@ -473,6 +508,24 @@ void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
            hdeltaRDpi-> Fill(dRDpi,ptD);
            hdeltaRDK->  Fill(dRDK,ptD);
            
+           if (isMCBkg){
+              fhImpParB->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
+              fhImpParB->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
+              fhImpParB->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
+              
+           } else {
+              fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
+              fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
+              fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
+           
+           }
+           
+           
+        } else{
+           fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
+           AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
+           fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
+           fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
         }
 
       } //Dstar specific
@@ -681,13 +734,18 @@ Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
    //
    
    // Statistics 
-   TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
+   TH1I* hstat = new TH1I("hstat","Statistics",6,-0.5,5.5);
    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
-   if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
-   else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
+   if(fUseMCInfo){
+      if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N D");
+      else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
+   } else hstat->GetXaxis()->SetBinLabel(3, "N Cand");
    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");  
+   if(fUseMCInfo) {
+      hstat->GetXaxis()->SetBinLabel(6, "N Background");
+   }
    hstat->SetNdivisions(1);
    fOutput->Add(hstat);
    
@@ -716,6 +774,9 @@ Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
       fOutput->Add(hPtPion);
    }
    
+   fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
+   fOutput->Add(fhImpPar);
+   
    const Int_t nbinsalpha=200;
    Float_t minalpha=-TMath::Pi(), maxalpha=TMath::Pi();
    const Int_t nbinsdeltaR= 200;
@@ -795,7 +856,21 @@ Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
         fOutput->Add(hdeltaRDKB);
 
       }
-   
+      fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
+      fOutput->Add(fhImpParB);
+      
+      fhInvMassS = new TH1F("hInvMassS", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
+      fhInvMassS->SetStats(kTRUE);
+      fhInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
+      fhInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+      fOutput->Add(fhInvMassS);
+      
+      fhInvMassB = new TH1F("hInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
+      fhInvMassB->SetStats(kTRUE);
+      fhInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
+      fhInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+      fOutput->Add(fhInvMassB);
+
    }
    return kTRUE; 
 }
index adb0ced..30cdaa8 100644 (file)
@@ -97,8 +97,13 @@ class AliAnalysisTaskSEDmesonsFilterCJ : public AliAnalysisTaskSE
 
   TClonesArray *fCandidateArray;   //! contains candidates selected by AliRDHFCuts
   TClonesArray *fSideBandArray;    //! contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar case only!!)
-
-  ClassDef(AliAnalysisTaskSEDmesonsFilterCJ,2); // class for charm-jet correlations
+  //Histograms
+  TH2F* fhImpPar;
+  TH2F* fhImpParB;
+  TH1F* fhInvMassS;
+  TH1F* fhInvMassB;
+  
+  ClassDef(AliAnalysisTaskSEDmesonsFilterCJ,3); // class for charm-jet correlations
 };
 
 #endif
index a04e33c..d8f3ae6 100644 (file)
@@ -58,7 +58,7 @@ void *AddTaskDFilterAndCorrelations(
 
   TString candname="DStar"; 
   if(cand==0)  candname="D0";
-  TString sR = Form("R%.1f",R);
+  TString sR = Form("R%.0f",R*10);
   
   TString taskFiltername="DmesonsFilterCJ";
   taskFiltername+=candname;
index 3c8aabb..d999f73 100644 (file)
@@ -9,7 +9,7 @@ void AddTasksFlavourJet(const Int_t iCandType = 1 /*0 = D0, 1=Dstar...*/,
    const Bool_t bIsMC = kFALSE,
    const Bool_t bIsReco = kFALSE,
    const Bool_t bIsMap = kFALSE,
-   TString sText=""/*completes the name of the candidate task lists*/
+   TString sText="",/*completes the name of the candidate task lists*/
    const Bool_t bUseEMCalTrig = kFALSE,
    const Bool_t bIsTrain=kFALSE
    )
@@ -111,6 +111,7 @@ void AddTasksFlavourJet(const Int_t iCandType = 1 /*0 = D0, 1=Dstar...*/,
         bIsReco,
         "",
         taskFJ->GetName(),
+        sInputTrk,
         //Form("JetR%s",sRadius[i].Data()),
         iLeading,
         leadHadType,
@@ -132,6 +133,7 @@ void AddTasksFlavourJet(const Int_t iCandType = 1 /*0 = D0, 1=Dstar...*/,
                                      kFALSE,
                                      sText,
                                      taskFJ->GetName(),
+                                     sInputTrk,
                                      iLeading,
                                      leadHadType,
                                      aRadius[i],