]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliAnalysisTaskFragmentationFunction.cxx
Up from Salvatore
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskFragmentationFunction.cxx
index 28540a876a91b5a1de52dddbaece41fa2748b3df..3eb357b17ee79c763f84788a9f99e227e520df29 100644 (file)
@@ -72,6 +72,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
    ,fUseAODInputJets(kTRUE)
    ,fFilterMask(0)
    ,fUsePhysicsSelection(kTRUE)
+   ,fEvtSelectionMask(0)
    ,fEventClass(0)
    ,fMaxVertexZ(10)
    ,fTrackPtCut(0)
@@ -166,7 +167,8 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
    ,fFFXiMax(0)        
    ,fFFNBinsZ(0)       
    ,fFFZMin(0)         
-   ,fFFZMax(0)         
+   ,fFFZMax(0)
+   ,fFFLogZBins(kTRUE)         
    ,fQAJetNBinsPt(0)   
    ,fQAJetPtMin(0)     
    ,fQAJetPtMax(0)     
@@ -376,6 +378,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fUseAODInputJets(kTRUE)
   ,fFilterMask(0)
   ,fUsePhysicsSelection(kTRUE)
+  ,fEvtSelectionMask(0)
   ,fEventClass(0)
   ,fMaxVertexZ(10)
   ,fTrackPtCut(0)
@@ -471,6 +474,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fFFNBinsZ(0)       
   ,fFFZMin(0)         
   ,fFFZMax(0)         
+  ,fFFLogZBins(kTRUE)         
   ,fQAJetNBinsPt(0)   
   ,fQAJetPtMin(0)     
   ,fQAJetPtMax(0)     
@@ -684,6 +688,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fUseAODInputJets(copy.fUseAODInputJets)
   ,fFilterMask(copy.fFilterMask)
   ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
+  ,fEvtSelectionMask(copy.fEvtSelectionMask)
   ,fEventClass(copy.fEventClass)
   ,fMaxVertexZ(copy.fMaxVertexZ)
   ,fTrackPtCut(copy.fTrackPtCut)
@@ -779,6 +784,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fFFNBinsZ(copy.fFFNBinsZ)       
   ,fFFZMin(copy.fFFZMin)         
   ,fFFZMax(copy.fFFZMax)         
+  ,fFFLogZBins(copy.fFFLogZBins)         
   ,fQAJetNBinsPt(copy.fQAJetNBinsPt)   
   ,fQAJetPtMin(copy.fQAJetPtMin)     
   ,fQAJetPtMax(copy.fQAJetPtMax)     
@@ -993,6 +999,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper
     fUseAODInputJets              = o.fUseAODInputJets;
     fFilterMask                   = o.fFilterMask;
     fUsePhysicsSelection          = o.fUsePhysicsSelection;
+    fEvtSelectionMask             = o.fEvtSelectionMask;
     fEventClass                   = o.fEventClass;
     fMaxVertexZ                   = o.fMaxVertexZ;
     fTrackPtCut                   = o.fTrackPtCut;
@@ -1093,6 +1100,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper
     fFFNBinsZ                     = o.fFFNBinsZ;       
     fFFZMin                       = o.fFFZMin;         
     fFFZMax                       = o.fFFZMax;         
+    fFFLogZBins                   = o.fFFLogZBins;         
     fQAJetNBinsPt                 = o.fQAJetNBinsPt;   
     fQAJetPtMin                   = o.fQAJetPtMin;     
     fQAJetPtMax                   = o.fQAJetPtMax;     
@@ -1307,7 +1315,7 @@ AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const
                                                         Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
                                                         Int_t nPt, Float_t ptMin, Float_t ptMax,
                                                         Int_t nXi, Float_t xiMin, Float_t xiMax,
-                                                        Int_t nZ , Float_t zMin , Float_t zMax )
+                                                        Int_t nZ , Float_t zMin , Float_t zMax, Bool_t useLogZBins)
   : TObject()
   ,fNBinsJetPt(nJetPt)
   ,fJetPtMin(jetPtMin)
@@ -1320,7 +1328,8 @@ AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const
   ,fXiMax(xiMax)   
   ,fNBinsZ(nZ)  
   ,fZMin(zMin)    
-  ,fZMax(zMax)    
+  ,fZMax(zMax)
+  ,fLogZBins(useLogZBins)
   ,fh2TrackPt(0)
   ,fh2Xi(0)
   ,fh2Z(0)
@@ -1345,7 +1354,8 @@ AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const
   ,fXiMax(copy.fXiMax)   
   ,fNBinsZ(copy.fNBinsZ)  
   ,fZMin(copy.fZMin)    
-  ,fZMax(copy.fZMax)    
+  ,fZMax(copy.fZMax)
+  ,fLogZBins(copy.fLogZBins)
   ,fh2TrackPt(copy.fh2TrackPt)
   ,fh2Xi(copy.fh2Xi)
   ,fh2Z(copy.fh2Z)
@@ -1374,6 +1384,7 @@ AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& AliAnalysisTaskFragment
     fNBinsZ     = o.fNBinsZ;  
     fZMin       = o.fZMin;    
     fZMax       = o.fZMax;    
+    fLogZBins   = o.fLogZBins;
     fh2TrackPt  = o.fh2TrackPt;
     fh2Xi       = o.fh2Xi;
     fh2Z        = o.fh2Z;
@@ -1403,7 +1414,17 @@ void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::DefineHistos()
   fh1JetPt   = new TH1F(Form("fh1FFJetPt%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
   fh2TrackPt = new TH2F(Form("fh2FFTrackPt%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsPt, fPtMin, fPtMax);
   fh2Xi      = new TH2F(Form("fh2FFXi%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsXi, fXiMin, fXiMax);
-  fh2Z       = new TH2F(Form("fh2FFZ%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, fZMin, fZMax);
+
+  if(!fLogZBins) fh2Z       = new TH2F(Form("fh2FFZ%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, fZMin, fZMax);
+  else{ // logarithmic z binning
+
+    fNBinsZ = fNBinsXi;
+    Double_t binLimsZ[fNBinsXi+1];     
+
+    CalcLogZBins(fNBinsXi,fXiMin,fXiMax,binLimsZ);
+
+    fh2Z    = new TH2F(Form("fh2FFZ%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, binLimsZ);
+  }
 
   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries"); 
   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPt,"jet p_{T} [GeV/c]","p_{T} [GeV/c]","entries");
@@ -1411,6 +1432,32 @@ void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::DefineHistos()
   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z,"jet p_{T} [GeV/c]","z","entries");
 }
 
+//_______________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::CalcLogZBins(const Int_t nBinsXi,const Double_t xiMin,const Double_t xiMax,Double_t* binLims){
+  
+  // calculate logarithmic binning corresponding to equidistant xi bins
+  // expect binLims vector of size nBinsXi+1
+
+  if(nBinsXi == 0){
+    Printf("%s:%d nBinsXi == 0",(char*)__FILE__,__LINE__); 
+    return;
+  }
+
+  Double_t step = (xiMax-xiMin)/nBinsXi;
+  
+  for(Int_t binZ = 0; binZ<nBinsXi; binZ++){
+   
+    Double_t xiUp = xiMax -  binZ*step;
+    Double_t xiLo = xiMax - (binZ+1)*step;
+      
+    Double_t zUp = TMath::Exp(-1*xiLo);
+    Double_t zLo = TMath::Exp(-1*xiUp);
+      
+    if(binZ == 0) binLims[0] = zLo;
+    binLims[binZ+1] = zUp;   
+  }
+}
+
 //_______________________________________________________________________________________________________________
 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::FillFF(Float_t trackPt, Float_t jetPt, Bool_t incrementJetPt, Float_t norm)
 {
@@ -2391,8 +2438,6 @@ Bool_t AliAnalysisTaskFragmentationFunction::Notify()
   return kTRUE;
 }
 
-
-
 //__________________________________________________________________
 void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
 {
@@ -2437,8 +2482,6 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   fJetsEmbedded = new TList();
   fJetsEmbedded->SetOwner(kFALSE);
 
-  // fJetsKine = new TList();
-  // fJetsKine->SetOwner(kTRUE); // delete AOD jets using mom from Kine Tree via TList::Clear()
 
   if(fBckgMode && 
      (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters ||  fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters ||
@@ -2461,7 +2504,8 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
 
   OpenFile(1);
   fCommonHistList = new TList();
-  
+  fCommonHistList->SetOwner(kTRUE);
+
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
   
@@ -2586,30 +2630,31 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   } // end: QA
 
   if(fFFMode){
+
     fFFHistosRecCuts                = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                     fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                     fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+                                                    fFFNBinsZ , fFFZMin , fFFZMax , fFFLogZBins);
     fFFHistosRecLeading        = new AliFragFuncHistos("RecLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                       fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                      fFFNBinsZ , fFFZMin , fFFZMax);
+                                                      fFFNBinsZ , fFFZMin , fFFZMax, fFFLogZBins);
     fFFHistosRecLeadingTrack   = new AliFragFuncHistos("RecLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                       fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                      fFFNBinsZ , fFFZMin , fFFZMax);
+                                                      fFFNBinsZ , fFFZMin , fFFZMax, fFFLogZBins);
     fFFHistosGen            = new AliFragFuncHistos("Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                     fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                     fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                    fFFNBinsZ , fFFZMin , fFFZMax);
+                                                    fFFNBinsZ , fFFZMin , fFFZMax, fFFLogZBins);
     fFFHistosGenLeading        = new AliFragFuncHistos("GenLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                       fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                      fFFNBinsZ , fFFZMin , fFFZMax);
+                                                      fFFNBinsZ , fFFZMin , fFFZMax, fFFLogZBins);
     fFFHistosGenLeadingTrack   = new AliFragFuncHistos("GenLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
                                                       fFFNBinsPt, fFFPtMin, fFFPtMax, 
                                                       fFFNBinsXi, fFFXiMin, fFFXiMax,  
-                                                      fFFNBinsZ , fFFZMin , fFFZMax);
+                                                      fFFNBinsZ , fFFZMin , fFFZMax, fFFLogZBins);
   } // end: FF
 
   if(fIJMode)
@@ -2771,14 +2816,36 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
      
       AliAnalysisTaskFragmentationFunction::SetProperties(fhnResponseJetTrackPt,3,labelsResponseJetTrackPt);
 
-      Int_t    nBinsResponseJetZ[3]     = {fFFNBinsJetPt, fFFNBinsZ,fFFNBinsZ};
-      Double_t binMinResponseJetZ[3]    = {fFFJetPtMin, fFFZMin, fFFZMin};
-      Double_t binMaxResponseJetZ[3]    = {fFFJetPtMax, fFFZMax, fFFZMax};
-      const char* labelsResponseJetZ[3] = { "jet p_{T} [GeV/c]","rec z","gen z"};
+      if(!fFFLogZBins){
 
-      fhnResponseJetZ  = new THnSparseF("fhnResponseJetZ","jet pt:track pt rec:track pt gen",3,
-                                       nBinsResponseJetZ,binMinResponseJetZ,binMaxResponseJetZ);
+       Int_t    nBinsResponseJetZ[3]     = {fFFNBinsJetPt, fFFNBinsZ,fFFNBinsZ};
+       Double_t binMinResponseJetZ[3]    = {fFFJetPtMin, fFFZMin, fFFZMin};
+       Double_t binMaxResponseJetZ[3]    = {fFFJetPtMax, fFFZMax, fFFZMax};
+
+       fhnResponseJetZ  = new THnSparseF("fhnResponseJetZ","jet pt:rec z rec:gen z",3,
+                                         nBinsResponseJetZ,binMinResponseJetZ,binMaxResponseJetZ);
+      }
+      else{
+       
+       Double_t binLims[fFFNBinsXi+1];
+       fFFHistosRecEffGen->CalcLogZBins(fFFNBinsXi,fFFXiMin,fFFXiMax,binLims);
+       
+       Int_t binsZ   = fFFNBinsXi;
+       Double_t zMin = binLims[0];
+       Double_t zMax = binLims[fFFNBinsXi];
+       
+       Int_t    nBinsResponseJetZ[3]     = {fFFNBinsJetPt, binsZ, binsZ};
+       Double_t binMinResponseJetZ[3]    = {fFFJetPtMin,   zMin,  zMin};
+       Double_t binMaxResponseJetZ[3]    = {fFFJetPtMax,   zMax,  zMax};
+       
+       fhnResponseJetZ  = new THnSparseF("fhnResponseJetZ","jet pt:rec z rec:gen z",3,
+                                         nBinsResponseJetZ,binMinResponseJetZ,binMaxResponseJetZ);
+       
+       fhnResponseJetZ->SetBinEdges(1,binLims);
+       fhnResponseJetZ->SetBinEdges(2,binLims);
+      }
       
+      const char* labelsResponseJetZ[3] = { "jet p_{T} [GeV/c]","rec z","gen z"};
       AliAnalysisTaskFragmentationFunction::SetProperties(fhnResponseJetZ,3,labelsResponseJetZ);
       
       Int_t    nBinsResponseJetXi[3]     = {fFFNBinsJetPt, fFFNBinsXi,fFFNBinsXi};
@@ -3183,13 +3250,11 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   fCommonHistList->Add(fh1VertexNContributors);
   fCommonHistList->Add(fh1VertexZ);    
   fCommonHistList->Add(fh1nRecJetsCuts);
-  if(genJets && genTracks){
-    fCommonHistList->Add(fh1Xsec);
-    fCommonHistList->Add(fh1Trials);
-    fCommonHistList->Add(fh1PtHard);
-    fCommonHistList->Add(fh1PtHardTrials);
-    if(genJets) fCommonHistList->Add(fh1nGenJets);
-  }
+  fCommonHistList->Add(fh1Xsec);
+  fCommonHistList->Add(fh1Trials);
+  fCommonHistList->Add(fh1PtHard);
+  fCommonHistList->Add(fh1PtHardTrials);
+  if(genJets) fCommonHistList->Add(fh1nGenJets);
 
   // FF histograms
   if(fFFMode){
@@ -3197,12 +3262,6 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
     fFFHistosRecLeading->AddToOutput(fCommonHistList);
     fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
     if(genJets && genTracks){
-      fCommonHistList->Add(fh1Xsec);
-      fCommonHistList->Add(fh1Trials);
-      fCommonHistList->Add(fh1PtHard);
-      fCommonHistList->Add(fh1PtHardTrials);
-      if(genJets) fCommonHistList->Add(fh1nGenJets);
-
       fFFHistosGen->AddToOutput(fCommonHistList);
       fFFHistosGenLeading->AddToOutput(fCommonHistList);
       fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
@@ -3414,6 +3473,7 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   TH1::AddDirectory(oldStatus);
 
   PostData(1, fCommonHistList);
+
 }
 
 //_______________________________________________
@@ -3433,19 +3493,18 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
        
 
   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
+
   // Trigger selection
   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
-  if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
-    if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
+
+  if(!(inputHandler->IsEventSelected() & fEvtSelectionMask)){
       fh1EvtSelection->Fill(1.);
       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
       PostData(1, fCommonHistList);
       return;
-    }
   }
-
+   
   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
   if(!fESD){
     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
@@ -3566,9 +3625,11 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
 
   //___ get MC information __________________________________________________________________
 
+  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
+
   Double_t ptHard = 0.;
   Double_t nTrials = 1; // trials for MC trigger weight for real data
-  
+
   if(fMCEvent){
     AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
 
@@ -3598,7 +3659,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
        }
       }
       
-      fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
+      //fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
     }
 
   }
@@ -3838,12 +3899,9 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
        }
        
        if(GetFFMinNTracks()>0 && jettracklist->GetSize()<=GetFFMinNTracks()) isBadJet = kTRUE;
-
-       std::cout<<" here, isBadJet "<<isBadJet<<std::endl;
        
        if(isBadJet) continue; 
 
-
        if(ptFractionEmbedded>=fCutFractionPtEmbedded){ // if no embedding: ptFraction = cutFraction = 0
          
          for(Int_t it=0; it<jettracklist->GetSize(); ++it){
@@ -3970,6 +4028,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
        for(Int_t it=0; it<jettracklist->GetSize(); ++it){
          
          AliVParticle*   trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
+         if(!trackVP)continue;
          TLorentzVector* trackV  = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P());
          
          Float_t jetPt   = jet->Pt();
@@ -4952,8 +5011,6 @@ void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist
   sumPt = 0;
   isBadMaxPt = kFALSE;
 
-  std::cout<<" enter GetJetTracksPointing, maxPt "<<maxPt<<std::endl;
-
   Double_t jetMom[3];
   jet->PxPyPz(jetMom);
   TVector3 jet3mom(jetMom);
@@ -4975,13 +5032,9 @@ void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist
       sumPt += track->Pt();
 
       if(maxPt>0 && track->Pt()>maxPt) isBadMaxPt = kTRUE;
-      if(maxPt>0 && track->Pt()>maxPt) std::cout<<" found bad jet, track pt "<<track->Pt()<<" max Pt "<<maxPt<<std::endl;
     }
   }
 
-  std::cout<<" leave GetJetTracksPointing, maxPt "<<maxPt<<" isBadMaxPt "<<isBadMaxPt<<std::endl;
-
-
   outputlist->Sort();
 }
 
@@ -5125,7 +5178,11 @@ void  AliAnalysisTaskFragmentationFunction::FillJetTrackHistosRecGen(TObject* hi
   Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
 
   if(!nTracksJet) return; 
-
+  if(jetPtGen<1e-03){ // check for 0 
+    if(fDebug>0) Printf("%s:%d gen jet pt 0 - return ",(char*)__FILE__,__LINE__);
+    return; 
+  }
+  
   Bool_t incrementJetPtGenFF = kTRUE; // needed in case we fill FFHistos
   Bool_t incrementJetPtRecFF = kTRUE; // needed in case we fill FFHistos