Changes to fast embedding and jetresponse: QA plot for delta AOD, Jet eta phi selecti...
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 27 Feb 2011 22:19:21 +0000 (22:19 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 27 Feb 2011 22:19:21 +0000 (22:19 +0000)
JETAN/AliAnalysisTaskFastEmbedding.cxx
JETAN/AliAnalysisTaskFastEmbedding.h
JETAN/AliAnalysisTaskJetCluster.cxx
PWG4/JetTasks/AliAnalysisTaskJetResponse.cxx
PWG4/JetTasks/AliAnalysisTaskJetResponse.h
PWG4/macros/AddTaskFastEmbedding.C
PWG4/macros/AddTaskJetResponse.C
PWG4/macros/AnalysisTrainPWG4Jets.C

index c4ec235..dfce88b 100644 (file)
@@ -37,6 +37,7 @@
 
 #include "AliAnalysisTaskFastEmbedding.h"
 #include "AliAnalysisManager.h"
+#include "AliESDtrack.h"
 #include "AliAODEvent.h"
 #include "AliAODTrack.h"
 #include "AliAODJet.h"
@@ -64,6 +65,10 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
       ,fEvtSelecMode(0)
       ,fEvtSelMinJetPt(-1)
       ,fEvtSelMaxJetPt(-1)
+      ,fEvtSelMinJetEta(-999.)
+      ,fEvtSelMaxJetEta( 999.)
+      ,fEvtSelMinJetPhi(0.)
+      ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
       ,fToyMinNbOfTracks(1)
       ,fToyMaxNbOfTracks(1)
       ,fToyMinTrackPt(50.)
@@ -78,9 +83,13 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
       ,fh1TrackPt(0)
       ,fh2TrackEtaPhi(0)
       ,fh1TrackN(0)
+      ,fh1JetPt(0)
+      ,fh2JetEtaPhi(0)
+      ,fh1JetN(0)
       ,fh1MCTrackPt(0)
       ,fh2MCTrackEtaPhi(0)
       ,fh1MCTrackN(0)
+      ,fh1AODfile(0)
 
 
 {
@@ -107,6 +116,10 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
       ,fEvtSelecMode(0)
       ,fEvtSelMinJetPt(-1)
       ,fEvtSelMaxJetPt(-1)
+      ,fEvtSelMinJetEta(-999.)
+      ,fEvtSelMaxJetEta( 999.)
+      ,fEvtSelMinJetPhi(0.)
+      ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
       ,fToyMinNbOfTracks(1)
       ,fToyMaxNbOfTracks(1)
       ,fToyMinTrackPt(50.)
@@ -121,9 +134,13 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
       ,fh1TrackPt(0)
       ,fh2TrackEtaPhi(0)
       ,fh1TrackN(0)
+      ,fh1JetPt(0)
+      ,fh2JetEtaPhi(0)
+      ,fh1JetN(0)
       ,fh1MCTrackPt(0)
       ,fh2MCTrackEtaPhi(0)
       ,fh1MCTrackN(0)
+      ,fh1AODfile(0)
 {
     // constructor
     DefineOutput(1, TList::Class());
@@ -148,6 +165,10 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTask
       ,fEvtSelecMode(copy.fEvtSelecMode)
       ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
       ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
+      ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
+      ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
+      ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
+      ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
       ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
       ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
       ,fToyMinTrackPt(copy.fToyMinTrackPt)
@@ -162,9 +183,13 @@ AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTask
       ,fh1TrackPt(copy.fh1TrackPt)
       ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
       ,fh1TrackN(copy.fh1TrackN)
+      ,fh1JetPt(copy.fh1JetPt)
+      ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
+      ,fh1JetN(copy.fh1JetN)
       ,fh1MCTrackPt(copy.fh1MCTrackPt)
       ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
       ,fh1MCTrackN(copy.fh1MCTrackN)
+      ,fh1AODfile(copy.fh1AODfile)
 {
     // copy constructor
 }
@@ -192,6 +217,10 @@ AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliA
         fEvtSelecMode      = o.fEvtSelecMode;
         fEvtSelMinJetPt    = o.fEvtSelMinJetPt;
         fEvtSelMaxJetPt    = o.fEvtSelMaxJetPt;
+        fEvtSelMinJetEta   = o.fEvtSelMinJetEta;
+        fEvtSelMaxJetEta   = o.fEvtSelMaxJetEta;
+        fEvtSelMinJetPhi   = o.fEvtSelMinJetPhi;
+        fEvtSelMaxJetPhi   = o.fEvtSelMaxJetPhi;
         fToyMinNbOfTracks  = o.fToyMinNbOfTracks;
         fToyMaxNbOfTracks  = o.fToyMaxNbOfTracks;
         fToyMinTrackPt     = o.fToyMinTrackPt;
@@ -206,9 +235,13 @@ AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliA
         fh1TrackPt         = o.fh1TrackPt;
         fh2TrackEtaPhi     = o.fh2TrackEtaPhi;
         fh1TrackN          = o.fh1TrackN;
+       fh1JetPt           = o.fh1JetPt;
+        fh2JetEtaPhi       = o.fh2JetEtaPhi;
+        fh1JetN            = o.fh1JetN;
         fh1MCTrackPt       = o.fh1MCTrackPt;
         fh2MCTrackEtaPhi   = o.fh2MCTrackEtaPhi;
         fh1MCTrackN        = o.fh1MCTrackN;
+        fh1AODfile         = o.fh1AODfile;
     }
 
     return *this;
@@ -229,28 +262,6 @@ void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
     // create output objects
     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
 
-    rndm = new TRandom3();
-    Int_t id = GetJobID();
-    if(id>-1) rndm->SetSeed(id);
-    else      rndm->SetSeed();   // a TTUID is generated and used for seed
-    AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
-
-
-    // embed mode with AOD
-    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-
-       // open input AOD
-       if(fAODPathArray){
-           Int_t rc = SelectAODfile();
-           if(rc<0) return;
-       }
-
-       Int_t rc = OpenAODfile();
-       if(rc<0) return;
-
-    } //end: embed mode with AOD
-
-
     // connect output aod
     // create a new branch for extra tracks
     fAODout = AODEvent();
@@ -278,34 +289,79 @@ void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
     //qa histograms
 
     OpenFile(1);
-    fHistList = new TList();
+    if(!fHistList) fHistList = new TList();
+    fHistList->SetOwner(kTRUE);
 
     Bool_t oldStatus = TH1::AddDirectoryStatus();
     TH1::AddDirectory(kFALSE);
 
-    fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 300, 0., 300.);
+    fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
     fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
-    fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 3000.);
+    fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
 
     fHistList->Add(fh1TrackPt);
     fHistList->Add(fh2TrackEtaPhi);
     fHistList->Add(fh1TrackN);
+       
+       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+       
+           fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 250, 0., 250.);
+           fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
+           fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
+       
+           fHistList->Add(fh1JetPt);
+        fHistList->Add(fh2JetEtaPhi);
+        fHistList->Add(fh1JetN);
+    }
 
  
     if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){ 
 
-       fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 300, 0., 300.);
+       fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
        fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
-       fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 3000.);
-
+       fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
+          
        fHistList->Add(fh1MCTrackPt);
        fHistList->Add(fh2MCTrackEtaPhi);
        fHistList->Add(fh1MCTrackN);
-
+         
+    }
+       
+    fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 100, 0, 100);
+    fHistList->Add(fh1AODfile);
+
+    // =========== Switch on Sumw2 for all histos ===========
+    for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
+        TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
+        if (h1){
+          h1->Sumw2();
+          continue;
+        }
     }
 
     TH1::AddDirectory(oldStatus);
 
+
+    // set seed
+    rndm = new TRandom3();
+    Int_t id = GetJobID();
+    if(id>-1) rndm->SetSeed(id);
+    else      rndm->SetSeed();   // a TTUID is generated and used for seed
+    AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
+
+    // embed mode with AOD
+    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+
+       // open input AOD
+       Int_t rc = OpenAODfile();
+       if(rc<0) return;
+       fh1AODfile->Fill(rc);
+
+    } //end: embed mode with AOD
+
+   
+    
+   PostData(1, fHistList);
 }
 
 
@@ -325,6 +381,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
 
     if(!fAODout){
        AliError("Need output AOD, but is not connected."); 
+       PostData(1, fHistList);
        return;
     }
 
@@ -332,6 +389,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
     TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
     if(!tracks){
         AliError("Extra track branch not found in output.");
+        PostData(1, fHistList);
         return;
     }
     tracks->Delete();
@@ -344,6 +402,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
     if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
        if(!fAODevent){
           AliError("Need input AOD, but is not connected."); 
+          PostData(1, fHistList);
           return;
        }
 
@@ -353,6 +412,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
        else                    aodJets = fAODevent->GetJets();
        if(!aodJets){
           AliError("Could not find jets in AOD. Check jet branch when indicated.");
+          PostData(1, fHistList);
           return;
        }
         
@@ -368,17 +428,20 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
               } 
               else {
                  AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
-                 Int_t rc = SelectAODfile();
-                 if(rc<0) return;
 
-                 rc = OpenAODfile();
-                 if(rc<0) return;
+                 Int_t rc = OpenAODfile();
+                 if(rc<0) {
+                    PostData(1, fHistList);
+                    return;
+                 }
+                 fh1AODfile->Fill(rc);
 
                 // new file => we must use the new jet array
                 if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
                 else                    aodJets = fAODevent->GetJets();
                 if(!aodJets){
                   AliError("Could not find jets in AOD. Check jet branch when indicated.");
+                   PostData(1, fHistList);
                   return;
                 }
               }
@@ -394,7 +457,9 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
                  if(!jet) continue;
 
                  if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
-                   && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)){
+                   && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+                    && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+                    && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
                     useEntry = kTRUE;
                     break;
                  } 
@@ -424,18 +489,65 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
        if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
           // loop over input tracks
           // add to output aod
-          Int_t nTracks = fAODevent->GetNumberOfTracks();
+          Int_t nTracks = 0;
+          Int_t nJets = aodJets->GetEntries();
+          Int_t nSelectedJets = 0;
+       AliAODJet *leadJet = 0x0; // used for jet tracks only
+           
+           if(fEmbedMode==kAODFull){
+                          nTracks = fAODevent->GetNumberOfTracks();
+                                  
+                              for(Int_t ij=0; ij<nJets; ++ij){
+                       AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+                       if(!jet) continue;
+                       if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+                          && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+                          && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+                          && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
+                                                 
+                                                     fh1JetPt->Fill(jet->Pt());
+                                                         fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
+                                                         nSelectedJets++;
+                                                         
+                                          }
+                   }                              
+                                  fh1JetN->Fill(nSelectedJets);
+                  }
+
+           if(fEmbedMode==kAODJetTracks){
+              // look for leading jet within selection
+              // get reference tracks for this jet
+              Float_t leadJetPt = 0.;
+              for(Int_t ij=0; ij<nJets; ++ij){
+                  AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+                  if(!jet) continue;
+                  if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+                     && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+                     && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+                     && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
+                     if(jet->Pt()>leadJetPt){
+                         leadJetPt = jet->Pt();
+                         leadJet = jet;
+                     } 
+                  }
+               }
+               if(leadJet){
+                          nTracks = leadJet->GetRefTracks()->GetEntriesFast();
+                                  fh1JetPt->Fill(leadJet->Pt());
+                   fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
+                   fh1JetN->Fill(1);                              
+                          }
+           }
            fh1TrackN->Fill((Float_t)nTracks);
 
           for(Int_t it=0; it<nTracks; ++it){
-              AliAODTrack *tmpTr = fAODevent->GetTrack(it);
-              // ?? test filter bit, or something else??
+              AliAODTrack *tmpTr = 0x0;
+               if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
+               if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
+               if(!tmpTr) continue; 
 
-               if(fEmbedMode==kAODJetTracks){
-                  // jet track? else continue
-                  // not implemented yet
-                  continue;
-               } 
+              tmpTr->SetStatus(AliESDtrack::kEmbedded);
 
               new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
               dummy = (*tracks)[nAODtracks-1];
@@ -482,6 +594,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
                AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
                if(!jet) continue;
               AliAODTrack *tmpTr = (AliAODTrack*)jet;
+              tmpTr->SetFlags(AliESDtrack::kEmbedded);
 
               new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
               dummy = (*tracks)[nAODtracks-1]; 
@@ -521,7 +634,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
           Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
           phi = TVector2::Phi_0_2pi(phi);
 
-          Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
+          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
 
           Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
           Float_t mom[3];
@@ -551,7 +664,7 @@ void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
                                              fToyFilterMap,  // select info
                                              -999.    // chi2 per NDF
                                            );
-          tmpTr->SetFlags(1<<27);
+          tmpTr->SetFlags(AliESDtrack::kEmbedded);
 
            new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
            dummy = (*tracks)[nAODtracks-1];
@@ -614,7 +727,18 @@ Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
      return n;
 }
 
-Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(){
+//__________________________________________________________________________
+
+Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
+
+    if(trial>3){
+        AliError("Could not open AOD files, give up ...");
+        return -1;
+    }
+       
+       Int_t rc = 0;
+       if(fAODPathArray) rc = SelectAODfile();
+       if(rc<0) return -1;
 
     TDirectory *owd = gDirectory;
     if (fAODfile)
@@ -622,8 +746,16 @@ Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(){
     fAODfile = TFile::Open(fAODPath.Data());
     owd->cd();
     if(!fAODfile){
-       AliError("Could not open AOD file.");
-       return -1;
+       
+       rc = -1;
+       if(fAODPathArray){
+           AliError(Form("Could not open AOD file %s, try another from the list ...", fAODPath.Data()));
+           rc = OpenAODfile(trial+1);
+       } else {
+              AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
+          }
+        
+       return rc;
     }
 
     fAODtree = (TTree*)fAODfile->Get("aodTree");
@@ -636,6 +768,6 @@ Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(){
     delete fAODevent;
     fAODevent = new AliAODEvent();
     fAODevent->ReadFromTree(fAODtree);
-    return 1;
+    return rc;  // file position in AOD path array, if array available
 }
 
index f3288bf..db40f3b 100644 (file)
@@ -46,6 +46,8 @@ class AliAnalysisTaskFastEmbedding : public AliAnalysisTaskSE {
        Int_t GetEvtSelecMode() {return fEvtSelecMode;}
 
        void SetEvtSelJetPtRange(Float_t minPt, Float_t maxPt) {fEvtSelMinJetPt = minPt; fEvtSelMaxJetPt = maxPt;}
+       void SetEvtSelJetEtaRange(Float_t minEta, Float_t maxEta) {fEvtSelMinJetEta = minEta; fEvtSelMaxJetEta = maxEta;}
+       void SetEvtSelJetPhiRange(Float_t minPhi, Float_t maxPhi) {fEvtSelMinJetPhi = minPhi; fEvtSelMaxJetPhi = maxPhi;}
        
         void SetToyNumberOfTrackRange(Int_t minN = 1, Int_t maxN = 1){ fToyMinNbOfTracks = minN, fToyMaxNbOfTracks = maxN; }
        void SetToyTrackRanges(Double_t minPt = 50., Double_t maxPt = 50., Double_t ptDistr=0,
@@ -87,7 +89,10 @@ class AliAnalysisTaskFastEmbedding : public AliAnalysisTaskSE {
        // event selection from AOD
        Float_t fEvtSelMinJetPt;       // minimum pt of the leading jet
        Float_t fEvtSelMaxJetPt;       // maximum pt of the leading jet
-       // ... todo: eta, phi, ...
+        Float_t fEvtSelMinJetEta;      // minimum eta of the leading jet
+        Float_t fEvtSelMaxJetEta;      // maximum eta of the leading jet
+        Float_t fEvtSelMinJetPhi;      // minimum phi of the leading jet
+        Float_t fEvtSelMaxJetPhi;      // maximum phi of the leading jet
         
          
         // settings for toy "track generation"
@@ -108,16 +113,21 @@ class AliAnalysisTaskFastEmbedding : public AliAnalysisTaskSE {
         TH1F  *fh1TrackPt;         //! track pt
         TH2F  *fh2TrackEtaPhi;     //! track eta-phi
         TH1F  *fh1TrackN;          //! nb. of tracks
+        TH1F  *fh1JetPt;           //! jet pt
+        TH2F  *fh2JetEtaPhi;       //! jet eta-phi
+        TH1F  *fh1JetN;            //! nb. of jets
         TH1F  *fh1MCTrackPt;       //! MC track pt
         TH2F  *fh2MCTrackEtaPhi;   //! MC track eta-phi
         TH1F  *fh1MCTrackN;        //! nb. of MC tracks
+        TH1I  *fh1AODfile;         //! used AOD files from AODPathArray
+               
 
        Int_t GetJobID();    // get job id (sub-job id on the GRID)
         Int_t SelectAODfile();
-        Int_t OpenAODfile();
+        Int_t OpenAODfile(Int_t trial = 0);
 
 
-       ClassDef(AliAnalysisTaskFastEmbedding, 3);
+       ClassDef(AliAnalysisTaskFastEmbedding, 4);
 };
 
 #endif
index 7cf425e..8b35028 100644 (file)
@@ -967,7 +967,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
    // Add the random cones...
    if(fNRandomCones>0&&fTCARandomConesOut){       
      // create a random jet within the acceptance
-     Double_t etaMax = 0.8 - fRparam;
+     Double_t etaMax = fTrackEtaWindow - fRparam;
      Int_t nCone = 0;
      Int_t nConeRan = 0;
      Double_t pTC = 1; // small number
index 909a4ac..a2d5495 100644 (file)
@@ -3,6 +3,7 @@
 #include "TMath.h"
 #include "TH1F.h"
 #include "TH2F.h"
+#include "TH3F.h"
 #include "TCanvas.h"
 
 #include "AliLog.h"
@@ -32,29 +33,38 @@ AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
   fMinContribVtx(1),
   fVtxZMin(-8.),
   fVtxZMax(8.),
-  fEvtClassMin(0),
-  fEvtClassMax(10),
+  fEvtClassMin(1),
+  fEvtClassMax(4),
   fCentMin(0.),
   fCentMax(100.),
   fJetEtaMin(-.5),
   fJetEtaMax(.5),
-  fJetDeltaEta(.4),
-  fJetDeltaPhi(.4),
+  fJetPtMin(20.),
+  fJetPtFractionMin(0.5),
+  fNMatchJets(4),
+  //fJetDeltaEta(.4),
+  //fJetDeltaPhi(.4),
   fkNbranches(2),
-  fkEvtClasses(10),
+  fkEvtClasses(5),
   fOutputList(0x0),
   fHistEvtSelection(0x0),
-  fHistPtLeadingJet(new TH1F*[fkNbranches]),
-  fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
-  fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
-  fHistDeltaEtaDeltaPhiLeadingJet(0x0),
-  fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
-  fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
+  fHistEvtClass(0x0),
+  fHistCentrality(0x0),
+  fHistPtJet(new TH1F*[fkNbranches]),
+  fHistEtaPhiJet(new TH2F*[fkNbranches]),
+  fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
+  fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
+  fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
+  fHistPtFraction(new TH2F*[fkEvtClasses]),
   fHistPtPtExtra(0x0),
   fHistPtResponse(new TH2F*[fkEvtClasses]),
   fHistPtSmearing(new TH2F*[fkEvtClasses]),
-  fHistdR(new TH2F*[fkEvtClasses]),
-  fHistArea(new TH2F*[fkEvtClasses])
+  fHistDeltaR(new TH2F*[fkEvtClasses]),
+  fHistArea(new TH2F*[fkEvtClasses]),
+  fHistDeltaArea(new TH2F*[fkEvtClasses])
 {
   // default Constructor
 
@@ -73,29 +83,38 @@ AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
   fMinContribVtx(1),
   fVtxZMin(-8.),
   fVtxZMax(8.),
-  fEvtClassMin(0),
-  fEvtClassMax(10),
+  fEvtClassMin(1),
+  fEvtClassMax(4),
   fCentMin(0.),
   fCentMax(100.),
   fJetEtaMin(-.5),
   fJetEtaMax(.5),
-  fJetDeltaEta(.4),
-  fJetDeltaPhi(.4),
+  fJetPtMin(20.),
+  fJetPtFractionMin(0.5),
+  fNMatchJets(4),
+  //fJetDeltaEta(.4),
+  //fJetDeltaPhi(.4),
   fkNbranches(2),
-  fkEvtClasses(10),
+  fkEvtClasses(5),
   fOutputList(0x0),
   fHistEvtSelection(0x0),
-  fHistPtLeadingJet(new TH1F*[fkNbranches]),
-  fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
-  fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
-  fHistDeltaEtaDeltaPhiLeadingJet(0x0),
-  fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
-  fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
+  fHistEvtClass(0x0),
+  fHistCentrality(0x0),
+  fHistPtJet(new TH1F*[fkNbranches]),
+  fHistEtaPhiJet(new TH2F*[fkNbranches]),
+  fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
+  fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
+  fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
+  fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
+  fHistPtFraction(new TH2F*[fkEvtClasses]),
   fHistPtPtExtra(0x0),
   fHistPtResponse(new TH2F*[fkEvtClasses]),
   fHistPtSmearing(new TH2F*[fkEvtClasses]),
-  fHistdR(new TH2F*[fkEvtClasses]),
-  fHistArea(new TH2F*[fkEvtClasses])
+  fHistDeltaR(new TH2F*[fkEvtClasses]),
+  fHistArea(new TH2F*[fkEvtClasses]),
+  fHistDeltaArea(new TH2F*[fkEvtClasses])
 {
   // Constructor
 
@@ -134,69 +153,111 @@ void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
 {
   // Create histograms
   // Called once
+  OpenFile(1);
+  if(!fOutputList) fOutputList = new TList;
+  fOutputList->SetOwner(kTRUE);
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
 
-  fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
+  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
+  
+  fHistEvtClass   = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
+  fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
 
   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
-    fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
-                                         Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
-                                         300, 0., 300.);
-    fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
-
-    fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
-                                             Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
-                                             100, -2., 2., 100, 0., 2*TMath::Pi());
+    fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
+                                         Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
+                                         250, 0., 250.);
+    fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
+
+    fHistEtaPhiJet[iBranch]    = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
+                                                        Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
+                                                        48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
+       fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
+                                                        Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
+                                                        48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
 
   }
-  fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
-                                            "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
-                                            100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());
+  
 
   fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
-                           300, 0., 300., 300, 0., 300.);
+                           250, 0., 250., 250, 0., 250.);
 
   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
+    fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
+                                                                     "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
+                                                                     101, -1.01, 1.01, 101, -1.01, 1.01);
+    fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
+                                                                     "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
+                                                                     101, -1.01, 1.01, 101, -1.01, 1.01);                                                                                                
+    fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
+                                                                                "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
+                                                                                100, -2., 2., 100, TMath::Pi(), TMath::Pi());
     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
-                                         300, 0., 300., 300, 0., 300.);
-    fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
-                                         200, -50., 150., 300, 0., 300.);
-    fHistdR[iEvtClass] = new TH2F(Form("hist_dR%i",iEvtClass), "", 200, -50.,150., 200, -1.,1.);
-    fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "", 200, -50.,150., 100, 0.,1.);
-
-    fHistEtaPhiLeadingJetCut[iEvtClass] = new TH2F(Form("fHistEtaPhiLeadingJetCut%i", iEvtClass),
-                                             Form("#eta - #phi of leading jet from event class %i;#eta;#phi", iEvtClass),
-                                             100, -2., 2., 100, 0., 2*TMath::Pi());
-    fHistDeltaEtaEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaLeadingJet%i", iEvtClass),
-                                                Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#eta", iEvtClass),
-                                                100, -1., 1., 100, -.5, .5);
-    fHistDeltaPtEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaLeadingJet%i", iEvtClass),
-                                                Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#p_{T}", iEvtClass),
-                                                100, -1., 1., 200, -50., 150);
+                                                             250, 0., 250., 250, 0., 250.);
+    fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
+                                                             200, -50., 150., 250, 0., 250.);
+    fHistDeltaR[iEvtClass]     = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
+    fHistArea[iEvtClass]       = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
+    fHistDeltaArea[iEvtClass]  = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3); 
+       
+    fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
+                                                Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
+                                                60, -.6, .6, 100, -.5, .5);
+    fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
+                                                Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
+                                                60, -.6, .6, 200, -50., 150);
+    fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
+                                             Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
+                                                                                 100, 0., 1., 250, 0., 250.);
   }
 
-  OpenFile(1);
-  fOutputList = new TList;
+  
   fOutputList->Add(fHistEvtSelection);
+  fOutputList->Add(fHistEvtClass);
+  fOutputList->Add(fHistCentrality);
+
   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
-    fOutputList->Add(fHistPtLeadingJet[iBranch]);
-    fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
+    fOutputList->Add(fHistPtJet[iBranch]);
+    fOutputList->Add(fHistEtaPhiJet[iBranch]);
+       fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
   }
-  fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
+  
   fOutputList->Add(fHistPtPtExtra);
+  
   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
+    fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
+    fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
+    fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
+       fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
+    fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
+       fOutputList->Add(fHistPtFraction[iEvtClass]);
+       
     fOutputList->Add(fHistPtResponse[iEvtClass]);
     fOutputList->Add(fHistPtSmearing[iEvtClass]);
-    fOutputList->Add(fHistdR[iEvtClass]);
+    fOutputList->Add(fHistDeltaR[iEvtClass]);
     fOutputList->Add(fHistArea[iEvtClass]);
-    fOutputList->Add(fHistEtaPhiLeadingJetCut[iEvtClass]);
-    fOutputList->Add(fHistDeltaEtaEtaLeadingJet[iEvtClass]);
-    fOutputList->Add(fHistDeltaPtEtaLeadingJet[iEvtClass]);
+    fOutputList->Add(fHistDeltaArea[iEvtClass]);
   }
+
+  // =========== Switch on Sumw2 for all histos ===========
+  for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
+    TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
+    if (h1){
+      h1->Sumw2();
+      continue;
+    }
+  }
+  TH1::AddDirectory(oldStatus);
+
+  PostData(1, fOutputList);
 }
 
 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
@@ -219,13 +280,14 @@ void AliAnalysisTaskJetResponse::UserExec(Option_t *)
     return;
   }
 
+  // -- event selection --
   fHistEvtSelection->Fill(1); // number of events before event selection
 
   // physics selection
   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
-    Printf(" Trigger Selection: event REJECTED ... ");
+    if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
     fHistEvtSelection->Fill(2);
     PostData(1, fOutputList);
     return;
@@ -237,7 +299,7 @@ void AliAnalysisTaskJetResponse::UserExec(Option_t *)
   if ((nTracksPrim < fMinContribVtx) ||
       (primVtx->GetZ() < fVtxZMin) ||
       (primVtx->GetZ() > fVtxZMax) ){
-    Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+    if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
     fHistEvtSelection->Fill(3);
     PostData(1, fOutputList);
     return;
@@ -245,7 +307,7 @@ void AliAnalysisTaskJetResponse::UserExec(Option_t *)
 
   // event class selection (from jet helper task)
   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
-  Printf("Event class %d", eventClass);
+  if(fDebug) Printf("Event class %d", eventClass);
   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
     fHistEvtSelection->Fill(4);
     PostData(1, fOutputList);
@@ -255,84 +317,131 @@ void AliAnalysisTaskJetResponse::UserExec(Option_t *)
   // centrality selection
   AliCentrality *cent = fESD->GetCentrality();
   Float_t centValue = cent->GetCentralityPercentile("TRK");
-  printf("centrality: %f\n", centValue);
+  if(fDebug) printf("centrality: %f\n", centValue);
   if (centValue < fCentMin || centValue > fCentMax){
-    fHistEvtSelection->Fill(5);
+    fHistEvtSelection->Fill(4);
     PostData(1, fOutputList);
     return;
   }
 
   fHistEvtSelection->Fill(0); // accepted events
+  fHistEvtClass->Fill(eventClass);
+  fHistCentrality->Fill(centValue);
+  // -- end event selection --
 
-
+  // fetch jets
   TClonesArray *aodJets[2];
   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
-  AliAODJet *leadingJet[2] = { 0x0, 0x0 };
 
   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
     fListJets[iJetType]->Clear();
-    if (!aodJets[iJetType])
-      continue;
+    if (!aodJets[iJetType]) continue;
+       
     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
-      if (jet)
-       fListJets[iJetType]->Add(jet);
+      if (jet) fListJets[iJetType]->Add(jet);
     }
     fListJets[iJetType]->Sort();
-    leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
-
-    if (leadingJet[iJetType])
-      fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
   }
-
-  // check if leading jets are close in eta-phi and compare their Pt
-  if (leadingJet[0] && leadingJet[1]) {
-    Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
-    Float_t deltaPhi = TVector2::Phi_mpi_pi(leadingJet[0]->Phi() - leadingJet[1]->Phi());
-
-    if (eventClass > -1 && eventClass < fkEvtClasses){
-
-       if(leadingJet[1]->Eta()<fJetEtaMax && leadingJet[1]->Eta()>fJetEtaMin){
-          fHistEtaPhiLeadingJetCut[eventClass]->Fill(leadingJet[1]->Eta(), leadingJet[1]->Phi());
-       }
-    }
-
-    // leading jets in "jet" acceptance
-    if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
-       leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
-       Printf("Jet not in eta acceptance.");
-    }
-    else{
-       // check association of jets
-       fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
-
-       if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
-          printf("leading jets two far apart\n");
-       else {
-          Float_t pt0 = leadingJet[0]->Pt();
-          Float_t pt1 = leadingJet[1]->Pt();
-          Float_t dPt = pt1-pt0;
-          Float_t jetarea = leadingJet[1]->EffectiveAreaCharged();
-
-          fHistPtLeadingJet[0]->Fill(pt0);
-          fHistPtLeadingJet[1]->Fill(pt1);
-
-          fHistPtPtExtra->Fill(pt0, pt1);
-
-          if (eventClass > -1 && eventClass < fkEvtClasses){
-             fHistPtResponse[eventClass]->Fill(pt1, pt0);
-              fHistPtSmearing[eventClass]->Fill(dPt, pt0);
-
-              Float_t dR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
-              fHistdR[eventClass]->Fill(dPt, dR);
-              fHistArea[eventClass]->Fill(dPt, jetarea);
-
-              fHistDeltaEtaEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), deltaEta);
-              fHistDeltaPtEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), dPt);
-          }
-       }
-    }
+  
+  // jet matching 
+  static TArrayI aMatchIndex(fListJets[0]->GetEntries());
+  static TArrayF aPtFraction(fListJets[0]->GetEntries());
+  if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
+  if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
+  
+  // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
+  AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
+                                            fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
+                                            aMatchIndex, aPtFraction, fDebug);
+                                                                                       
+  // loop over matched jets
+  Int_t ir = -1; // index of matched reconstruced jet
+  Float_t fraction = 0.;
+  AliAODJet *jet[2]  = { 0x0, 0x0 };
+  Float_t jetEta[2]  = { 0., 0. };
+  Float_t jetPhi[2]  = { 0., 0. };
+  Float_t jetPt[2]   = { 0., 0. };
+  Float_t jetArea[2] = { 0., 0. };
+   
+  for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
+     ir = aMatchIndex[ig];
+        if(ir<0) continue;
+        fraction = aPtFraction[ig];
+        
+        // fetch jets
+        jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
+        jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
+        if(!jet[0] || !jet[1]) continue;
+        
+        for(Int_t i=0; i<fkNbranches; ++i){
+           jetEta[i]  = jet[i]->Eta();
+               jetPhi[i]  = jet[i]->Phi();
+               jetPt[i]   = jet[i]->Pt();
+               jetArea[i] = jet[i]->EffectiveAreaCharged();
+        }
+        if(eventClass > -1 && eventClass < fkEvtClasses){
+        fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
+        }
+        
+        if(fraction<fJetPtFractionMin) continue;
+        
+        // calculate parameters of associated jets
+        Float_t deltaPt    = jetPt[1]-jetPt[0];
+        Float_t deltaEta   = jetEta[1]-jetEta[0];
+        Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
+        Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
+        Float_t deltaArea  = jetArea[1]-jetArea[0];
+        
+        // fill histograms before acceptance cut
+        for(Int_t i=0; i<fkNbranches; ++i){
+            fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
+         }
+        if(eventClass > -1 && eventClass < fkEvtClasses){
+             fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
+        }
+        
+        // jet acceptance + minimum pT check
+        if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
+           jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
+               
+               if(fDebug){
+               Printf("Jet not in eta acceptance.");
+                       Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
+                       Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
+               }
+               continue;
+     }
+        if(jetPt[1] < fJetPtMin){
+           if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
+           continue;
+        }
+        
+        
+        
+        // fill histograms
+        for(Int_t i=0; i<fkNbranches; ++i){
+            fHistPtJet[i]        -> Fill(jetPt[i]);
+             fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
+        }
+        
+        fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
+        
+        if(eventClass > -1 && eventClass < fkEvtClasses){
+                fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
+                
+                fHistPtResponse[eventClass]             -> Fill(jetPt[1], jetPt[0]);
+                fHistPtSmearing[eventClass]             -> Fill(deltaPt,  jetPt[0]);
+                
+                fHistDeltaPtEtaJet[eventClass]          -> Fill(jetEta[0], deltaPt);
+                 fHistDeltaEtaEtaJet[eventClass]         -> Fill(jetEta[0], deltaEta);
+                
+                fHistDeltaR[eventClass]                 -> Fill(deltaPt, deltaR);
+                fHistArea[eventClass]                   -> Fill(deltaPt, jetArea[1]);
+                 fHistDeltaArea[eventClass]              -> Fill(deltaPt, deltaArea);
+                
+        }
   }
 
   PostData(1, fOutputList);
@@ -346,3 +455,4 @@ void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
   if (!GetOutputData(1))
     return;
 }
+
index 759eaf8..9b6d539 100644 (file)
@@ -3,6 +3,7 @@
 
 class TH1F;
 class TH2F;
+class TH3F;
 class AliESDEvent;
 class AliAODEvent;
 
@@ -23,6 +24,7 @@ class AliAnalysisTaskJetResponse : public AliAnalysisTaskSE {
 
   virtual AliVEvent::EOfflineTriggerTypes GetOfflineTrgMask() const { return fOfflineTrgMask; }
   virtual void     GetBranchNames(TString &branch1, TString &branch2) const { branch1 = fJetBranchName[0]; branch2 = fJetBranchName[1]; }
+  virtual Int_t    GetMinContribVtx() const { return fMinContribVtx; };
   virtual Float_t  GetVtxZMin() const { return fVtxZMin; }
   virtual Float_t  GetVtxZMax() const { return fVtxZMax; }
   virtual Int_t    GetEvtClassMin() const { return fEvtClassMin; }
@@ -31,11 +33,15 @@ class AliAnalysisTaskJetResponse : public AliAnalysisTaskSE {
   virtual Float_t  GetCentMax() const { return fCentMax; }
   virtual Float_t  GetJetEtaMin() const { return fJetEtaMin; }
   virtual Float_t  GetJetEtaMax() const { return fJetEtaMax; }
-  virtual Float_t  GetJetDeltaEta() const { return fJetDeltaEta; }
-  virtual Float_t  GetJetDeltaPhi() const { return fJetDeltaPhi; }
+  virtual Float_t  GetJetPtMin() const { return fJetPtMin; }
+  virtual Float_t  GetJetPtFractionMin() const { return fJetPtFractionMin; }
+  virtual Int_t    GetNMatchJets() const { return fNMatchJets; }
+  //virtual Float_t  GetJetDeltaEta() const { return fJetDeltaEta; }
+  //virtual Float_t  GetJetDeltaPhi() const { return fJetDeltaPhi; }
 
   virtual void     SetBranchNames(const TString &branch1, const TString &branch2);
   virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+  virtual void     SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
   virtual void     SetVtxZMin(Float_t z) { fVtxZMin = z; }
   virtual void     SetVtxZMax(Float_t z) { fVtxZMax = z; }
   virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
@@ -44,8 +50,11 @@ class AliAnalysisTaskJetResponse : public AliAnalysisTaskSE {
   virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
   virtual void     SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; }
   virtual void     SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; }
-  virtual void     SetJetDeltaEta(Float_t eta) { fJetDeltaEta = eta; }
-  virtual void     SetJetDeltaPhi(Float_t phi) { fJetDeltaPhi = phi; }
+  virtual void     SetJetPtMin(Float_t pt) { fJetPtMin = pt; }
+  virtual void     SetJetPtFractionMin(Float_t pt) { fJetPtFractionMin = pt; }
+  virtual void     SetNMatchJets(Int_t n) { fNMatchJets = n; }
+  //virtual void     SetJetDeltaEta(Float_t eta) { fJetDeltaEta = eta; }
+  //virtual void     SetJetDeltaPhi(Float_t phi) { fJetDeltaPhi = phi; }
 
  private:
   // ESD/AOD events
@@ -67,30 +76,39 @@ class AliAnalysisTaskJetResponse : public AliAnalysisTaskSE {
   Float_t fCentMax;      // upper bound on centrality
   Float_t fJetEtaMin;     // lower bound on eta for found jets
   Float_t fJetEtaMax;     // upper bound on eta for found jets
-  Float_t fJetDeltaEta;   // max difference in eta to match leading jets
-  Float_t fJetDeltaPhi;   // max difference in phi to match leading jets
+  Float_t fJetPtMin;      // minimum jet pT
+  Float_t fJetPtFractionMin; // minimum fraction for positiv match of jets
+  Int_t   fNMatchJets;       // maximal nb. of jets taken for matching
+  //Float_t fJetDeltaEta;   // max difference in eta to match leading jets
+  //Float_t fJetDeltaPhi;   // max difference in phi to match leading jets
 
   // output objects
   const Int_t fkNbranches;                   //! number of branches to be read
   const Int_t fkEvtClasses;                  //! number of event classes
   TList *fOutputList;                        //! output data container
-  TH1F *fHistEvtSelection;                   //! event selection statistic
-  TH1F **fHistPtLeadingJet;                  //! pt of leading jet
-  TH2F **fHistEtaPhiLeadingJet;              //! eta-phi of leading jet
-  TH2F **fHistEtaPhiLeadingJetCut;           //! eta-phi of leading jet in eta acceptace per event class
-  TH2F  *fHistDeltaEtaDeltaPhiLeadingJet;    //! delta eta vs. delta phi of leading jets
-  TH2F **fHistDeltaEtaEtaLeadingJet;         //! delta eta vs. eta of leading jet per event class
-  TH2F **fHistDeltaPtEtaLeadingJet;          //! delta eta vs. eta of leading jet per event class
+  TH1I  *fHistEvtSelection;                  //! event selection statistic
+  TH1I  *fHistEvtClass;                      //! event classes (from helper task)
+  TH1F  *fHistCentrality;                    //! centrality of the event  
+  TH1F **fHistPtJet;                         //! pt distribution of jets
+  TH2F **fHistEtaPhiJet;                     //! eta-phi distribution of jets (before acceptance cuts)
+  TH2F **fHistEtaPhiJetCut;                  //! eta-phi distribution of jets in eta acceptace per event class
+  TH2F **fHistDeltaEtaDeltaPhiJet;           //! delta eta vs. delta phi of matched jets (before acceptance cuts)
+  TH2F **fHistDeltaEtaDeltaPhiJetCut;        //! delta eta vs. delta phi of matched jets
+  TH2F **fHistDeltaEtaDeltaPhiJetNOMatching; //! delta eta vs. delta phi of jets which do not match
+  TH2F **fHistDeltaEtaEtaJet;                //! delta eta vs. eta of matched jets per event class
+  TH2F **fHistDeltaPtEtaJet;                 //! delta eta vs. eta of matched jets per event class
+  TH2F **fHistPtFraction;                    //! fraction from embedded jet in reconstructed jet per event class
   TH2F  *fHistPtPtExtra;                     //! jet pt response
-  TH2F **fHistPtResponse;                    //! histograms per event class
+  TH2F **fHistPtResponse;                    //! jet pt response per event class
   TH2F **fHistPtSmearing;                    //! emb-jet pt vs (emb+UE - emb) pt
-  TH2F **fHistdR;                            //! shift dR of jets vs (emb+UE - emb) pt
+  TH2F **fHistDeltaR;                        //! shift dR of jets vs (emb+UE - emb) pt
   TH2F **fHistArea;                          //! area of jets vs (emb+UE - emb) pt
+  TH2F **fHistDeltaArea;                     //! delta area of jets vs (emb+UE - emb) pt
 
   AliAnalysisTaskJetResponse(const AliAnalysisTaskJetResponse&); // not implemented
   AliAnalysisTaskJetResponse& operator=(const AliAnalysisTaskJetResponse&); // not implemented
 
-  ClassDef(AliAnalysisTaskJetResponse, 1);
+  ClassDef(AliAnalysisTaskJetResponse, 2);
 };
 
 #endif
index 28ee916..86a5ef7 100644 (file)
@@ -19,12 +19,12 @@ AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(){
 
     // ## set ranges for toy ##
     //SetToyTrackRanges(
-    Double_t minPt = 30.;   Double_t maxPt = 300.;
+    Double_t minPt = 20.;   Double_t maxPt = 200.;
     Double_t minEta = -0.5; Double_t maxEta = 0.5;
     Double_t minPhi = 0.;   Double_t maxPhi = 2*TMath::Pi();
     //fToyDistributionTrackPt: 0 = uniform distribution
     //                         else = exponential / power law (not implemented yet)
-    //task->SetToyNumberOfTrackRange(5,700);
+    //task->SetToyNumberOfTrackRange(4,4);
     //task->SetToyTrackRanges(0.15, 300., 5,-.9, .9, 0., 2*TMath::Pi());
     task->SetToyTrackRanges(minPt,maxPt,0.,minEta,maxEta,minPhi,maxPhi);
     task->SetToyFilterMap((1<<32)-1);
@@ -74,6 +74,7 @@ AliAnalysisTaskFastEmbedding* AddTaskFastEmbedding(const char* filepath, Int_t m
           task->SetAODPath(filepath);
        }
        if(mode==1){ // path to text file with list of paths of multiple AODs
+           Printf("Read aod paths from file %s", filepath);
            TObjArray* array = new TObjArray();
            TObjString* ostr = 0;
            TString line;
index b034a92..3a1ba8f 100644 (file)
@@ -1,10 +1,10 @@
-AliAnalysisTaskJetResponse* AddTaskJetResponse(Char_t* type = "clusters", Char_t* jf = "FASTKT", Float_t radius = 0.4, UInt_t filterMask = 256 , Float_t ptTrackMin = 0.15, Int_t iBack = 1, Int_t eventClassMin = 1, Int_t eventClassMax = 5){
+AliAnalysisTaskJetResponse* AddTaskJetResponse(Char_t* type = "clusters", Char_t* jf = "FASTKT", Float_t radius = 0.4, UInt_t filterMask = 256 , Float_t ptTrackMin = 0.15, Int_t iBack = 1, Int_t eventClassMin = 1, Int_t eventClassMax = 4){
 
   Printf("adding task jet response\n");
 
     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
     if(!mgr){
-       ::Error("AddTaskJetResponse", "No analysis manager to connect ot.");
+       ::Error("AddTaskJetResponse", "No analysis manager to connect to.");
        return NULL;
     }
     if(!mgr->GetInputEventHandler()){
@@ -43,6 +43,8 @@ AliAnalysisTaskJetResponse* AddTaskJetResponse(Char_t* type = "clusters", Char_t
     task->SetEvtClassMax(eventClassMax);
     task->SetCentMin(0.);
     task->SetCentMax(100.);
+    //task->SetJetDeltaEta(0.2);
+    //task->SetJetDeltaPhi(0.2);
 
 
     mgr->AddTask(task);
index d077438..26d4a21 100644 (file)
@@ -593,7 +593,7 @@ void AnalysisTrainPWG4Jets(const char *analysis_mode="local",
 
 
        if(kDeltaAODJetName.Length()==0&&kFilterAOD){
-        if(kIsPbPb)taskCl->SetJetTriggerPtCut(40.);
+        if(kIsPbPb)taskCl->SetJetTriggerPtCut(0.0001);
         else taskCl->SetJetTriggerPtCut(20.);
        }
 
@@ -722,7 +722,6 @@ void AnalysisTrainPWG4Jets(const char *analysis_mode="local",
 
        taskjetServ->SetZVertexCut(8.);
      }
-       taskjetServ->SetDebugLevel(3);
      if(iAODanalysis){
        //       taskjetServ->SetDebugLevel(3);
        taskjetServ->SetAODInput(kTRUE);