]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCore.cxx
select good quality high pt tracks via the filterbit
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
index 394b798faa68700c3ee7df10e70c8b0829ba66e1..09e8cc5f7124aea9554a1243351900152b52fa9f 100644 (file)
@@ -45,6 +45,7 @@
 #include "AliAnalysisHelperJetTasks.h"
 #include "AliInputEventHandler.h"
 #include "AliAODJetEventBackground.h"
+#include "AliAODMCParticle.h"
 #include "AliAnalysisTaskFastEmbedding.h"
 #include "AliAODEvent.h"
 #include "AliAODHandler.h"
@@ -52,9 +53,6 @@
 
 #include "AliAnalysisTaskJetCore.h"
 
-using std::cout;
-using std::endl;
-
 ClassImp(AliAnalysisTaskJetCore)
 
 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
@@ -73,6 +71,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(16),
+fFilterType(2),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -84,9 +84,12 @@ fCheckMethods(0),
 fDoEventMixing(0), 
 fFlagPhiBkg(0),
 fFlagEtaBkg(0),
+fFlagJetHadron(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(100),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -116,7 +119,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -134,10 +136,12 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),
 fh2Ntriggers(0x0),
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
+fh2Ntriggers2(0x0), 
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
 fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
+fh2RPT(0x0),
+fh3spectriggeredC20RP(0x0),
 fh3spectriggeredC20(0x0),
 fh3spectriggeredC3060(0x0)
 
@@ -180,6 +184,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(16),
+fFilterType(2),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -191,9 +197,12 @@ fCheckMethods(0),
 fDoEventMixing(0),
 fFlagPhiBkg(0),
 fFlagEtaBkg(0),
+fFlagJetHadron(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(100),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -223,7 +232,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -241,10 +249,12 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),    
 fh2Ntriggers(0x0),
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
+fh2Ntriggers2(0x0),
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
 fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
+fh2RPT(0x0),
+fh3spectriggeredC20RP(0x0),
 fh3spectriggeredC20(0x0),
 fh3spectriggeredC3060(0x0)
 
@@ -349,7 +359,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);
@@ -359,9 +368,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.);
@@ -383,12 +391,15 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
    
 
     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.);
+    fh2Ntriggers2=new TH2F("# of triggers2","",50,0.,50.,50,0.,50.);
+
+    fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
+    fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
+    fh2RPJets=new TH2F("RPJet","",fNRPBins,0.,fNRPBins,150,0.,150.);
+    fh2RPT=new TH2F("RPTrigger","",fNRPBins,0.,fNRPBins,150,0.,150.); 
+    fh3spectriggeredC20RP = new TH3F("Triggered spectrumC20RP","",3,0.,3.,140,-80.,200.,10,0.,50.);
+    fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,10,0.,50.);
+    fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
 
     
     
@@ -413,7 +424,7 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
       
       fOutputList->Add(fh3JetTrackC3060);
       fOutputList->Add(fh3JetTrackC20);
-      fOutputList->Add(fh3JetTrackC4080); 
+     
             
      
 
@@ -441,10 +452,12 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 
  
        fOutputList->Add(fh2Ntriggers);
-        fOutputList->Add(fh2JetDensity);
-        fOutputList->Add(fh2JetDensityA4);
+        fOutputList->Add(fh2Ntriggers2);
+        fOutputList->Add(fh3JetDensity);
+        fOutputList->Add(fh3JetDensityA4);
         fOutputList->Add(fh2RPJets);
-       fOutputList->Add(fh3spectriggeredC4080);
+        fOutputList->Add(fh2RPT);
+        fOutputList->Add(fh3spectriggeredC20RP); 
        fOutputList->Add(fh3spectriggeredC20); 
        fOutputList->Add(fh3spectriggeredC3060);   
 
@@ -499,9 +512,6 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
     }}
     
 
-
-
-
    // -- event selection --
    fHistEvtSelection->Fill(1); // number of events before event selection
 
@@ -618,10 +628,7 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
      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);
@@ -645,12 +652,15 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    Int_t trigBBTrack=-1;
    Int_t trigInTrack=-1;
    fRPAngle = aod->GetHeader()->GetEventplane();     
-
    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
    if(!partback){  
    PostData(1, fOutputList);
    return;}
+
    fh2Ntriggers->Fill(centValue,partback->Pt());
+   Int_t phiBinT = GetPhiBin(RelativePhi(partback->Phi(),fRPAngle));
+   if(centValue<20.) fh2RPT->Fill(phiBinT,partback->Pt()); 
+
    Double_t accep=2.*TMath::Pi()*1.8;
    Int_t injet4=0;
    Int_t injet=0; 
@@ -660,54 +670,69 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
            phibig  = jetbig->Phi();
            ptbig   = jetbig->Pt();
            if(ptbig==0) continue; 
-           Int_t phiBin = GetPhiBin(phibig-fRPAngle);       
+           Int_t phiBin = GetPhiBin(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(fFlagEtaBkg!=0){
+           if(fFlagEtaBkg==1){
           Double_t etadif= partback->Eta()-etabig;
            if(TMath::Abs(etadif)<=0.5){             
-           if(centValue>40. && centValue<80.) fh3JetTrackC4080->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
+          
            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>40. && centValue<80.) fh3JetTrackC4080->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
            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(fFlagJetHadron==0){
+           if(fFlagPhiBkg!=0) 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(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 
 
-           if(fFlagPhiBkg!=0) 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);
+                
                    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.)  fh3spectriggeredC20RP->Fill(phiBin,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());
 
                    if(ptcorr<=0) continue;
-                       AliAODTrack* leadtrack; 
+
+                  
+                       AliAODTrack* leadtrack=0; 
                        Int_t ippt=0;
-                       Double_t ppt=-10;   
+                       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;
+                      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){                        
@@ -746,19 +771,26 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
                  if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
        
                   
-         if(fDoEventMixing==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();
-         
+          if(centValue<20.) fh2Ntriggers2->Fill(partback->Pt(),part->Pt());
+          if(fDoEventMixing==0 && fFlagOnlyRecoil==0){ 
           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
          
 
@@ -770,8 +802,8 @@ 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
 
 
@@ -803,6 +835,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.
@@ -940,26 +973,60 @@ void AliAnalysisTaskJetCore::Terminate(const Option_t *)
 
 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
 
-     Int_t iCount = 0;
-     AliAODEvent *aod = 0;
+      Int_t iCount = 0;
+      AliAODEvent *aod = 0;
+
      if(!fESD)aod = fAODIn;
      else aod = fAODOut;   
-     Int_t index=-1;
+      Int_t index=-1;
      Double_t ptmax=-10;
-    for(int it = 0;it < aod->GetNumberOfTracks();++it){
+
+
+    
+     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(bGood==false) 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(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;
+       }
+      }
+     }
   
    
-  return index;
+    // 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;
  
 }
 
@@ -980,7 +1047,7 @@ 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;  }}
@@ -1004,7 +1071,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;