Use Can define number of pt bins and their size; LS histos binned in pt (Renu)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 10 Dec 2009 22:28:07 +0000 (22:28 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 10 Dec 2009 22:28:07 +0000 (22:28 +0000)
PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
PWG3/vertexingHF/AliAnalysisTaskSEDplus.h

index 8acfa3a..de16c0a 100644 (file)
@@ -50,19 +50,15 @@ AliAnalysisTaskSE(),
 fOutput(0), 
 fHistNEvents(0),
 fNtupleDplus(0),
-fHistLS(0),
-fHistLStrip(0),
-fHistLStripcut(0),
-fHistOS(0),
-fHistOSbkg(0),
-fHistLSnoweight(0),
-fHistLSsinglecut(0),
+fUpmasslimit(1.965),
+fLowmasslimit(1.765),
+fNPtBins(0),
 fFillNtuple(kFALSE),
 fReadMC(kFALSE),
 fDoLS(kFALSE),
 fVHF(0)
 {
-  // Default constructor
+   // Default constructor
 }
 
 //________________________________________________________________________
@@ -71,21 +67,18 @@ AliAnalysisTaskSE(name),
 fOutput(0),
 fHistNEvents(0),
 fNtupleDplus(0),
-fHistLS(0),
-fHistLStrip(0),
-fHistLStripcut(0),
-fHistOS(0),
-fHistOSbkg(0),
-fHistLSnoweight(0),
-fHistLSsinglecut(0),
+fUpmasslimit(1.965),
+fLowmasslimit(1.765),
+fNPtBins(0),
 fFillNtuple(fillNtuple),
 fReadMC(kFALSE),
 fDoLS(kFALSE),
 fVHF(0)
 {
+  Double_t ptlim[5]={0.,2.,3.,5,9999999.};
+  SetPtBinLimit(5, ptlim);
   // Default constructor
-
-  // Output slot #1 writes into a TList container
+   // Output slot #1 writes into a TList container
   DefineOutput(1,TList::Class());  //My private output
 
   if(fFillNtuple){
@@ -102,16 +95,214 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     delete fOutput;
     fOutput = 0;
   }
-  if(fNtupleDplus){
-    delete fNtupleDplus;
-    fNtupleDplus=0;
-  }
   if (fVHF) {
     delete fVHF;
     fVHF = 0;
   }
-
+  
+  // if(fArrayBinLimits) {
+  //delete fArrayBinLimits;
+  //fArrayBinLimits= 0;
+  //} 
+  
 }  
+//_________________________________________________________________
+void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
+  fUpmasslimit = 1.865+range;
+  fLowmasslimit = 1.865-range;
+}
+//_________________________________________________________________
+void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
+  if(uplimit>lowlimit)
+    {
+      fUpmasslimit = lowlimit;
+      fLowmasslimit = uplimit;
+    }
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Double_t* lim){
+  // define pt bins for analysis
+  if(n>kMaxPtBins){
+    printf("Max. number of Pt bins = %d\n",kMaxPtBins);
+    fNPtBins=kMaxPtBins;
+    fArrayBinLimits[0]=0.;
+    fArrayBinLimits[1]=2.;
+    fArrayBinLimits[2]=3.;
+    fArrayBinLimits[3]=5.;
+    for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
+  }else{
+    fNPtBins=n-1;
+    fArrayBinLimits[0]=lim[0];
+    for(Int_t i=1; i<fNPtBins+1; i++) 
+      if(lim[i]>fArrayBinLimits[i-1]){
+       fArrayBinLimits[i]=lim[i];
+      }
+      else {
+       fArrayBinLimits[i]=fArrayBinLimits[i-1];
+      }
+    for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
+  }
+  if(fDebug > 1){
+    printf("Number of Pt bins = %d\n",fNPtBins);
+    for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);    
+  }
+}
+//_________________________________________________________________
+Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
+  if(ibin>fNPtBins)return -1;
+  return fArrayBinLimits[ibin];
+} 
+
+//_________________________________________________________________
+void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
+
+/*
+ * Fill the Like Sign histograms
+ */
+
+  Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
+
+  //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++;
+      }
+    }
+    if(track->Charge()<0)
+      {
+       nNegTrks++;
+       if(track->Pt()>=0.4){
+         nspcminus++;
+       }
+      }
+  }
+
+  Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
+
+  Int_t nDplusLS=0;
+  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
+  Int_t index; 
+
+  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 like 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())){
+
+      //set tight cuts values
+      Int_t iPtBin=-1;
+      Double_t ptCand = d->Pt();
+      if(ptCand<2.){
+       //iPtBin=0;
+       cutsDplus[7]=0.08;
+       cutsDplus[8]=0.5;
+       cutsDplus[9]=0.970;
+       cutsDplus[10]=0.0055;
+      }
+      else if(ptCand>2. && ptCand<3){ 
+       //iPtBin=1;
+       cutsDplus[7]=0.08;
+       cutsDplus[8]=0.5;
+       cutsDplus[9]=0.987;
+       cutsDplus[10]=0.005;
+      }else if(ptCand>3. && ptCand<5){ 
+       //iPtBin=2;
+       cutsDplus[7]=0.1;
+       cutsDplus[8]=0.5;
+       cutsDplus[9]=0.990;
+       cutsDplus[10]=0.0035;
+      }else{
+       //iPtBin=3;
+       cutsDplus[7]=0.1;
+       cutsDplus[8]=0.5;
+       cutsDplus[9]=0.995;
+       cutsDplus[10]=0.001;
+      }
+      
+      for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
+       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
+      }
+      
+      if(iPtBin<0){
+       return;
+      }
+
+      Bool_t passTightCuts=d->SelectDplus(cutsDplus);
+
+      Int_t sign= d->GetCharge();
+      Float_t wei=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
+       
+       wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
+       wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
+      }
+      
+      if(sign<0&&nNegTrks>2&&nspcminus>2){     
+       wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
+       wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
+
+      }
+
+      Float_t invMass = d->InvMassDplus();
+      
+      
+      index=GetLSHistoIndex(iPtBin);
+      fMassHistLS[index]->Fill(invMass,wei);
+      fMassHistLS[index+1]->Fill(invMass);
+      fMassHistLS[index+2]->Fill(invMass,wei2);
+      fMassHistLS[index+3]->Fill(invMass,wei3);
+      fMassHistLS[index+4]->Fill(invMass,wei4);
+
+      
+      if(passTightCuts){
+       fMassHistLSTC[index]->Fill(invMass,wei);
+       fMassHistLSTC[index+1]->Fill(invMass);
+       fMassHistLSTC[index+2]->Fill(invMass,wei2);
+       fMassHistLSTC[index+3]->Fill(invMass,wei3);
+       fMassHistLSTC[index+4]->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::Init()
@@ -140,72 +331,87 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
   fOutput->SetOwner();
   fOutput->SetName("OutputHistos");
 
-  Int_t nPtBins=4;
   TString hisname;
-  for(Int_t i=0;i<nPtBins;i++){
-    hisname.Form("hMassPt%d",i);
-    TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hm->Sumw2();
-    hisname.Form("hSigPt%d",i);
-    TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hs->Sumw2();
-    hisname.Form("hBkgPt%d",i);
-    TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hb->Sumw2();
+  Int_t index=0;
+  Int_t indexLS=0;
+  for(Int_t i=0;i<fNPtBins;i++){
 
+    index=GetHistoIndex(i);
+    indexLS=GetLSHistoIndex(i);
+    hisname.Form("hMassPt%d",i);
+    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHist[index]->Sumw2();
     hisname.Form("hMassPt%dTC",i);
-    TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hmtc->Sumw2();
+    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistTC[index]->Sumw2();
+    hisname.Form("hLSPt%dTC",i);
+    fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLSTC[indexLS]->Sumw2();
+    hisname.Form("hLSPt%dLC",i);
+    fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLS[indexLS]->Sumw2();
+    
+    index=GetSignalHistoIndex(i);    
+    indexLS++;
+    hisname.Form("hSigPt%d",i);
+    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHist[index]->Sumw2();
     hisname.Form("hSigPt%dTC",i);
-    TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hstc->Sumw2();
+    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistTC[index]->Sumw2();
+    hisname.Form("hLSPt%dLCnw",i);
+    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLS[indexLS]->Sumw2();
+    hisname.Form("hLSPt%dTCnw",i);
+    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLSTC[indexLS]->Sumw2();
+
+    index=GetBackgroundHistoIndex(i); 
+    indexLS++;
+    hisname.Form("hBkgPt%d",i);
+    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHist[index]->Sumw2();
     hisname.Form("hBkgPt%dTC",i);
-    TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hbtc->Sumw2();
-
-
-    fOutput->Add(hm);
-    fOutput->Add(hs);
-    fOutput->Add(hb);
-    fOutput->Add(hmtc);
-    fOutput->Add(hstc);
-    fOutput->Add(hbtc);
+    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistTC[index]->Sumw2();
+    hisname.Form("hLSPt%dLCntrip",i);
+    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLS[indexLS]->Sumw2();
+    hisname.Form("hLSPt%dTCntrip",i);
+    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLSTC[indexLS]->Sumw2();
+
+    indexLS++;
+    hisname.Form("hLSPt%dLCntripsinglecut",i);
+    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLS[indexLS]->Sumw2();
+    hisname.Form("hLSPt%dTCntripsinglecut",i);
+    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLSTC[indexLS]->Sumw2();
+
+    indexLS++;
+    hisname.Form("hLSPt%dLCspc",i);
+    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLS[indexLS]->Sumw2();
+    hisname.Form("hLSPt%dTCspc",i);
+    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
+    fMassHistLSTC[indexLS]->Sumw2();
   }
-  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);
+  for(Int_t i=0; i<3*fNPtBins; i++){
+    fOutput->Add(fMassHist[i]);
+    fOutput->Add(fMassHistTC[i]);
+  }
+  for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
+    fOutput->Add(fMassHistLS[i]);
+    fOutput->Add(fMassHistLSTC[i]);
+  }
 
-    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);
-  }
+  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
+ fHistNEvents->Sumw2();
+ fHistNEvents->SetMinimum(0);
+ fOutput->Add(fHistNEvents);
 
   if(fFillNtuple){
     OpenFile(2); // 2 is the slot number of the ntuple
@@ -221,14 +427,11 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   // Execute analysis for current event:
   // heavy flavor candidates association to MC truth
 
-  
-  
-  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
-  fHistNEvents->Fill(0); // count event
+   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
+ fHistNEvents->Fill(0); // count event
   // Post the data already here
   PostData(1,fOutput);
   
-
   TClonesArray *array3Prong = 0;
   TClonesArray *arrayLikeSign =0;
   if(!aod && AODEvent() && IsStandardAOD()) {
@@ -256,7 +459,7 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   }
   if(!arrayLikeSign) {
     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
-    // do in any case the OS analysis.
+    return;
   }
 
  
@@ -273,14 +476,14 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
     if(!arrayMC) {
       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
-      return;
+      //    return;
     }
     
   // load MC header
     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
     if(!mcHeader) {
-      printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
-      return;
+    printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
+    return;
     }
   }
   
@@ -289,6 +492,7 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   
   
   Int_t nOS=0;
+  Int_t index;
   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++) {
@@ -302,34 +506,40 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
     }
 
     if(d->SelectDplus(fVHF->GetDplusCuts())) {
-      Int_t iPtBin=0;
+      Int_t iPtBin = -1;
       Double_t ptCand = d->Pt();
+      
       if(ptCand<2.){
-       iPtBin=0;
+       //iPtBin=0;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;
-       cutsDplus[9]=0.979;
-       cutsDplus[10]=0.0000055;
+       cutsDplus[9]=0.970;
+       cutsDplus[10]=0.0055;
       }
       else if(ptCand>2. && ptCand<3){ 
-       iPtBin=1;
+       //iPtBin=1;
        cutsDplus[7]=0.08;
        cutsDplus[8]=0.5;
-       cutsDplus[9]=0.991;
-       cutsDplus[10]=0.000005;
+       cutsDplus[9]=0.987;
+       cutsDplus[10]=0.005;
       }else if(ptCand>3. && ptCand<5){ 
-       iPtBin=2;
+       //iPtBin=2;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
-       cutsDplus[9]=0.9955;
-       cutsDplus[10]=0.0000035;
+       cutsDplus[9]=0.990;
+       cutsDplus[10]=0.0035;
       }else{
-       iPtBin=3;
+       //iPtBin=3;
        cutsDplus[7]=0.1;
        cutsDplus[8]=0.5;
-       cutsDplus[9]=0.997;
-       cutsDplus[10]=0.000001;
+       cutsDplus[9]=0.995;
+       cutsDplus[10]=0.001;
+      }
+      
+      for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
+       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
       }
+      
       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
       Int_t labDp=-1;
       Float_t deltaPx=0.;
@@ -359,18 +569,6 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       }
       Double_t invMass=d->InvMassDplus();
 
-      TString hisNameA(Form("hMassPt%d",iPtBin));
-      TString hisNameS(Form("hSigPt%d",iPtBin));
-      TString hisNameB(Form("hBkgPt%d",iPtBin));
-      TString hisNameATC(Form("hMassPt%dTC",iPtBin));
-      TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
-      TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
-      
-      ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
-      if(passTightCuts){
-       ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
-      }
-
       Float_t tmp[22];
       if(fFillNtuple){           
        tmp[0]=pdgCode;
@@ -399,23 +597,41 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        PostData(2,fNtupleDplus);
       }
 
-      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);
+      if(iPtBin>=0){
+      
+       index=GetHistoIndex(iPtBin);
+       fMassHist[index]->Fill(invMass);
+       if(passTightCuts){
+         fMassHistTC[index]->Fill(invMass);
+
+       }
+       
+       if(fReadMC){
+         if(labDp>=0) {
+           index=GetSignalHistoIndex(iPtBin);
+           fMassHist[index]->Fill(invMass);
+
+           if(passTightCuts){
+             fMassHistTC[index]->Fill(invMass);
+
+           }
+           
+         }else{
+           index=GetBackgroundHistoIndex(iPtBin);
+           fMassHist[index]->Fill(invMass);
+
+           if(passTightCuts){
+             fMassHistTC[index]->Fill(invMass);
+
+           }   
          }
-         fHistOSbkg->Fill(invMass);
        }
       }
-
-      fHistOS->Fill(invMass);
+      /*
+      //start OS analysis
+      if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
+      fHistOS->Fill(d->InvMassDplus());
+      */
       nOS++;
     }
     if(unsetvtx) d->UnsetOwnPrimaryVtx();
@@ -423,129 +639,97 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
  
   //start LS analysis
   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
-   
+  
   PostData(1,fOutput);    
   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");
+
   fOutput = dynamic_cast<TList*> (GetOutputData(1));
   if (!fOutput) {     
     printf("ERROR: fOutput not available\n");
     return;
   }
-  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
+ fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
+
+ TString hisname;
+ Int_t index=0;
+ Int_t indexLS=0;
+ for(Int_t i=0;i<fNPtBins;i++){
+    index=GetHistoIndex(i);
+    if(fDoLS)indexLS=GetLSHistoIndex(i);
+    hisname.Form("hMassPt%d",i);
+    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    hisname.Form("hMassPt%dTC",i);
+    fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    if(fDoLS){
+      hisname.Form("hLSPt%dTC",i);
+      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      hisname.Form("hLSPt%dLC",i);
+      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    } 
+    
+    index=GetSignalHistoIndex(i);    
+    if(fDoLS)indexLS++;
+    hisname.Form("hSigPt%d",i);
+    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    hisname.Form("hSigPt%dTC",i);
+    fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    if(fDoLS){
+      hisname.Form("hLSPt%dLCnw",i);
+      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      hisname.Form("hLSPt%dTCnw",i);
+      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    }
+
+    index=GetBackgroundHistoIndex(i); 
+    if(fDoLS)indexLS++;
+    hisname.Form("hBkgPt%d",i);
+    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    hisname.Form("hBkgPt%dTC",i);
+    fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    if(fDoLS){
+      hisname.Form("hLSPt%dLCntrip",i);
+      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      hisname.Form("hLSPt%dTCntrip",i);
+      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    
+
+      indexLS++;
+      hisname.Form("hLSPt%dLCntripsinglecut",i);
+      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      
+      hisname.Form("hLSPt%dTCntripsinglecut",i);
+      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      
+      
+      indexLS++;
+      hisname.Form("hLSPt%dLCspc",i);
+      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+      
+      hisname.Form("hLSPt%dTCspc",i);
+      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
+    }
+  }
+
   if(fFillNtuple){
     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
   }
-  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;
index 66b65a3..28ab2b2 100644 (file)
@@ -16,6 +16,7 @@
 #include <TSystem.h>
 #include <TNtuple.h>
 #include <TH1F.h>
+#include <TArrayD.h>
 
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisVertexingHF.h"
@@ -27,34 +28,52 @@ class AliAnalysisTaskSEDplus : public AliAnalysisTaskSE
   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;}
+  void SetMassLimits(Float_t range);
+  void SetMassLimits(Float_t lowlimit, Float_t uplimit);
+  void SetPtBinLimit(Int_t n, Double_t *limitarray);
+  
+  Float_t GetUpperMassLimit(){return fUpmasslimit;}
+  Float_t GetLowerMassLimit(){return fLowmasslimit;}
+  Int_t GetNBinsPt(){return fNPtBins;}
+  Double_t GetPtBinLimit(Int_t ibin);
+  
+  void LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS);
+
   // 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); 
+  Int_t GetHistoIndex(Int_t iPtBin) const { return iPtBin*3;}
+  Int_t GetSignalHistoIndex(Int_t iPtBin) const { return iPtBin*3+1;}
+  Int_t GetBackgroundHistoIndex(Int_t iPtBin) const { return iPtBin*3+2;}
+  Int_t GetLSHistoIndex(Int_t iPtBin)const { return iPtBin*5;}
+  enum {kMaxPtBins=10};
+
   TList   *fOutput; //! list send on output slot 0
   TH1F    *fHistNEvents; //!hist. for No. of events
+  TH1F *fMassHist[3*kMaxPtBins]; //!hist. for inv mass (LC)
+  TH1F *fMassHistTC[3*kMaxPtBins]; //!hist. for inv mass (TC)
+  TH1F *fMassHistLS[5*kMaxPtBins];//!hist. for LS inv mass (LC)
+  TH1F *fMassHistLSTC[5*kMaxPtBins];//!hist. for LS inv mass (TC)
   TNtuple *fNtupleDplus; //! 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
-  
+  Float_t fUpmasslimit;  //upper inv mass limit for histos
+  Float_t fLowmasslimit; //lower inv mass limit for histos
+  Int_t fNPtBins; //number of bins in Pt for histograms
+  Double_t fArrayBinLimits[kMaxPtBins+1]; //limits for the Pt bins
   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
+  Bool_t fReadMC;    //flag for access to MC
+  Bool_t fDoLS;      //flag to do LS analysis
   AliAnalysisVertexingHF *fVHF;  // Vertexer heavy flavour (used to pass the cuts)
   
   ClassDef(AliAnalysisTaskSEDplus,4); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates