]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskJetCore.cxx
added method to propagate parameters only
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
index 4d39131ef0abdf00df6c6e47dec60f75ef4584af..22f15f249966d311c309e6e5c0c32f13756bf1dd 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"
@@ -66,6 +67,7 @@ fAODExtension(0x0),
 fBackgroundBranch(""),
 fNonStdFile(""),
 fIsPbPb(kTRUE),
+fDebug(0), 
 fOfflineTrgMask(AliVEvent::kAny),
 fMinContribVtx(1),
 fVtxZMin(-10.),
@@ -73,6 +75,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -86,8 +90,11 @@ fFlagPhiBkg(0),
 fFlagEtaBkg(0),
 fFlagJetHadron(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fFlagOnlyHardest(1),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(50),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -117,7 +124,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -135,11 +141,15 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),
 fh2Ntriggers(0x0),
-fh2Ntriggers2(0x0), 
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0), 
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fh3spectriggeredC10(0x0),
 fh3spectriggeredC20(0x0),
 fh3spectriggeredC3060(0x0)
 
@@ -175,6 +185,7 @@ fAODExtension(0x0),
 fBackgroundBranch(""),
 fNonStdFile(""),
 fIsPbPb(kTRUE),
+fDebug(0),
 fOfflineTrgMask(AliVEvent::kAny),
 fMinContribVtx(1),
 fVtxZMin(-10.),
@@ -182,6 +193,8 @@ fVtxZMax(10.),
 fEvtClassMin(0),
 fEvtClassMax(4),
 fFilterMask(0),
+fFilterMaskBestPt(0),
+fFilterType(0),
 fRadioFrac(0.2),
 fMinDist(0.1),
 fCentMin(0.),
@@ -195,8 +208,11 @@ fFlagPhiBkg(0),
 fFlagEtaBkg(0),
 fFlagJetHadron(0),
 fFlagRandom(0),
+fFlagOnlyRecoil(0),
+fFlagOnlyHardest(1),
+fTrackTypeRec(kTrackUndef),
 fRPAngle(0),
-fNRPBins(3),
+fNRPBins(50),
 fJetEtaMin(-.5),
 fJetEtaMax(.5),
 fNevents(0),
@@ -226,7 +242,6 @@ fh2JetCoreMethod1C60(0x0),
 fh2JetCoreMethod2C60(0x0),
 fh3JetTrackC3060(0x0),
 fh3JetTrackC20(0x0),
-fh3JetTrackC4080(0x0), 
 fh2AngStructpt1C10(0x0),
 fh2AngStructpt2C10(0x0),
 fh2AngStructpt3C10(0x0),
@@ -244,11 +259,15 @@ fh2AngStructpt2C60(0x0),
 fh2AngStructpt3C60(0x0),
 fh2AngStructpt4C60(0x0),    
 fh2Ntriggers(0x0),
-fh2Ntriggers2(0x0),
-fh2JetDensity(0x0),
-fh2JetDensityA4(0x0),
-fh2RPJets(0x0),
-fh3spectriggeredC4080(0x0),
+fh2Ntriggers2C10(0x0),
+fh2Ntriggers2C20(0x0),
+fh3JetDensity(0x0),
+fh3JetDensityA4(0x0),
+fh2RPJetsC10(0x0),
+fh2RPJetsC20(0x0),
+fh2RPTC10(0x0),
+fh2RPTC20(0x0), 
+fh3spectriggeredC10(0x0),
 fh3spectriggeredC20(0x0),
 fh3spectriggeredC3060(0x0)
 
@@ -353,7 +372,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);
@@ -363,9 +381,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.);
@@ -387,14 +404,17 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
    
 
     fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
-    fh2Ntriggers2=new TH2F("# of triggers2","",100,0.,4000.,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.);
+    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.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.);
+    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.);  
+    fh3spectriggeredC10 = new TH3F("Triggered spectrumC10","",100,0.,1.,140,-80.,200.,50,0.,50.);
+    fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",100,0.,1.,140,-80.,200.,50,0.,50.);
+    fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",100,0.,1.,140,-80.,200.,10,0.,50.);
 
     
     
@@ -419,7 +439,7 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
       
       fOutputList->Add(fh3JetTrackC3060);
       fOutputList->Add(fh3JetTrackC20);
-      fOutputList->Add(fh3JetTrackC4080); 
+     
             
      
 
@@ -447,13 +467,17 @@ void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 
  
        fOutputList->Add(fh2Ntriggers);
-        fOutputList->Add(fh2Ntriggers2);
-        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);
+        fOutputList->Add(fh3spectriggeredC10); 
+        fOutputList->Add(fh3spectriggeredC20); 
+        fOutputList->Add(fh3spectriggeredC3060);   
 
      
    // =========== Switch on Sumw2 for all histos ===========
@@ -506,16 +530,13 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
     }}
     
 
-
-
-
    // -- event selection --
    fHistEvtSelection->Fill(1); // number of events before event selection
 
    // physics 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);
@@ -625,10 +646,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);
@@ -653,44 +671,49 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
    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());
-   fh2Ntriggers2->Fill(ParticleList.GetEntries(),partback->Pt());
 
+
+   //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
+   //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
+   //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
+
+   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());
    Double_t accep=2.*TMath::Pi()*1.8;
    Int_t injet4=0;
    Int_t injet=0; 
+
    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(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;}
@@ -698,18 +721,19 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 
 
-                   if(centValue<20.) fh2RPJets->Fill(phiBin, ptcorr);
+          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<10.)  fh3spectriggeredC10->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());
 
                    if(ptcorr<=0) continue;
 
-                  
                        AliAODTrack* leadtrack=0; 
                        Int_t ippt=0;
                        Double_t ppt=-10;
@@ -770,14 +794,13 @@ 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) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());                  
+                  if(centValue<20) 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();
@@ -787,8 +810,11 @@ void AliAnalysisTaskJetCore::UserExec(Option_t *)
           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
          
 
@@ -800,11 +826,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){
@@ -833,6 +859,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.
@@ -970,33 +997,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;
+     Double_t ptmax=-10;
 
-     if(!aod)return iCount;
 
-     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);
-      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
+      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(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;
+       }
+      }
+     }
   
    
-  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;
  
 }