]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetSpectrum2.cxx
fixing a very nasty bug, and cleaning up
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetSpectrum2.cxx
index 67cdd07da77f72a4dca0c22a1fbeaf6f8b03d8bb..d5c624df7fbc3e444d6700606d8e8c225c1b6280 100644 (file)
@@ -61,7 +61,7 @@
 #include "AliJetKineReaderHeader.h"
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
-
+#include "AliCFContainer.h"
 
 #include "AliAnalysisHelperJetTasks.h"
 
@@ -100,6 +100,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
   fNTrigger(0),
   fTriggerBit(0x0),
   fNAcceptance(0),
+  fNBinsLeadingTrackPt(10),
+  fNBinsMult(20),
   fAnalysisType(0),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),
@@ -123,6 +125,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
   fTriggerName(0x0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
+  fh1AvgTrials(0x0),
   fh1PtHard(0x0),
   fh1PtHardNoW(0x0),  
   fh1PtHardTrials(0x0),
@@ -204,6 +207,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fNTrigger(0),
   fTriggerBit(0x0),
   fNAcceptance(0),
+  fNBinsLeadingTrackPt(10),
+  fNBinsMult(20),
   fAnalysisType(0),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),
@@ -227,6 +232,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fTriggerName(0x0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
+  fh1AvgTrials(0x0),
   fh1PtHard(0x0),
   fh1PtHardNoW(0x0),  
   fh1PtHardTrials(0x0),
@@ -326,6 +332,7 @@ Bool_t AliAnalysisTaskJetSpectrum2::Notify()
     // construct a poor man average trials 
     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
+    fh1Trials->Fill("#sum{ntrials}",ftrials); 
   }  
 
   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
@@ -373,14 +380,14 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
     
 
 
-  const Int_t nBinPt = 120;
+  const Int_t nBinPt = 150;
   Double_t binLimitsPt[nBinPt+1];
   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
     if(iPt == 0){
       binLimitsPt[iPt] = -50.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
     }
   }
   const Int_t nBinPhi = 90;
@@ -412,6 +419,9 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
   fHistList->Add(fh1Trials);
+  fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
+  fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
+  fHistList->Add(fh1AvgTrials);
   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
   fHistList->Add(fh1PtHard);
   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
@@ -503,18 +513,21 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
     // Bins:  Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
     const Int_t nBinsSparse1 = 9;
     const Int_t nBinsLeadingTrackPt = 10;
-    Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,120, 10,  20,    fNRPBins, 5,fNTrigger,nBinsLeadingTrackPt,fNAcceptance+1};
+    const Int_t nBinsArea = 10;
+    Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,120, 10,  fNBinsMult,    fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
     if(cJetBranch.Contains("RandomCone")){
       nBins1[1] = 600;
       nBins1[5] = 1;
     }
-    const Double_t xmin1[nBinsSparse1]  = {        -0.5,-50,  0,   0,        -0.5, 0.,-0.5,0.,-0.5,};
+    const Double_t xmin1[nBinsSparse1]  = {        -0.5,-50,  0,   0,        -0.5, 0.,         -0.5,  0.,           -0.5,};
     const Double_t xmax1[nBinsSparse1]  = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
     
+    const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.05,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
     const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {xmin1[7],1.,2.,3.,4.,5.,6.,8.,10.,12.,xmax1[7]}; //store pT of leading track in jet
 
     fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
-    fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
+    fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
+    if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
     fHistList->Add(fhnJetPt[ij]);
 
 
@@ -543,9 +556,10 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
     
     // Bins:track number  pTtrack, cent, mult, RP.   total bins = 224 k
     const Int_t nBinsSparse3 = 6;
-    const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    fNRPBins,fNTrigger};
+    const Int_t nRPBinsSparse3 = 3;
+    const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    nRPBinsSparse3,fNTrigger};
     const Double_t xmin3[nBinsSparse3]  = { -0.5,      0,   0,      0,        -0.5,-0.5};
-    const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,fNRPBins-0.5,fNTrigger-0.5};  
+    const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,nRPBinsSparse3-0.5,fNTrigger-0.5};  
     
       // change the binning ot the pT axis:
     Double_t *xPt3 = new Double_t[nBins3[1]+1];
@@ -803,7 +817,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
   Double_t eventW = 1;
   Double_t ptHard = 0; 
   Double_t nTrials = 1; // Trials for MC trigger 
-  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
+  fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); 
 
   // Getting some global properties
   fCentrality = GetCentrality();
@@ -1014,11 +1028,12 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles
       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
       //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
       Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
+      Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
 
       var1[1] = ptJet;
       var1[4] = phiBin;
       var1[5] = jet->EffectiveAreaCharged();
-      var1[7] = (leadTrack?leadTrack->Pt():0);//pT of leading jet
+      var1[7] = ptLead;
       var1[8] = CheckAcceptance(phiJet,etaJet);
 
       var2[1] = ptJet;
@@ -1026,10 +1041,10 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles
       var2[4] = phiJet;
       var2[5] = jet->EffectiveAreaCharged();
       var2[7] = var1[8];
-      var2[8] = (leadTrack?leadTrack->Charge()*leadTrack->Pt():0);//pT of leading jet
+      var2[8] = (leadTrack?leadTrack->Charge()*leadTrack->Pt():0);//pT of leading jet x charge
 
       if(ij<kMaxJets){
-       if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
+       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
        var1[0] = ij;
        var2[0] = ij;
        for(int it = 0;it <fNTrigger;it++){
@@ -1051,7 +1066,7 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles
        }
       }
 
-      if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
+      fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
 
       if(particlesList.GetSize()&&ij<kMaxJets){
        // Particles... correlated with jets...
@@ -1259,7 +1274,9 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ
     }
   }
 
-  Double_t container[6];
+  Double_t container[8];
+  Double_t containerRec[4];
+  Double_t containerGen[4];
 
   // loop over generated jets
   // consider the 
@@ -1269,18 +1286,29 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ
     Double_t phiGen = genJet->Phi();
     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
     Double_t etaGen = genJet->Eta();
-    container[3] = ptGen;
-    container[4] = etaGen;
-    container[5] = phiGen;
-    fhnJetContainer->Fill(&container[3],kStep0);
-    if(JetSelected(genJet)){
-      fhnJetContainer->Fill(&container[3],kStep1);
-      Int_t ir = aRecIndex[ig];
-      if(ir>=0&&ir<recJetsList.GetEntries()){   
-       fhnJetContainer->Fill(&container[3],kStep2);
+    containerGen[0] = ptGen;
+    containerGen[1] = etaGen;
+    containerGen[2] = phiGen;
+    containerGen[3] = genJet->GetPtLeading();
+
+    fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
+
+    Int_t ir = aRecIndex[ig];
+    if(ir>=0&&ir<recJetsList.GetEntries()){   
+      fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
+      
+      if(JetSelected(genJet)){
+       fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
+
        AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir); 
-       if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
-       if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
+
+       fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
+       if(JetSelected(recJet)) {
+         fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
+
+         containerGen[3] = recJet->GetPtLeading();
+         fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner in eta window with leading track on reconstructed level
+       }
       }
     }
   }// loop over generated jets used for matching...
@@ -1294,37 +1322,44 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ
     Double_t ptRec = recJet->Pt();
     Double_t phiRec = recJet->Phi();
     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
-    // do something with dijets...
     
+    containerRec[0] = ptRec;
+    containerRec[1] = etaRec;
+    containerRec[2] = phiRec;
+    containerRec[3] = recJet->GetPtLeading();
+
     container[0] = ptRec;
     container[1] = etaRec;
     container[2] = phiRec;
+    container[3] = recJet->GetPtLeading();
 
-    fhnJetContainer->Fill(container,kStep0+kMaxStep);
     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
   
     if(JetSelected(recJet)){
-      fhnJetContainer->Fill(container,kStep1+kMaxStep);
+      fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
       // Fill Correlation
       Int_t ig = aGenIndex[ir];
       if(ig>=0 && ig<genJetsList.GetEntries()){
-       fhnJetContainer->Fill(container,kStep2+kMaxStep);
        if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
        AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
        Double_t ptGen  = genJet->Pt();
        Double_t phiGen = genJet->Phi();
        if(phiGen<0)phiGen+=TMath::Pi()*2.; 
        Double_t etaGen = genJet->Eta();
-      
-       container[3] = ptGen;
-       container[4] = etaGen;
-       container[5] = phiGen;
-       // 
-       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
-       // 
-       if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
-       fhnJetContainer->Fill(container,kStep3+kMaxStep);
-       fhnCorrelation->Fill(container,0);
+
+       containerGen[0] = ptGen;
+       containerGen[1] = etaGen;
+       containerGen[2] = phiGen;
+       containerGen[3] = genJet->GetPtLeading();
+
+       container[4] = ptGen;
+       container[5] = etaGen;
+       container[6] = phiGen;
+       container[7] = genJet->GetPtLeading();
+
+       fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
+
+       fhnCorrelation->Fill(container);
        if(ptGen>0){
          Float_t delta = (ptRec-ptGen)/ptGen;
          fh2RelPtFGen->Fill(ptGen,delta);
@@ -1342,10 +1377,11 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
   // Create the particle container for the correction framework manager and 
   // link it
   //
-  const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
+  const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
-  const Double_t kEtamin = -3.0, kEtamax = 3.0;
+  const Double_t kEtamin = -1.0, kEtamax = 1.0;
   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
+  const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
 
   // can we neglect migration in eta and phi?
   // phi should be no problem since we cover full phi and are phi symmetric
@@ -1355,8 +1391,9 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
   //arrays for the number of bins in each dimension
   Int_t iBin[kNvar];
   iBin[0] = 125; //bins in pt
-  iBin[1] =  1; //bins in eta 
-  iBin[2] = 1; // bins in phi
+  iBin[1] =  4; //bins in eta 
+  iBin[2] =  1; // bins in phi
+  iBin[3] =  10; // bins in leading track Pt
 
   //arrays for lower bounds :
   Double_t* binEdges[kNvar];
@@ -1368,9 +1405,19 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
-
-
-  fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin);
+  binEdges[3][0]= kPtLeadingTrackPtMin;
+  binEdges[3][1]= 1.;
+  binEdges[3][2]= 2.;
+  binEdges[3][3]= 3.;
+  binEdges[3][4]= 4.;
+  binEdges[3][5]= 5.;
+  binEdges[3][6]= 6.;
+  binEdges[3][7]= 8.;
+  binEdges[3][8]= 10.;
+  binEdges[3][9]= 12.;
+  binEdges[3][10]= kPtLeadingTrackPtMax;
+
+  fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
   for (int k=0; k<kNvar; k++) {
     fhnJetContainer->SetBinLimits(k,binEdges[k]);
   }
@@ -1384,10 +1431,10 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
     thnDim[k+kNvar] = iBin[k];
   }
 
-  fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
+  fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
   for (int k=0; k<kNvar; k++) {
-    fhnCorrelation->SetBinLimits(k,binEdges[k]);
-    fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
+    fhnCorrelation->SetBinEdges(k,binEdges[k]);
+    fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
   }
 
   for(Int_t ivar = 0; ivar < kNvar; ivar++)