]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCore.cxx
Add classimp
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
index 6f7b635efad146837b4142de72127ea8b74284ed..66e1c556bc23e69f1e83af466e06342f86990e16 100644 (file)
@@ -1,4 +1,3 @@
-
 // ******************************************
 // This task computes several jet observables like 
 // the fraction of energy in inner and outer coronnas,
@@ -32,7 +31,7 @@
 #include "TH3F.h"
 #include "THnSparse.h"
 #include "TCanvas.h"
-
+#include "TRandom3.h"
 #include "AliLog.h"
 
 #include "AliAnalysisTask.h"
@@ -45,6 +44,7 @@
 #include "AliAnalysisHelperJetTasks.h"
 #include "AliInputEventHandler.h"
 #include "AliAODJetEventBackground.h"
+#include "AliAODMCParticle.h"
 #include "AliAnalysisTaskFastEmbedding.h"
 #include "AliAODEvent.h"
 #include "AliAODHandler.h"
@@ -52,6 +52,9 @@
 
 #include "AliAnalysisTaskJetCore.h"
 
+using std::cout;
+using std::endl;
+
 ClassImp(AliAnalysisTaskJetCore)
 
 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
@@ -70,19 +73,37 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
 fCentMax(100.),
 fNInputTracksMin(0),
 fNInputTracksMax(-1),
+fRequireITSRefit(0),
+fApplySharedClusterCut(0),
 fAngStructCloseTracks(0),
 fCheckMethods(0),
 fDoEventMixing(0), 
 fFlagPhiBkg(0),
+fFlagEtaBkg(0),
+fFlagJetHadron(0),
+fFrac(0.8),
+fTTLowRef(11),
+fTTUpRef(13),
+fTTLowSig(15),
+fTTUpSig(19),
+fHardest(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fFlagOnlyHardest(1),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(50),
+fSemigoodCorrect(0),
+fHolePos(4.71),
+fHoleWidth(0.2),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -95,6 +116,7 @@ fJetPtFractionMin(0.5),
 fNMatchJets(4),
 fMatchMaxDist(0.8),
 fKeepJets(kFALSE),
+fRunAnaAzimuthalCorrelation(kFALSE),
 fkNbranches(2),
 fkEvtClasses(12),
 fOutputList(0x0),
@@ -112,7 +134,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -129,14 +150,21 @@ fh2AngStructpt1C60(0x0),
 fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),
+fh1TrigRef(0x0),
+fh1TrigSig(0x0),
 fh2Ntriggers(0x0),
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
-fh3spectriggeredC20(0x0),
-fh3spectriggeredC3060(0x0)
-
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0), 
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fHJetSpec(0x0),
+fhTTPt(0x0),
+fHJetPhiCorr(0x0),
+fRandom(0x0)
  
 {
    // default Constructor
@@ -176,19 +204,37 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
 fCentMax(100.),
 fNInputTracksMin(0),
 fNInputTracksMax(-1),
+fRequireITSRefit(0),
+fApplySharedClusterCut(0),
 fAngStructCloseTracks(0),
 fCheckMethods(0),
 fDoEventMixing(0),
 fFlagPhiBkg(0),
+fFlagEtaBkg(0),
+fFlagJetHadron(0),
+fFrac(0.8),
+fTTLowRef(11),
+fTTUpRef(13),
+fTTLowSig(15),
+fTTUpSig(19),
+fHardest(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fFlagOnlyHardest(1),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(50),
+fSemigoodCorrect(0),
+fHolePos(4.71),
+fHoleWidth(0.2),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -201,6 +247,7 @@ fJetPtFractionMin(0.5),
 fNMatchJets(4),
 fMatchMaxDist(0.8),
 fKeepJets(kFALSE),
+fRunAnaAzimuthalCorrelation(kFALSE),
 fkNbranches(2),
 fkEvtClasses(12),
 fOutputList(0x0),
@@ -218,7 +265,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -235,13 +281,21 @@ fh2AngStructpt1C60(0x0),
 fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),    
+fh1TrigRef(0x0),
+fh1TrigSig(0x0),
 fh2Ntriggers(0x0),
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
-fh3spectriggeredC20(0x0),
-fh3spectriggeredC3060(0x0)
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0),
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fHJetSpec(0x0),
+fhTTPt(0x0),
+fHJetPhiCorr(0x0),
+fRandom(0x0)
 
  {
    // Constructor
@@ -298,6 +352,12 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
    TH1::AddDirectory(kFALSE);
 
 
+        // set seed
+   fRandom = new TRandom3(0);
+  
+
+
+
    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
@@ -344,7 +404,6 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
 
     if(fCheckMethods){
-
     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
@@ -354,9 +413,8 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
 
-    fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,34,0,3.4);
-    fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,34,0,3.4);
-    fh3JetTrackC4080=new TH3F("JetTrackC4080","",50,0,50,150,0.,150.,34,0,3.4);
+    fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,35,0.,3.5);
+    fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,35,0.,3.5);
     if(fAngStructCloseTracks>0){
     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
@@ -376,14 +434,18 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
 
    
+    fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
+    fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);  
+    fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
+    fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
+    fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
+    fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.07","",100,0.,4000.,100,0.,5.,25,0.,50.);
+    fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,25,0.,50.);
+    fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
+    fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.); 
+    fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.); 
+    fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);  
 
-    fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
-    fh2JetDensity=new TH2F("Jet density vs centrality A>0.4","",100,0.,4000.,100,0.,5.);
-    fh2JetDensityA4=new TH2F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.);
-    fh2RPJets=new TH2F("RPJet","",3,0.,3.,150,0.,150.); 
-    fh3spectriggeredC4080 = new TH3F("Triggered spectrumC4080","",5,0.,1.,140,-80.,200.,50,0.,50.);
-    fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",5,0.,1.,140,-80.,200.,50,0.,50.);
-    fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",5,0.,1.,140,-80.,200.,50,0.,50.);
 
     
     
@@ -408,7 +470,7 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
       
       fOutputList->Add(fh3JetTrackC3060);
       fOutputList->Add(fh3JetTrackC20);
-      fOutputList->Add(fh3JetTrackC4080); 
+     
             
      
 
@@ -434,15 +496,58 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 
 
 
+        fOutputList->Add(fh1TrigRef);
+        fOutputList->Add(fh1TrigSig); 
        fOutputList->Add(fh2Ntriggers);
-        fOutputList->Add(fh2JetDensity);
-        fOutputList->Add(fh2JetDensityA4);
-        fOutputList->Add(fh2RPJets);
-       fOutputList->Add(fh3spectriggeredC4080);
-       fOutputList->Add(fh3spectriggeredC20); 
-       fOutputList->Add(fh3spectriggeredC3060);   
+        fOutputList->Add(fh2Ntriggers2C10);
+        fOutputList->Add(fh2Ntriggers2C20); 
+        fOutputList->Add(fh3JetDensity);
+        fOutputList->Add(fh3JetDensityA4);
+        fOutputList->Add(fh2RPJetsC10);
+        fOutputList->Add(fh2RPJetsC20);
+         fOutputList->Add(fh2RPTC10);
+        fOutputList->Add(fh2RPTC20);
+
+        const Int_t dimSpec = 5;
+       const Int_t nBinsSpec[dimSpec]     = {100,6, 140, 50, fNRPBins};
+       const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
+       const Double_t hiBinSpec[dimSpec]  = {100,1, 200, 50,  static_cast<Double_t>(fNRPBins)};
+       fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
+
+             //change binning in jet area
+     Double_t *xPt6 = new Double_t[7];
+     xPt6[0] = 0.;
+     xPt6[1]=0.07;
+     xPt6[2]=0.2;
+     xPt6[3]=0.4;
+     xPt6[4]=0.6;
+     xPt6[5]=0.8; 
+     xPt6[6]=1;
+    fHJetSpec->SetBinEdges(1,xPt6);
+    delete [] xPt6;
+
+
+
+
+
+       fOutputList->Add(fHJetSpec);  
+
+
+       if(fRunAnaAzimuthalCorrelation)
+         {
+           fhTTPt = new TH2F("fhTTPt","Trigger track p_{T} vs centrality",10,0,100,100,0,100);
+           fOutputList->Add(fhTTPt);
+
+           const Int_t dimCor = 5;
+           const Int_t nBinsCor[dimCor]     = {50, 200, 100,              8,   10};
+           const Double_t lowBinCor[dimCor] = {0,  -50, -0.5*TMath::Pi(), 0,   0};
+           const Double_t hiBinCor[dimCor]  = {50, 150, 1.5*TMath::Pi(),  0.8, 100};
+           fHJetPhiCorr = new THnSparseF("fHJetPhiCorr","TT p_{T} vs jet p_{T} vs dPhi vs area vs centrality",dimCor,nBinsCor,lowBinCor,hiBinCor);
+           fOutputList->Add(fHJetPhiCorr);
+         }
+
 
+        
      
    // =========== Switch on Sumw2 for all histos ===========
    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
@@ -494,16 +599,23 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
     }}
     
 
+   // -- event selection --
+   fHistEvtSelection->Fill(1); // number of events before event selection
 
 
+       Bool_t selected=kTRUE;
+       selected = AliAnalysisHelperJetTasks::Selected();
+       if(!selected){
+      // no selection by the service task, we continue
+       PostData(1,fOutputList);
+       return;}
+    
 
-   // -- event selection --
-   fHistEvtSelection->Fill(1); // number of events before event selection
 
-   // physics selection
+  // physics selection: this is now redundant, all should appear as accepted after service task selection
    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
-   cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
+        std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
       fHistEvtSelection->Fill(2);
@@ -511,6 +623,8 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
       return;
    }
 
+
+      
    // vertex selection
    if(!aod){
      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
@@ -585,9 +699,10 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
    if(fFlagRandom==0){
      if(externalBackground)rho = externalBackground->GetBackground(0);}
-   if(fFlagRandom==1){
+   if(fFlagRandom==2){
       if(externalBackground)rho = externalBackground->GetBackground(2);}
-
+   if(fFlagRandom==3){
+      if(externalBackground)rho = externalBackground->GetBackground(3);}
    // fetch jets
    TClonesArray *aodJets[2];
    aodJets[0]=0;
@@ -602,21 +717,33 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  } 
 
 
-   //Double_t ptsub[aodJets[0]->GetEntriesFast()];
-   //Int_t inord[aodJets[0]->GetEntriesFast()];
-   //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
-   //  ptsub[n]=0;
-   //  inord[n]=0;}   
 
+   Int_t nT=0;
    TList ParticleList;
-   Int_t nT = GetListOfTracks(&ParticleList);
+   Double_t minT=0;
+   Double_t maxT=0;
+   Int_t number=0;
+   Double_t dice=fRandom->Uniform(0,1);
+   if(dice>fFrac){ minT=fTTLowRef;
+                   maxT=fTTUpRef;}
+   if(dice<=fFrac){minT=fTTLowSig;
+                   maxT=fTTUpSig;} 
+
+
+
+   if(fHardest==1 || fHardest==2) nT = GetListOfTracks(&ParticleList);
+   if(fHardest==0) nT=SelectTrigger(&ParticleList,minT,maxT,number);
+   if(nT<0){  
+   PostData(1, fOutputList);
+   return;}   
+
+      if(dice>fFrac) fh1TrigRef->Fill(number);
+      if(dice<=fFrac)fh1TrigSig->Fill(number);
+
      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
       fListJets[iJetType]->Clear();
       if (!aodJets[iJetType]) continue;
-
       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
-      
-   
       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
          if (jet) fListJets[iJetType]->Add(jet);
@@ -630,7 +757,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Double_t phibig=0.;
    Double_t etasmall=0;
    Double_t ptsmall=0;
-   Double_t areasmall=0;
+   //   Double_t areasmall=0;
    Double_t phismall=0.;
          
   
@@ -638,80 +765,148 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Int_t iCount=0; 
    Int_t trigJet=-1;
    Int_t trigBBTrack=-1;
-   Int_t trigInTrack=-1;
+   //   Int_t trigInTrack=-1;
    fRPAngle = aod->GetHeader()->GetEventplane();     
 
+   if(fHardest==0 || fHardest==1){
    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
    if(!partback){  
    PostData(1, fOutputList);
    return;}
-   fh2Ntriggers->Fill(centValue,partback->Pt());
+
+    if(fSemigoodCorrect){
+   Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
+   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){ 
+   PostData(1, fOutputList);
+   return;}} 
+
+
+    }
+
+
+   for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
+     if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
+     AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
+     if(!partback) continue;
+     if(partback->Pt()<8) continue;
+
    Double_t accep=2.*TMath::Pi()*1.8;
    Int_t injet4=0;
    Int_t injet=0; 
+
+   if(fSemigoodCorrect){
+   Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
+   if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6) continue;}
+
+   
+
+   fh2Ntriggers->Fill(centValue,partback->Pt());
+   Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
+   if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
+   if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
+
+
+
    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
            etabig  = jetbig->Eta();
            phibig  = jetbig->Phi();
            ptbig   = jetbig->Pt();
            if(ptbig==0) continue; 
-           Int_t phiBin = GetPhiBin(phibig-fRPAngle);       
+           Double_t phiBin = RelativePhi(phibig,fRPAngle);       
            areabig = jetbig->EffectiveAreaCharged();
            Double_t ptcorr=ptbig-rho*areabig;
           if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
-       
-           if(areabig>=0.2) injet=injet+1;
+           if(areabig>=0.07) injet=injet+1;
            if(areabig>=0.4) injet4=injet4+1;   
-           Double_t dphi=RelativePhi(partback->Phi(),phibig);  
-           if(centValue>40. && centValue<80.) fh3JetTrackC4080->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
+           Double_t dphi=RelativePhi(partback->Phi(),phibig); 
+
+           if(fFlagEtaBkg==1){
+          Double_t etadif= partback->Eta()-etabig;
+           if(TMath::Abs(etadif)<=0.5){             
+          
+           if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
+           if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
+           if(fFlagEtaBkg==0){
            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
-           if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
+           if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
 
-           if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
+
+           if(fFlagJetHadron==0){
+           if(fFlagPhiBkg==1) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
            if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
-                   if(centValue<20.) fh2RPJets->Fill(phiBin, ptcorr);
+           if(fFlagPhiBkg==2) if(TMath::Abs(dphi)<TMath::Pi()-0.7) continue;
+           if(fFlagPhiBkg==3) if(TMath::Abs(dphi)<TMath::Pi()-0.5) continue;}
+           if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
+
+
+          if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
+          if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
                    Double_t dismin=100.;
                    Double_t ptmax=-10.; 
                    Int_t index1=-1;
                    Int_t index2=-1;
-           if(centValue>40. && centValue<80.)  fh3spectriggeredC4080->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
-           if(centValue<20.)  fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
-           if(centValue>30. && centValue<60.)  fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
+         
+                  Float_t phitt=partback->Phi();
+                   if(phitt<0)phitt+=TMath::Pi()*2.; 
+                   Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
+
+                  Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(), static_cast<Double_t>(phiBintt)};
+                 fHJetSpec->Fill(fillspec);
+           
+          
 
                    if(ptcorr<=0) continue;
-                       AliAODTrack* leadtrack; 
+
+                       AliAODTrack* leadtrack=0; 
                        Int_t ippt=0;
-                       Double_t ppt=-10;   
-                      TRefArray *genTrackList = jetbig->GetRefTracks();
-                       Int_t nTracksGenJet = genTrackList->GetEntriesFast();
-                       AliAODTrack* genTrack;
-                       for(Int_t ir=0; ir<nTracksGenJet; ++ir){
-                       genTrack = (AliAODTrack*)(genTrackList->At(ir));
-                      if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
-                        ippt=ir;}}
-                        leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
-                        if(!leadtrack) continue;
-                       
-                       //store one trigger info                   
-                        if(iCount==0){                        
-                         trigJet=i;
-                          trigBBTrack=nT;
-                          trigInTrack=ippt;
-                          iCount=iCount+1;} 
+                       Double_t ppt=-10;
+                       if(fFlagJetHadron==0){   
+                        TRefArray *genTrackList = jetbig->GetRefTracks();
+                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
+                        AliAODTrack* genTrack;
+                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
+                          genTrack = (AliAODTrack*)(genTrackList->At(ir));
+                          if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
+                            ippt=ir;}}
+                        leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
+                        if(!leadtrack) continue;
+                      }
 
+
+
+                      AliVParticle* leadtrackb=0;
+                       if(fFlagJetHadron!=0){
+                        Int_t nTb = GetHardestTrackBackToJet(jetbig);
+                         leadtrackb = (AliVParticle*)ParticleList.At(nTb);
+                         if(!leadtrackb) continue;  
+                      }
+
+
+
+
+                       
+                      //store one trigger info                   
+                      if(iCount==0){                        
+                        trigJet=i;
+                        trigBBTrack=nT;
+                        //                      trigInTrack=ippt;
+                        iCount=iCount+1;} 
+                      
    
-                 if(fCheckMethods){
-                  for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
-                  AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
-                  etasmall  = jetsmall->Eta();
-                  phismall = jetsmall->Phi();
-                  ptsmall   = jetsmall->Pt();
-                  areasmall = jetsmall->EffectiveAreaCharged();
-                  Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
-                 tmpDeltaR=TMath::Sqrt(tmpDeltaR);
-                    //Fraction in the jet core  
-                    if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
-                   index2=j;}  
+                      if(fCheckMethods){
+                        for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
+                          AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
+                          etasmall  = jetsmall->Eta();
+                          phismall = jetsmall->Phi();
+                          ptsmall   = jetsmall->Pt();
+                          //                      areasmall = jetsmall->EffectiveAreaCharged();
+                          Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
+                          tmpDeltaR=TMath::Sqrt(tmpDeltaR);
+                          //Fraction in the jet core  
+                          if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
+                            index2=j;}  
                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
                      index1=j;}} //en of loop over R=0.2 jets
                 //method1:most concentric jet=core 
@@ -727,21 +922,27 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
                  if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
                  if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
-       
-                  
-         if(fDoEventMixing==0){
+                 if(centValue<10&&leadtrack) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());                  
+                  if(centValue<20&&leadtrack) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());  
+         if(fDoEventMixing==0 && fFlagOnlyRecoil==0){ 
         for(int it = 0;it<ParticleList.GetEntries();++it){
          AliVParticle *part = (AliVParticle*)ParticleList.At(it);
                  Double_t deltaR = jetbig->DeltaR(part);
           Double_t deltaEta = etabig-part->Eta();
          
           Double_t deltaPhi=phibig-part->Phi();
           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
-                 Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()};                     fhnDeltaR->Fill(jetEntries);}
+         Double_t pTcont=0;
+          if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
+          if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt(); 
+          Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};  
+           fhnDeltaR->Fill(jetEntries);}
+
+
+         }
+        
         
-        }
           //end of track loop, we only do it if EM is switched off
          
 
@@ -753,11 +954,11 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
 
    }
-   if(injet>0) fh2JetDensity->Fill(ParticleList.GetEntries(),injet/accep);
-   if(injet4>0)fh2JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep);
+   if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
+   if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
           //end of jet loop
 
-
+   //}
 
 
           if(fDoEventMixing>0){
@@ -786,6 +987,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                         fNevents=fNevents+1;  
                         if(fNevents==10) fTindex=fTindex+1; 
            }}}
+
               if(fTindex==10&&fNevents==10) fCountAgain=0;
 
                // Copy the triggers from the current event into the buffer.
@@ -807,6 +1009,33 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
          
          }
 
+         /////////////////////////////////////////////////////////////////////////////
+         ////////////////////// Rongrong's analysis //////////////////////////////////
+         if(fRunAnaAzimuthalCorrelation)
+           {
+             fhTTPt->Fill(centValue,partback->Pt());
+             for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
+               {
+                 AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
+                 Double_t jetPt   = jet->Pt();
+                 Double_t jetEta  = jet->Eta();
+                 Double_t jetPhi  = jet->Phi();
+                 if(jetPt==0) continue; 
+                 if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
+                 Double_t jetArea = jet->EffectiveAreaCharged();
+                 Double_t jetPtCorr=jetPt-rho*jetArea;
+                 Double_t dPhi=jetPhi-partback->Phi();
+                 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
+                 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
+                 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
+                 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
+
+                 Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
+                 fHJetPhiCorr->Fill(fill);
+               }
+           }
+         /////////////////////////////////////////////////////////////////////////////
+         /////////////////////////////////////////////////////////////////////////////
 
 
      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
@@ -900,7 +1129,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
 
     
 
-
+   }
 
 
 
@@ -923,29 +1152,135 @@ void AliAnalysisTaskJetCore::Terminate(const Option_t *)
 
 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
 
+      Int_t iCount = 0;
+      AliAODEvent *aod = 0;
+
+     if(!fESD)aod = fAODIn;
+     else aod = fAODOut;   
+
+     if(!aod)return 0;
+
+     Int_t index=-1;
+     Double_t ptmax=-10;
+
+
+    
+     for(int it = 0;it < aod->GetNumberOfTracks();++it){
+      AliAODTrack *tr = aod->GetTrack(it);
+      Bool_t bGood = false;
+      if(fFilterType == 0)bGood = true;
+      else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
+      else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();    
+      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+      if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
+      if(bGood==false) continue;
+      if (fApplySharedClusterCut) {
+           Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
+           if (frac > 0.4) continue;
+      }
+     if(TMath::Abs(tr->Eta())>0.9)continue;
+      if(tr->Pt()<0.15)continue;
+      list->Add(tr);
+      iCount++;
+      if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
+       if(tr->TestFilterBit(fFilterMaskBestPt)){
+         if(tr->Pt()>ptmax){ 
+           ptmax=tr->Pt();     
+           index=iCount-1;
+         }
+       }
+      }
+      else{
+       if(tr->Pt()>ptmax){ 
+         ptmax=tr->Pt();       
+         index=iCount-1;
+       }
+      }
+     }
+  
+   
+    // else if (type == kTrackAODMCCharged) {
+    // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
+    // if(!tca)return iCount;
+    // for(int it = 0;it < tca->GetEntriesFast();++it){
+    //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
+    //   if(!part)continue;
+    //   if(part->Pt()<0.15)continue;
+    //   if(!part->IsPhysicalPrimary())continue;
+    //   if(part->Charge()==0)continue;
+    //   if(TMath::Abs(part->Eta())>0.9)continue;
+    //   list->Add(part);
+    //   iCount++;
+    //   if(part->Pt()>ptmax){ ptmax=part->Pt();
+    //         index=iCount-1;}}}
+      return index;
+}
+
+
+
+Int_t  AliAnalysisTaskJetCore::SelectTrigger(TList *list,Double_t minT,Double_t maxT,Int_t &number){
      Int_t iCount = 0;
      AliAODEvent *aod = 0;
      if(!fESD)aod = fAODIn;
      else aod = fAODOut;   
+     if(!aod)return 0;
      Int_t index=-1;
-     Double_t ptmax=-10;
-    for(int it = 0;it < aod->GetNumberOfTracks();++it){
+     Int_t triggers[100];
+     for(Int_t cr=0;cr<100;cr++){triggers[cr]=-1;}
+     Int_t im=0;
+     for(int it = 0;it < aod->GetNumberOfTracks();++it){
       AliAODTrack *tr = aod->GetTrack(it);
+      Bool_t bGood = false;
+      if(fFilterType == 0)bGood = true;
+      else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
+      else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+       if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
+      if(bGood==false) continue;
+      if (fApplySharedClusterCut) {
+           Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
+           if (frac > 0.4) continue;
+      }
       if(TMath::Abs(tr->Eta())>0.9)continue;
       if(tr->Pt()<0.15)continue;
       list->Add(tr);
       iCount++;
-      if(tr->Pt()>ptmax){ ptmax=tr->Pt();
-      index=iCount-1;}
       
-    }
+      if(tr->Pt()>=minT && tr->Pt()<maxT){
+       triggers[im]=iCount-1;
+        im=im+1;}
+
+     }
+      number=im;
+      Int_t rd=0;
+      if(im==0) rd=0;
+      if(im>0) rd=fRandom->Integer(im);
+      index=triggers[rd];
+
+     
   
    
-  return index;
+  
+      return index;
  
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
  
     AliAODEvent *aod = 0;
@@ -954,7 +1289,7 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
     Int_t index=-1;
     Double_t ptmax=-10;
     Double_t dphi=0;
-    Double_t dif=0;
+    // Double_t dif=0;
     Int_t iCount=0;
     for(int it = 0;it < aod->GetNumberOfTracks();++it){
       AliAODTrack *tr = aod->GetTrack(it);
@@ -963,10 +1298,11 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
       if(tr->Pt()<0.15)continue;
       iCount=iCount+1;
       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
-      if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
+      if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
       index=iCount-1;
-      dif=dphi;  }}
+      //  dif=dphi; 
+      }}
   
       return index;
 
@@ -987,7 +1323,7 @@ Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
      if(!fESD)aod = fAODIn;
      else aod = fAODOut;   
   
-    for(int it = 0;it < aod->GetNumberOfTracks();++it){
+      for(int it = 0;it < aod->GetNumberOfTracks();++it){
       AliAODTrack *tr = aod->GetTrack(it);
       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
       if(TMath::Abs(tr->Eta())>0.9)continue;