]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliAnalysisTaskJetCluster.cxx
Add print on the found particle tag
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
index 2265781b56778070b9c6a987a6dfeb5a537b8bb9..1b27611e03ad0eaef7c61bdf5cb3306359db1308 100644 (file)
@@ -84,26 +84,33 @@ AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
   delete fAODJetBackgroundOut;
 }
 
-AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
+AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): 
+  AliAnalysisTaskSE(),
   fAOD(0x0),
   fAODExtension(0x0),
   fRef(new TRefArray),
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
-  fUseGlobalSelection(kFALSE),
   fUseBackgroundCalc(kFALSE),
+  fEventSelection(kFALSE),     
   fFilterMask(0),
+  fFilterType(0),
+  fJetTypes(kJet),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),  
   fNSkipLeadingRan(0),
+  fNSkipLeadingCone(0),
   fNRandomCones(0),
   fAvgTrials(1),
   fExternalWeight(1),
   fTrackEtaWindow(0.9),    
   fRecEtaWindow(0.5),
   fTrackPtCut(0.),                                                     
-  fJetOutputMinPt(1),
+  fJetOutputMinPt(0.150),
+  fMaxTrackPtInJet(100.),
   fJetTriggerPtCut(0),
+  fVtxZCut(8),
+  fVtxR2Cut(1),
   fCentCutUp(0),
   fCentCutLo(0),
   fNonStdBranch(""),
@@ -198,20 +205,26 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fRef(new TRefArray),
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
-  fUseGlobalSelection(kFALSE),
   fUseBackgroundCalc(kFALSE),
+  fEventSelection(kFALSE),                                                     
   fFilterMask(0),
+  fFilterType(0),
+  fJetTypes(kJet),
   fTrackTypeRec(kTrackUndef),
   fTrackTypeGen(kTrackUndef),
   fNSkipLeadingRan(0),
+  fNSkipLeadingCone(0),
   fNRandomCones(0),
   fAvgTrials(1),
   fExternalWeight(1),    
   fTrackEtaWindow(0.9),    
   fRecEtaWindow(0.5),
   fTrackPtCut(0.),                                                     
-  fJetOutputMinPt(1),
+  fJetOutputMinPt(0.150),
+  fMaxTrackPtInJet(100.),
   fJetTriggerPtCut(0),
+  fVtxZCut(8),
+  fVtxR2Cut(1),
   fCentCutUp(0),
   fCentCutLo(0),
   fNonStdBranch(""),
@@ -335,16 +348,18 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
       //  -> cleared in the UserExec....
       // here we can also have the case that the brnaches are written to a separate file
       
-      fTCAJetsOut = new TClonesArray("AliAODJet", 0);
-      fTCAJetsOut->SetName(fNonStdBranch.Data());
-      AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
-      
+      if(fJetTypes&kJet){
+       fTCAJetsOut = new TClonesArray("AliAODJet", 0);
+       fTCAJetsOut->SetName(fNonStdBranch.Data());
+       AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
+      }
    
-      
-      fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
-      fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
-      AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
-      
+      if(fJetTypes&kJetRan){
+       fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
+       fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
+       AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
+      }
+
       if(fUseBackgroundCalc){
        if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
          fAODJetBackgroundOut = new AliAODJetEventBackground();
@@ -353,20 +368,24 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
        }
       }
       // create the branch for the random cones with the same R 
-      TString cName = Form("%sRandomCone",fNonStdBranch.Data());
-
+      TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
+  
       if(fNRandomCones>0){
-       if(!AODEvent()->FindListObject(cName.Data())){
-         fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
-         fTCARandomConesOut->SetName(cName.Data());
-         AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
+       if(fJetTypes&kRC){
+         if(!AODEvent()->FindListObject(cName.Data())){
+           fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
+           fTCARandomConesOut->SetName(cName.Data());
+           AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
+         }
        }
        // create the branch with the random for the random cones on the random event
-       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
-       if(!AODEvent()->FindListObject(cName.Data())){
-         fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
-         fTCARandomConesOutRan->SetName(cName.Data());
-         AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
+       if(fJetTypes&kRCRan){
+         cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+         if(!AODEvent()->FindListObject(cName.Data())){
+           fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
+           fTCARandomConesOutRan->SetName(cName.Data());
+           AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
+         }
        }
       }
     
@@ -391,14 +410,14 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
   //
   //  Histogram
     
-  const Int_t nBinPt = 200;
+  const Int_t nBinPt = 100;
   Double_t binLimitsPt[nBinPt+1];
   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
     if(iPt == 0){
       binLimitsPt[iPt] = 0.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
     }
   }
   
@@ -426,7 +445,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
     }
   }
 
-  const int nChMax = 4000;
+  const int nChMax = 5000;
 
   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
@@ -634,15 +653,6 @@ void AliAnalysisTaskJetCluster::Init()
 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 {
 
-  if(fUseGlobalSelection){
-    // no selection by the service task, we continue
-    if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
-    PostData(1, fHistList);
-    return;
-  }
-
-
-
   // handle and reset the output jet branch 
 
   if(fTCAJetsOut)fTCAJetsOut->Delete();
@@ -702,24 +712,28 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       fh1ZPhySel->Fill(zVtx);
     }
 
-
-    if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
+    if(fEventSelection){
+      if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
        Float_t yvtx = vtxAOD->GetY();
        Float_t xvtx = vtxAOD->GetX();
        Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
-       if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
+       if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
          if(physicsSelection){
            selectEvent = true;
          }
        }
-    }
-    if(fCentCutUp>0){
-      if(cent<fCentCutLo||cent>fCentCutUp){
-       selectEvent = false;
       }
+      if(fCentCutUp>0){
+       if(cent<fCentCutLo||cent>fCentCutUp){
+         selectEvent = false;
+       }
+      }
+    }else{
+      selectEvent = true;
     }
-
   }
+  
+
   if(!selectEvent){
     PostData(1, fHistList);
     return;
@@ -803,17 +817,18 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
  
   // now create the object that holds info about ghosts                        
+  /*
   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
     // reduce CPU time...
     fGhostArea = 0.5; 
     fActiveAreaRepeats = 0; 
   }
-   
- fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
+  */
+
+  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
   fastjet::AreaType areaType =   fastjet::active_area;
   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
-  vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
   
   //range where to compute background
@@ -825,9 +840,9 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
  
 
+  const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
+  const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
 
-  inclusiveJets = clustSeq.inclusive_jets();
-  sortedJets    = sorted_by_pt(inclusiveJets);
  
   fh1NJetsRec->Fill(sortedJets.size());
 
@@ -917,8 +932,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
       
       fh1PtJetsRecIn->Fill(tmpPt);
-      // Fill Spectra with constituents
-      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
+      // Fill Spectra with constituentsemacs
+      const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
 
       fh1NConstRec->Fill(constituents.size());
       fh2PtNch->Fill(nCh,tmpPt);
@@ -931,6 +946,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        fh1PtJetConstRec->Fill(part->Pt());
        if(aodOutJet){
          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
+         if(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
        }
        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
       }
@@ -983,7 +999,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
        AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
        bool skip = false;
-       for(int jj = 0; jj < TMath::Min(nRec,2);jj++){// test for overlap with leading jets
+       for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
         AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
         if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
           skip = true;
@@ -1002,8 +1018,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        }
        if(skip)continue;
        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
-       if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nConeRan++]) AliAODJet(tmpRecC);
-       if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nCone++]) AliAODJet(tmpRecC);
+       if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
+       if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
      }// loop over random cones creation
 
   
@@ -1017,6 +1033,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
         for(int ir = 0;ir < fNRandomCones;ir++){
           AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);  
           if(jC&&jC->DeltaR(vp)<fRparam){
+            if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
             jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
           }
         }  
@@ -1039,6 +1056,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
           for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
             AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);  
             if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
+              if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
               jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
             }
           }  
@@ -1057,7 +1075,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
           Double_t phiC = rC->Phi();
           // massless jet, unit vector
           pTC = rC->ChargedBgEnergy();
-          if(pTC<=0)pTC = 0.1; // for almost empty events
+          if(pTC<=0)pTC = 0.001; // for almost empty events
           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
           Double_t pZC = pTC/TMath::Tan(thetaC);
           Double_t pXC = pTC * TMath::Cos(phiC);
@@ -1069,7 +1087,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
         }
        }
      }
-     if(!fTCARandomConesOutRan){
+     if(fTCARandomConesOutRan){
        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
         AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
         // same wit random
@@ -1078,7 +1096,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
           Double_t phiC = rC->Phi();
           // massless jet, unit vector
           pTC = rC->ChargedBgEnergy();
-          if(pTC<=0)pTC = 0.1;// for almost empty events
+          if(pTC<=0)pTC = 0.001;// for almost empty events
           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
           Double_t pZC = pTC/TMath::Tan(thetaC);
           Double_t pXC = pTC * TMath::Cos(phiC);
@@ -1181,13 +1199,12 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
  }
 
  // find the random jets
- vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
+
  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
   
- inclusiveJetsRan = clustSeqRan.inclusive_jets();
- sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
-
  // fill the jet information from random track
+ const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
+ const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
 
   fh1NJetsRecRan->Fill(sortedJetsRan.size());
 
@@ -1241,7 +1258,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       Float_t tmpPt = tmpRec.Pt();
       fh1PtJetsRecInRan->Fill(tmpPt);
       // Fill Spectra with constituents
-      vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
+      const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
       fh1NConstRecRan->Fill(constituents.size());
       fh2NConstPtRan->Fill(tmpPt,constituents.size());
       fh2PtNchRan->Fill(nCh,tmpPt);
@@ -1314,7 +1331,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        }
      }
    }   
-
    if(select){
      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
      fh1CentralitySelect->Fill(cent);
@@ -1322,8 +1339,12 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
      aodH->SetFillAOD(kTRUE);
    }
  }
-
- if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+ if (fDebug > 2){
+   if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
+   if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
+   if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
+   if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
+ }
  PostData(1, fHistList);
 }
 
@@ -1346,13 +1367,28 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
       else aod = AODEvent();
       if(!aod){
+       if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
        return iCount;
       }
       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())>fTrackEtaWindow)continue;
-       if(tr->Pt()<fTrackPtCut)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)||(!bGood)))){
+         if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
+         continue;
+       }
+       if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
+         if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
+         continue;
+       }
+       if(tr->Pt()<fTrackPtCut){
+         if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
+         continue;
+       }
+       if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
        list->Add(tr);
        iCount++;
       }
@@ -1373,6 +1409,12 @@ Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
 
        AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
        if(!trackAOD)continue;
+       Bool_t bGood = false;
+       if(fFilterType == 0)bGood = true;
+       else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
+       else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
+       if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
+        if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
        if(trackAOD->Pt()<fTrackPtCut) continue;
        list->Add(trackAOD);
        iCount++;