]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
Adding some histograms for jet shape
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
index 67ce106e39e832985665db2dcb9f24b1398df055..8cd6b4f1d8fa7badeab2674ca2d0710f0c1e0a42 100644 (file)
 ClassImp(AliAnalysisTaskJetSpectrum2)
 
 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
-  fJetHeaderRec(0x0),
-  fJetHeaderGen(0x0),
-  fAOD(0x0),
-  fhnCorrelation(0x0),
-  fBranchRec("jets"),
-  fBranchGen(""),
-  fUseAODJetInput(kFALSE),
-  fUseAODTrackInput(kFALSE),
-  fUseAODMCInput(kFALSE),
-  fUseGlobalSelection(kFALSE),
-  fUseExternalWeightOnly(kFALSE),
-  fLimitGenJetEta(kFALSE),
-  fFilterMask(0),
-  fAnalysisType(0),
+                                                           fJetHeaderRec(0x0),
+                                                           fJetHeaderGen(0x0),
+                                                           fAOD(0x0),
+                                                           fhnCorrelation(0x0),
+  fhnCorrelationPhiZRec(0x0),
+  f1PtScale(0x0),
+                                                           fBranchRec("jets"),
+                                                           fBranchGen(""),
+                                                           fUseAODJetInput(kFALSE),
+                                                           fUseAODTrackInput(kFALSE),
+                                                           fUseAODMCInput(kFALSE),
+                                                           fUseGlobalSelection(kFALSE),
+                                                           fUseExternalWeightOnly(kFALSE),
+                                                           fLimitGenJetEta(kFALSE),
+                                                           fFilterMask(0),
+                                                           fAnalysisType(0),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),
   fAvgTrials(1),
   fExternalWeight(1),    
-  fRecEtaWindow(0.5),
+                                                           fRecEtaWindow(0.5),
+  fMinJetPt(0),
+                                                           fDeltaPhiWindow(20./180.*TMath::Pi()),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -92,6 +96,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
   fh1PtTrackRec(0x0),   
   fh1SumPtTrackRec(0x0),  
   fh1SumPtTrackAreaRec(0x0),  
+  fh1TmpRho(0x0),
   fh1PtJetsRecIn(0x0),
   fh1PtJetsLeadingRecIn(0x0),
   fh1PtTracksRecIn(0x0),
@@ -104,13 +109,23 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
   fh2TracksLeadingPhiEta(0x0),
   fh2TracksLeadingPhiPt(0x0),
   fh2TracksLeadingJetPhiPt(0x0),
+  fh2JetPtJetPhi(0x0),
+  fh2TrackPtTrackPhi(0x0),
+  fh2DijetDeltaPhiPt(0x0),      
+  fh2DijetAsymPt(0x0),          
+  fh2DijetAsymPtCut(0x0),       
+  fh2DijetDeltaPhiDeltaEta(0x0),
+  fh2DijetPt2vsPt1(0x0),        
+  fh2DijetDifvsSum(0x0),        
+  fh1DijetMinv(0x0),            
+  fh1DijetMinvCut(0x0),         
   fHistList(0x0)  
 {
   for(int i = 0;i < kMaxStep*2;++i){
     fhnJetContainer[i] = 0;
   }
   for(int i = 0;i < kMaxJets;++i){
-    fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
+    fh2PhiPt[i] = fh2PhiEta[i] = fh2RhoPtRec[i] = fh2PsiPtRec[i] = fh2FragRec[i] = fh2PsiPtGen[i] = fh2FragGen[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
   }  
 
@@ -122,6 +137,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fJetHeaderGen(0x0),
   fAOD(0x0),
   fhnCorrelation(0x0),
+  fhnCorrelationPhiZRec(0x0),
+  f1PtScale(0x0),
   fBranchRec("jets"),
   fBranchGen(""),
   fUseAODJetInput(kFALSE),
@@ -137,6 +154,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fAvgTrials(1),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
+  fMinJetPt(0),
+  fDeltaPhiWindow(20./180.*TMath::Pi()),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -147,6 +166,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fh1PtTrackRec(0x0),   
   fh1SumPtTrackRec(0x0),  
   fh1SumPtTrackAreaRec(0x0),  
+  fh1TmpRho(0x0),
   fh1PtJetsRecIn(0x0),
   fh1PtJetsLeadingRecIn(0x0),
   fh1PtTracksRecIn(0x0),
@@ -159,6 +179,16 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
   fh2TracksLeadingPhiEta(0x0),
   fh2TracksLeadingPhiPt(0x0),
   fh2TracksLeadingJetPhiPt(0x0),
+  fh2JetPtJetPhi(0x0),
+  fh2TrackPtTrackPhi(0x0),
+  fh2DijetDeltaPhiPt(0x0),      
+  fh2DijetAsymPt(0x0),          
+  fh2DijetAsymPtCut(0x0),       
+  fh2DijetDeltaPhiDeltaEta(0x0),
+  fh2DijetPt2vsPt1(0x0),        
+  fh2DijetDifvsSum(0x0),        
+  fh1DijetMinv(0x0),            
+  fh1DijetMinvCut(0x0),         
   fHistList(0x0)
 {
 
@@ -166,7 +196,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
     fhnJetContainer[i] = 0;
   }  
   for(int i = 0;i < kMaxJets;++i){
-    fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
+    fh2PhiPt[i] = fh2PhiEta[i] = fh2RhoPtRec[i] = fh2PsiPtRec[i] = fh2FragRec[i] = fh2PsiPtGen[i] =  fh2RhoPtGen[i] = fh2FragGen[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
   }
   DefineOutput(1, TList::Class());  
@@ -230,7 +260,7 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
   //
   //  Histogram
     
-  const Int_t nBinPt = 100;
+  const Int_t nBinPt = 240;
   Double_t binLimitsPt[nBinPt+1];
   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
     if(iPt == 0){
@@ -253,6 +283,18 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
   }
 
 
+  const Int_t nBinPhi2 = 360;
+  Double_t binLimitsPhi2[nBinPhi2+1];
+  for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
+    if(iPhi2==0){
+      binLimitsPhi2[iPhi2] = 0.;
+    }
+    else{
+      binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
+    }
+  }
+
+
 
   const Int_t nBinEta = 40;
   Double_t binLimitsEta[nBinEta+1];
@@ -304,6 +346,10 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
+  fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
+  fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
+
+
   for(int ij = 0;ij < kMaxJets;++ij){
     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
@@ -314,6 +360,17 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
                              nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 
+    fh2RhoPtRec[ij] =  new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
+                               20,0.,1.,nBinPt,binLimitsPt);
+    fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
+                               20,0.,1.,nBinPt,binLimitsPt);
+
+    fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
+                               20,0.,1.,nBinPt,binLimitsPt);
+    fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
+                               20,0.,1.,nBinPt,binLimitsPt);
+    if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",20,0.,1);
+
 
     fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
                           nBinFrag,0.,1.,nBinPt,binLimitsPt);
@@ -321,11 +378,24 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
                             nBinFrag,0.,10.,nBinPt,binLimitsPt);
 
     fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
-                          nBinFrag,0.,1.,nBinPt,binLimitsPt);
+                             nBinFrag,0.,1.,nBinPt,binLimitsPt);
     fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
-                            nBinFrag,0.,10.,nBinPt,binLimitsPt);
+                               nBinFrag,0.,10.,nBinPt,binLimitsPt);
   }
 
+  // Dijet histograms
+
+  fh2DijetDeltaPhiPt       = new TH2F("fh2DeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
+  fh2DijetAsymPt            = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
+  fh2DijetAsymPtCut         = new TH2F("fh2DijetAsymCut","Pt asymmetry after delta phi cut;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
+  fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
+  fh2DijetPt2vsPt1          = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
+  fh2DijetDifvsSum          = new TH2F("fh2DijetDifvsSum","Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
+  fh1DijetMinv               = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
+  fh1DijetMinvCut           = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
+
+
+
 
   const Int_t saveLevel = 3; // large save level more histos
   if(saveLevel>0){
@@ -355,17 +425,35 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
     for(int ij = 0;ij<kMaxJets;++ij){
       fHistList->Add( fh1PtRecIn[ij]);
+
       if(fBranchGen.Length()>0){
        fHistList->Add( fh1PtGenIn[ij]);
        fHistList->Add( fh2FragGen[ij]);
        fHistList->Add( fh2FragLnGen[ij]);
+       fHistList->Add(fh2RhoPtGen[ij]);
+       fHistList->Add(fh2PsiPtGen[ij]);
+       fHistList->Add( fh2FragGen[ij]);
       }
       fHistList->Add( fh2PhiPt[ij]);
       fHistList->Add( fh2PhiEta[ij]);
+      fHistList->Add(fh2RhoPtRec[ij]);
+      fHistList->Add(fh2PsiPtRec[ij]);
       fHistList->Add( fh2FragRec[ij]);
       fHistList->Add( fh2FragLnRec[ij]);
     }
     fHistList->Add(fhnCorrelation);
+    fHistList->Add(fhnCorrelationPhiZRec);
+    fHistList->Add(fh2JetPtJetPhi);
+    fHistList->Add(fh2TrackPtTrackPhi);
+
+    fHistList->Add(fh2DijetDeltaPhiPt);       
+    fHistList->Add(fh2DijetAsymPt);       
+    fHistList->Add(fh2DijetAsymPtCut);               
+    fHistList->Add(fh2DijetDeltaPhiDeltaEta);        
+    fHistList->Add(fh2DijetPt2vsPt1);                
+    fHistList->Add(fh2DijetDifvsSum);                
+    fHistList->Add(fh1DijetMinv);                    
+    fHistList->Add(fh1DijetMinvCut);                 
   }
 
   // =========== Switch on Sumw2 for all histos ===========
@@ -387,7 +475,6 @@ void AliAnalysisTaskJetSpectrum2::Init()
   // Initialization
   //
 
-  Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
 
 }
@@ -500,19 +587,19 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
        pythiaGenHeader->TriggerJet(ip,p);
        pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
        
-       /*
-       if(fLimitGenJetEta){
-         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
-            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
-       }
-       */
-
        if(fBranchGen.Length()==0){
+         /*    
+         if(fLimitGenJetEta){
+           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
+              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
+         }
+         */
+         if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
          // if we have MC particles and we do not read from the aod branch
          // use the pythia jets
          genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
+         iCount++;
        }
-       iCount++;
       }
     }
     if(fBranchGen.Length()==0)nGenJets = iCount;    
@@ -553,6 +640,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
             tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
        }
        */
+       if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
        genJets[iCount] = *tmp;
        iCount++;
       }
@@ -584,7 +672,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
     Float_t tmpPt = tmp->Pt();
     if(tmpPt>ptOld){
-      Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
+      Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
     }
     ptOld = tmpPt;
   }
@@ -694,14 +782,17 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
     delete genIter;
   }
   
-  fh1NRecJets->Fill(nRecJets);
   nRecJets = TMath::Min(nRecJets,kMaxJets);
-
+  
+  Int_t iCountRec = 0;
   for(int ir = 0;ir < nRecJets;++ir){
     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
     if(!tmp)continue;
+    if(tmp->Pt()<fMinJetPt)continue;
     recJets[ir] = *tmp;
+    iCountRec++;
   }
+  nRecJets = iCountRec;
 
 
   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
@@ -729,6 +820,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 
 
   Double_t container[6];
+  Double_t containerPhiZ[6];
 
   // loop over generated jets
 
@@ -750,30 +842,38 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
     Int_t ir = iRecIndex[ig];
 
     if(TMath::Abs(etaGen)<fRecEtaWindow){
+      fh1TmpRho->Reset();
+
       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
       fh1PtGenIn[ig]->Fill(ptGen,eventW);
       // fill the fragmentation function
       for(int it = 0;it<genParticles.GetEntries();++it){
        AliVParticle *part = (AliVParticle*)genParticles.At(it);
-       if(genJets[ig].DeltaR(part)<radiusGen){
+       Float_t deltaR = genJets[ig].DeltaR(part);
+       fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
+       if(deltaR<radiusGen){
          Float_t z = part->Pt()/ptGen;
          Float_t lnz =  -1.*TMath::Log(z);
          fh2FragGen[ig]->Fill(z,ptGen,eventW);
          fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
        }
+
       }
-      if(ir>=0&&ir<nRecJets){
-       fhnJetContainer[kStep3]->Fill(&container[3],eventW);
+      Float_t rhoSum = 0;
+      for(int ibx = 1;ibx<fh2RhoPtGen[ir]->GetNbinsX();ibx++){
+       Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
+       Float_t rho = fh1TmpRho->GetBinContent(ibx);
+       rhoSum += rho;
+       fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
+       fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
       }
-    }    
-
+    }
     if(ir>=0&&ir<nRecJets){   
       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
       Double_t etaRec = recJets[ir].Eta();
       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
     }
   }// loop over generated jets
-
   
   Float_t sumPt = 0;
   for(int it = 0;it<recParticles.GetEntries();++it){
@@ -782,6 +882,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
     if(TMath::Abs(part->Eta())<0.9){
       Float_t pt = part->Pt();
       fh1PtTrackRec->Fill(pt,eventW);
+      fh2TrackPtTrackPhi->Fill(pt,part->Phi());
       sumPt += pt;
     }
   }
@@ -795,14 +896,56 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
     Double_t ptRec = recJets[ir].Pt();
     Double_t phiRec = recJets[ir].Phi();
     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
+
+
+  // do something with dijets...
+  if(ir==1){
+    Double_t etaRec1 = recJets[0].Eta();
+    Double_t ptRec1 = recJets[0].Pt();
+    Double_t phiRec1 = recJets[0].Phi();
+    if(phiRec1<0)phiRec1+=TMath::Pi()*2.;    
+    
+  
+    if(TMath::Abs(etaRec1)<fRecEtaWindow
+      &&TMath::Abs(etaRec)<fRecEtaWindow){
+  
+      Float_t deltaPhi = phiRec1 - phiRec;
+      
+      if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
+      if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();      
+      deltaPhi = TMath::Abs(deltaPhi);
+      fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);      
+      Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
+      fh2DijetAsymPt->Fill(asym,ptRec1);
+      fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
+      fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);        
+      fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);        
+      Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
+                        recJets[0].Px()*recJets[1].Px()- 
+                        recJets[0].Py()*recJets[1].Py()- 
+                        recJets[0].Pz()*recJets[1].Py());
+      minv = TMath::Sqrt(minv);
+      // with mass == 0;
+      
+      fh1DijetMinv->Fill(minv);            
+      if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
+       fh1DijetMinvCut->Fill(minv);         
+       fh2DijetAsymPtCut->Fill(asym,ptRec1);      
+      }
+    }
+  }
+
+
     container[0] = ptRec;
     container[1] = etaRec;
     container[2] = phiRec;
-
-    if(ptRec>10.&&fDebug>0){
+    containerPhiZ[0] = ptRec;
+    containerPhiZ[1] = phiRec;
+    if(ptRec>30.&&fDebug>0){
       // need to cast to int, otherwise the printf overwrites
       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
-      Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
+      Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
+      if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
       //  aodH->SetFillAOD(kTRUE);
       fAOD->GetHeader()->Print();
       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
@@ -822,10 +965,16 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
 
     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+    Float_t zLeading = -1;
     if(TMath::Abs(etaRec)<fRecEtaWindow){
+      fh2JetPtJetPhi->Fill(ptRec,phiRec);
       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
       fh1PtRecIn[ir]->Fill(ptRec,eventW);
       // fill the fragmentation function
+      
+      fh1TmpRho->Reset();
+
       for(int it = 0;it<recParticles.GetEntries();++it){
        AliVParticle *part = (AliVParticle*)recParticles.At(it);
        Float_t eta = part->Eta();
@@ -834,20 +983,38 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
          if(phi<0)phi+=TMath::Pi()*2.;    
          Float_t dPhi = phi - phiRec;
          Float_t dEta = eta - etaRec;
-         if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;    
+         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
+         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
          fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
          fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
        }
-       if(recJets[ir].DeltaR(part)<radiusRec){
+
+       Float_t deltaR = recJets[ir].DeltaR(part);
+       fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
+
+
+       if(deltaR<radiusRec){
          Float_t z = part->Pt()/ptRec;
+         if(z>zLeading)zLeading=z;
          Float_t lnz =  -1.*TMath::Log(z);
          fh2FragRec[ir]->Fill(z,ptRec,eventW);
          fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
        }
       }
+      // fill the jet shapes
+      Float_t rhoSum = 0;
+      for(int ibx = 1;ibx<fh2RhoPtRec[ir]->GetNbinsX();ibx++){
+       Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
+       Float_t rho = fh1TmpRho->GetBinContent(ibx);
+       rhoSum += rho;
+       fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
+       fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
+      }
     }
 
 
+    containerPhiZ[2] = zLeading;
+
     // Fill Correlation
     Int_t ig = iGenIndex[ir];
     if(ig>=0 && ig<nGenJets){
@@ -861,7 +1028,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
       container[3] = ptGen;
       container[4] = etaGen;
       container[5] = phiGen;
-
+      containerPhiZ[3] = ptGen;
       // 
       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
       // 
@@ -870,12 +1037,14 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
       if(TMath::Abs(etaRec)<fRecEtaWindow){
        fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
        fhnCorrelation->Fill(container);
+       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
+
       }// if etarec in window
 
     } 
-    ////////////////////////////////////////////////////  
     else{
-
+      containerPhiZ[3] = 0;
+      if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
     }
   }// loop over reconstructed jets
 
@@ -892,19 +1061,19 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
   // link it
   //
   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
-  const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
+  const Double_t kPtmin = 0.0, kPtmax = 160.; // we do not want to have empty bins at the beginning...
   const Double_t kEtamin = -3.0, kEtamax = 3.0;
   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
+  const Double_t kZmin = 0., kZmax = 1;
 
   // can we neglect migration in eta and phi?
   // phi should be no problem since we cover full phi and are phi symmetric
   // eta migration is more difficult  due to needed acceptance correction
   // in limited eta range
 
-
   //arrays for the number of bins in each dimension
   Int_t iBin[kNvar];
-  iBin[0] = 100; //bins in pt
+  iBin[0] = 160; //bins in pt
   iBin[1] =  1; //bins in eta 
   iBin[2] = 1; // bins in phi
 
@@ -942,14 +1111,42 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
   }
   fhnCorrelation->Sumw2();
 
+  // for second correlation histogram
+
+
+  const Int_t kNvarPhiZ   = 4; 
+  //arrays for the number of bins in each dimension
+  Int_t iBinPhiZ[kNvarPhiZ];
+  iBinPhiZ[0] = 80; //bins in pt
+  iBinPhiZ[1] = 72; //bins in phi 
+  iBinPhiZ[2] = 20; // bins in Z
+  iBinPhiZ[3] = 80; //bins in ptgen
+
+  //arrays for lower bounds :
+  Double_t* binEdgesPhiZ[kNvarPhiZ];
+  for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
+    binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
+
+  for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
+  for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
+  for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
+  for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
+
+  fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
+  for (int k=0; k<kNvarPhiZ; k++) {
+    fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
+  }
+  fhnCorrelationPhiZRec->Sumw2();
+
+
   // Add a histogram for Fake jets
-  //  thnDim[3] = AliPID::kSPECIES;
-  //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
-  // for(Int_t idim = 0; idim < kNvar; idim++)
-  //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
+
   for(Int_t ivar = 0; ivar < kNvar; ivar++)
     delete [] binEdges[ivar];
 
+  for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
+    delete [] binEdgesPhiZ[ivar];
+
 }
 
 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)