X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliAnalysisTaskJetCluster.cxx;h=8b350289fd14357e5b7b1ad881ca7db256683407;hb=46465e39caa4a61b4b739532745c982016c1d230;hp=4cbca532054fad4a0a02209b28a01ea49a704eee;hpb=e2ca75196931e8169f577df9c3fc1a8b7ec282dd;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index 4cbca532054..8b350289fd1 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -70,6 +70,17 @@ ClassImp(AliAnalysisTaskJetCluster) AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){ delete fRef; delete fRandom; + + if(fTCAJetsOut)fTCAJetsOut->Delete(); + delete fTCAJetsOut; + if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); + delete fTCAJetsOutRan; + if(fTCARandomConesOut)fTCARandomConesOut->Delete(); + delete fTCARandomConesOut; + if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); + delete fTCARandomConesOutRan; + if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); + delete fAODJetBackgroundOut; } AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), @@ -84,9 +95,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), fNSkipLeadingRan(0), - fNRandomCones(5), + fNRandomCones(0), fAvgTrials(1), - fExternalWeight(1), + fExternalWeight(1), + fTrackEtaWindow(0.9), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(1), @@ -104,6 +116,11 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fGhostArea(0.01), fActiveAreaRepeats(1), fGhostEtamax(1.5), + fTCAJetsOut(0x0), + fTCAJetsOutRan(0x0), + fTCARandomConesOut(0x0), + fTCARandomConesOutRan(0x0), + fAODJetBackgroundOut(0x0), fRandom(0), fh1Xsec(0x0), fh1Trials(0x0), @@ -186,9 +203,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), fNSkipLeadingRan(0), - fNRandomCones(5), + fNRandomCones(0), fAvgTrials(1), fExternalWeight(1), + fTrackEtaWindow(0.9), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(1), @@ -206,6 +224,11 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fGhostArea(0.01), fActiveAreaRepeats(1), fGhostEtamax(1.5), + fTCAJetsOut(0x0), + fTCAJetsOutRan(0x0), + fTCARandomConesOut(0x0), + fTCARandomConesOutRan(0x0), + fAODJetBackgroundOut(0x0), fRandom(0), fh1Xsec(0x0), fh1Trials(0x0), @@ -310,67 +333,56 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() // Create a new branch for jets... // -> cleared in the UserExec.... // here we can also have the case that the brnaches are written to a separate file - - TClonesArray *tca = new TClonesArray("AliAODJet", 0); - tca->SetName(fNonStdBranch.Data()); - AddAODBranch("TClonesArray",&tca,fNonStdFile.Data()); - + + fTCAJetsOut = new TClonesArray("AliAODJet", 0); + fTCAJetsOut->SetName(fNonStdBranch.Data()); + AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); + - TClonesArray *tcaran = new TClonesArray("AliAODJet", 0); - tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data()); - - + + 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()))){ - AliAODJetEventBackground* evBkg = new AliAODJetEventBackground(); - evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); - AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data()); + fAODJetBackgroundOut = new AliAODJetEventBackground(); + fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); + AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data()); } - // create the branch for the random cones with the same R - TString cName = Form("%sRandomCone",fNonStdBranch.Data()); + } + // create the branch for the random cones with the same R + TString cName = Form("%sRandomCone",fNonStdBranch.Data()); + + if(fNRandomCones>0){ if(!AODEvent()->FindListObject(cName.Data())){ - TClonesArray *tcaC = new TClonesArray("AliAODJet", 0); - tcaC->SetName(cName.Data()); - AddAODBranch("TClonesArray",&tcaC,fNonStdFile.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())){ - TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0); - tcaCran->SetName(cName.Data()); - AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data()); + fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0); + fTCARandomConesOutRan->SetName(cName.Data()); + AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data()); } } - + if(fNonStdFile.Length()!=0){ // // case that we have an AOD extension we need to fetch the jets from the extended output - // we identifay the extension aod event by looking for the branchname + // we identify the extension aod event by looking for the branchname AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); - TObjArray* extArray = aodH->GetExtensions(); - if (extArray) { - TIter next(extArray); - while ((fAODExtension=(AliAODExtension*)next())){ - TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()); - if(fDebug>10){ - Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__); - fAODExtension->GetAOD()->Dump(); - } - if(obj){ - if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data()); - break; - } - fAODExtension = 0; - } - } + // case that we have an AOD extension we need can fetch the background maybe from the extended output + fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0); } } if(!fHistList)fHistList = new TList(); fHistList->SetOwner(); + PostData(1, fHistList); // post data in any case once Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); @@ -441,9 +453,9 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt); - fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt); - fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt); - fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt); + fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); + fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); + fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt); fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5); fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5); @@ -513,12 +525,13 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); - if(fUseBackgroundCalc){ + if(fNRandomCones>0&&fUseBackgroundCalc){ for(int i = 0;i<3;i++){ fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100); } } + for(int i = 0;i < kMaxCent;i++){ fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1)); fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1)); @@ -552,7 +565,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh1Z); fHistList->Add(fh1ZSelect); fHistList->Add(fh1ZPhySel); - if(fUseBackgroundCalc){ + if(fNRandomCones>0&&fUseBackgroundCalc){ for(int i = 0;i<3;i++){ fHistList->Add(fh1BiARandomCones[i]); fHistList->Add(fh1BiARandomConesRan[i]); @@ -592,7 +605,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh2NConstLeadingPtRan); fHistList->Add(fh2TracksLeadingJetPhiPtRan); fHistList->Add(fh2TracksLeadingJetPhiPtWRan); - } + } // =========== Switch on Sumw2 for all histos =========== for (Int_t i=0; iGetEntries(); ++i) { @@ -630,40 +643,18 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // handle and reset the output jet branch - // only need this once - TClonesArray* jarray = 0; - TClonesArray* jarrayran = 0; - TClonesArray* rConeArray = 0; - TClonesArray* rConeArrayRan = 0; - AliAODJetEventBackground* evBkg = 0; - if(fNonStdBranch.Length()!=0) { - if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data())); - if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data())); - if(jarray)jarray->Delete(); // this is our responsibility, clear before filling again - if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random"))); - if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random"))); - if(jarrayran)jarrayran->Delete(); // this is our responsibility, clear before filling again - - if(fUseBackgroundCalc){ - if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))); - if(!evBkg) evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))); - if(evBkg)evBkg->Reset(); - - TString cName = Form("%sRandomCone",fNonStdBranch.Data()); - if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data())); - if(!rConeArray)rConeArray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data())); - if(rConeArray)rConeArray->Delete(); - - cName = Form("%sRandomCone_random",fNonStdBranch.Data()); - if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data())); - if(!rConeArrayRan)rConeArrayRan = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data())); - if(rConeArrayRan)rConeArrayRan->Delete(); - } - } - static AliAODJetEventBackground* externalBackground = 0; + if(fTCAJetsOut)fTCAJetsOut->Delete(); + if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); + if(fTCARandomConesOut)fTCARandomConesOut->Delete(); + if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); + if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); + + AliAODJetEventBackground* externalBackground = 0; if(!externalBackground&&fBackgroundBranch.Length()){ externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data())); + if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); + if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());; } // // Execute analysis for current event @@ -715,7 +706,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Float_t yvtx = vtxAOD->GetY(); Float_t xvtx = vtxAOD->GetX(); Float_t r2 = yvtx*yvtx+xvtx*xvtx; - if(TMath::Abs(zVtx)<20.&&r2<1.){ // apply vertex cut later on + if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on if(physicsSelection){ selectEvent = true; } @@ -759,31 +750,6 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) vector inputParticlesRec; vector inputParticlesRecRan; - - - // create a random jet within the acceptance - - if(fUseBackgroundCalc){ - Double_t etaMax = 0.9 - fRparam; - Int_t nCone = 0; - Int_t nConeRan = 0; - Double_t pT = 1; // small number - for(int ir = 0;ir < fNRandomCones;ir++){ - Double_t eta = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax - Double_t phi = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi - // massless jet - Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); - Double_t pZ = pT/TMath::Tan(theta); - Double_t pX = pT * TMath::Cos(phi); - Double_t pY = pT * TMath::Sin(phi); - Double_t p = TMath::Sqrt(pT*pT+pZ*pZ); - AliAODJet tmpRec (pX,pY,pZ, p); - tmpRec.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below - if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRec); - if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRec); - } - } - // Generate the random cones @@ -796,19 +762,10 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) jInp.set_user_index(i); inputParticlesRec.push_back(jInp); - if(fUseBackgroundCalc&&rConeArray){ - for(int ir = 0;ir < fNRandomCones;ir++){ - AliAODJet *jC = (AliAODJet*)rConeArray->At(ir); - if(jC&&jC->DeltaR(vp)SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); - } - } - } - // the randomized input changes eta and phi, but keeps the p_T if(i>=fNSkipLeadingRan){// eventually skip the leading particles Double_t pT = vp->Pt(); - Double_t eta = 1.8 * fRandom->Rndm() - 0.9; + Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow; Double_t phi = 2.* TMath::Pi() * fRandom->Rndm(); Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); @@ -822,19 +779,10 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) jInpRan.set_user_index(i); inputParticlesRecRan.push_back(jInpRan); vTmpRan.SetPxPyPzE(pX,pY,pZ,p); - - if(fUseBackgroundCalc&&rConeArrayRan){ - for(int ir = 0;ir < fNRandomCones;ir++){ - AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir); - if(jC&&jC->DeltaR(&vTmpRan)SetBgEnergy(jC->ChargedBgEnergy()+vTmpRan.Pt(),0); - } - } - } } // fill the tref array, only needed when we write out jets - if(jarray){ + if(fTCAJetsOut){ if(i == 0){ fRef->Delete(); // make sure to delete before placement new... new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); @@ -842,47 +790,6 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fRef->Add(vp); } }// recparticles - if(fUseBackgroundCalc){ - Float_t jetArea = fRparam*fRparam*TMath::Pi(); - for(int ir = 0;ir < fNRandomCones;ir++){ - // rescale the momntum vectors for the random cones - if(!rConeArray)continue; - AliAODJet *rC = (AliAODJet*)rConeArray->At(ir); - if(rC){ - Double_t eta = rC->Eta(); - Double_t phi = rC->Phi(); - // massless jet, unit vector - Double_t pT = rC->ChargedBgEnergy(); - if(pT<=0)pT = 0.1; // for almost empty events - Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); - Double_t pZ = pT/TMath::Tan(theta); - Double_t pX = pT * TMath::Cos(phi); - Double_t pY = pT * TMath::Sin(phi); - Double_t p = TMath::Sqrt(pT*pT+pZ*pZ); - rC->SetPxPyPzE(pX,pY,pZ, p); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - rC = (AliAODJet*)rConeArrayRan->At(ir); - // same wit random - if(rC){ - Double_t eta = rC->Eta(); - Double_t phi = rC->Phi(); - // massless jet, unit vector - Double_t pT = rC->ChargedBgEnergy(); - if(pT<=0)pT = 0.1;// for almost empty events - Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta)); - Double_t pZ = pT/TMath::Tan(theta); - Double_t pX = pT * TMath::Cos(phi); - Double_t pY = pT * TMath::Sin(phi); - Double_t p = TMath::Sqrt(pT*pT+pZ*pZ); - rC->SetPxPyPzE(pX,pY,pZ, p); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - } - } - if(inputParticlesRec.size()==0){ if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); @@ -956,7 +863,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(externalBackground){ // carefull has to be filled in a task before // todo, ReArrange to the botom - pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); + pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged(); } pt = leadingJet.Pt() - pTback; // correlation of leading jet with tracks @@ -981,55 +888,56 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } - - for(int j = 0; j < nRec;j++){ - AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E()); - aodOutJet = 0; - nAodOutTracks = 0; - Float_t tmpPt = tmpRec.Pt(); - Float_t tmpPtBack = 0; - if(externalBackground){ - // carefull has to be filled in a task before + TLorentzVector vecareab; + for(int j = 0; j < nRec;j++){ + AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E()); + aodOutJet = 0; + nAodOutTracks = 0; + Float_t tmpPt = tmpRec.Pt(); + + if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted... + aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec); + aodOutJet->GetRefTracks()->Clear(); + Double_t area1 = clustSeq.area(sortedJets[j]); + aodOutJet->SetEffArea(area1,0); + fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]); + vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e()); + aodOutJet->SetVectorAreaCharged(&vecareab); + } + + + Float_t tmpPtBack = 0; + if(externalBackground){ + // carefull has to be filled in a task before // todo, ReArrange to the botom - tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged(); - } - tmpPt = tmpPt - tmpPtBack; - if(tmpPt<0)tmpPt = 0; // avoid negative weights... - - fh1PtJetsRecIn->Fill(tmpPt); - // Fill Spectra with constituents - vector constituents = clustSeq.constituents(sortedJets[j]); - fh1NConstRec->Fill(constituents.size()); - fh2PtNch->Fill(nCh,tmpPt); - fh2PtNchN->Fill(nCh,tmpPt,constituents.size()); - fh2NConstPt->Fill(tmpPt,constituents.size()); - // loop over constiutents and fill spectrum - - // Add the jet information and the track references to the output AOD - - if(tmpPt>fJetOutputMinPt&&jarray){ - aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec); - Double_t area1 = clustSeq.area(sortedJets[j]); - aodOutJet->SetEffArea(area1,0); - } + tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged(); + } + tmpPt = tmpPt - tmpPtBack; + if(tmpPt<0)tmpPt = 0; // avoid negative weights... + + fh1PtJetsRecIn->Fill(tmpPt); + // Fill Spectra with constituents + vector constituents = clustSeq.constituents(sortedJets[j]); - - for(unsigned int ic = 0; ic < constituents.size();ic++){ - AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); - fh1PtJetConstRec->Fill(part->Pt()); - if(aodOutJet){ - aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); - } - if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt()); - } - + fh1NConstRec->Fill(constituents.size()); + fh2PtNch->Fill(nCh,tmpPt); + fh2PtNchN->Fill(nCh,tmpPt,constituents.size()); + fh2NConstPt->Fill(tmpPt,constituents.size()); + // loop over constiutents and fill spectrum + + for(unsigned int ic = 0; ic < constituents.size();ic++){ + AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); + fh1PtJetConstRec->Fill(part->Pt()); + if(aodOutJet){ + aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); + } + if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt()); + } + // correlation Float_t tmpPhi = tmpRec.Phi(); Float_t tmpEta = tmpRec.Eta(); - if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; - - - + if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; if(j==0){ fh1PtJetsLeadingRecIn->Fill(tmpPt); fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); @@ -1051,12 +959,145 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); } fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); - } + }// loop over reconstructed jets delete recIter; + + + + // Add the random cones... + if(fNRandomCones>0&&fTCARandomConesOut){ + // create a random jet within the acceptance + Double_t etaMax = fTrackEtaWindow - fRparam; + Int_t nCone = 0; + Int_t nConeRan = 0; + Double_t pTC = 1; // small number + for(int ir = 0;ir < fNRandomCones;ir++){ + Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax + Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi + // massless jet + Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); + Double_t pZC = pTC/TMath::Tan(thetaC); + Double_t pXC = pTC * TMath::Cos(phiC); + Double_t pYC = pTC * TMath::Sin(phiC); + 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 + AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E()); + if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){ + skip = true; + break; + } + } + // test for overlap with previous cones to avoid double counting + for(int iic = 0;iicAt(iic); + if(iicone){ + if(iicone->DeltaR(&tmpRecC)<2.*fRparam){ + skip = true; + break; + } + } + } + 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); + }// loop over random cones creation + + + // loop over the reconstructed particles and add up the pT in the random cones + // maybe better to loop over randomized particles not in the real jets... + // but this by definition brings dow average energy in the whole event + AliAODJet vTmpRanR(1,0,0,1); + for(int i = 0; i < recParticles.GetEntries(); i++){ + AliVParticle *vp = (AliVParticle*)recParticles.At(i); + if(fTCARandomConesOut){ + for(int ir = 0;ir < fNRandomCones;ir++){ + AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir); + if(jC&&jC->DeltaR(vp)SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); + } + } + }// add up energy in cone + + // the randomized input changes eta and phi, but keeps the p_T + if(i>=fNSkipLeadingRan){// eventually skip the leading particles + Double_t pTR = vp->Pt(); + Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow; + Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm(); + + Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR)); + Double_t pZR = pTR/TMath::Tan(thetaR); + + Double_t pXR = pTR * TMath::Cos(phiR); + Double_t pYR = pTR * TMath::Sin(phiR); + Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR); + vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR); + if(fTCARandomConesOutRan){ + for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ + AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir); + if(jC&&jC->DeltaR(&vTmpRanR)SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); + } + } + } + } + }// loop over recparticles + + Float_t jetArea = fRparam*fRparam*TMath::Pi(); + if(fTCARandomConesOut){ + for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ + // rescale the momntum vectors for the random cones + + AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir); + if(rC){ + Double_t etaC = rC->Eta(); + Double_t phiC = rC->Phi(); + // massless jet, unit vector + pTC = rC->ChargedBgEnergy(); + if(pTC<=0)pTC = 0.1; // 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); + Double_t pYC = pTC * TMath::Sin(phiC); + Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); + rC->SetPxPyPzE(pXC,pYC,pZC, pC); + rC->SetBgEnergy(0,0); + rC->SetEffArea(jetArea,0); + } + } + } + if(!fTCARandomConesOutRan){ + for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ + AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir); + // same wit random + if(rC){ + Double_t etaC = rC->Eta(); + Double_t phiC = rC->Phi(); + // massless jet, unit vector + pTC = rC->ChargedBgEnergy(); + if(pTC<=0)pTC = 0.1;// 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); + Double_t pYC = pTC * TMath::Sin(phiC); + Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); + rC->SetPxPyPzE(pXC,pYC,pZC, pC); + rC->SetBgEnergy(0,0); + rC->SetEffArea(jetArea,0); + } + } + } + }// if(fNRandomCones //background estimates:all bckg jets(0) & wo the 2 hardest(1) - if(evBkg){ + + + + + if(fAODJetBackgroundOut){ vector jets2=sortedJets; if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); Double_t bkg1=0; @@ -1066,14 +1107,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Double_t sigma2=0.; Double_t meanarea2=0.; - clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true); - evBkg->SetBackground(0,bkg1,sigma1,meanarea1); + clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true); + fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1); // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone)); // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone)); clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); - evBkg->SetBackground(1,bkg2,sigma2,meanarea2); + fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); @@ -1153,11 +1194,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Int_t nRecOverRan = inclusiveJetsRan.size(); Int_t nRecRan = inclusiveJetsRan.size(); + if(inclusiveJetsRan.size()>0){ AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E()); Float_t pt = leadingJet.Pt(); Int_t iCount = 0; + TLorentzVector vecarearanb; + for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); while(ptFill(nCh,tmpPt,constituents.size()); - if(tmpPt>fJetOutputMinPt&&jarrayran){ - aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec); + if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ + aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); - - aodOutJetRan->SetEffArea(arearan,0); } + aodOutJetRan->GetRefTracks()->Clear(); + aodOutJetRan->SetEffArea(arearan,0); + fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]); + vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e()); + aodOutJetRan->SetVectorAreaCharged(&vecarearanb); - - - - for(unsigned int ic = 0; ic < constituents.size();ic++){ - // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); - // fh1PtJetConstRec->Fill(part->Pt()); - if(aodOutJetRan){ - aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index())); - } } - - - // correlation Float_t tmpPhi = tmpRec.Phi(); @@ -1236,12 +1271,12 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } - if(evBkg){ + if(fAODJetBackgroundOut){ Double_t bkg3=0.; Double_t sigma3=0.; Double_t meanarea3=0.; clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); - evBkg->SetBackground(2,bkg3,sigma3,meanarea3); + fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3); // float areaRandomCone = rRandomCone2 *TMath::Pi(); /* fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); @@ -1268,9 +1303,9 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) */ float rho = 0; if(externalBackground)rho = externalBackground->GetBackground(2); - if(jarray){ - for(int i = 0;i < jarray->GetEntriesFast();i++){ - AliAODJet *jet = (AliAODJet*)jarray->At(i); + if(fTCAJetsOut){ + for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){ + AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i); Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged(); if(ptSub>=minPt){ select = true; @@ -1304,20 +1339,43 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); Int_t iCount = 0; - if(type==kTrackAOD){ - AliAODEvent *aod = 0; - if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); - else aod = AODEvent(); - if(!aod){ - return iCount; + if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){ + if(type!=kTrackAODextraonly) { + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod){ + 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()Add(tr); + 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())>0.9)continue; - if(tr->Pt()Add(tr); - iCount++; + if(type==kTrackAODextra || type==kTrackAODextraonly) { + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + + if(!aod){ + return iCount; + } + TClonesArray *aodExtraTracks = dynamic_cast(aod->FindListObject("aodExtraTracks")); + if(!aodExtraTracks)return iCount; + for(int it =0; itGetEntries(); it++) { + AliVParticle *track = dynamic_cast ((*aodExtraTracks)[it]); + if (!track) continue; + + AliAODTrack *trackAOD = dynamic_cast (track); + if(!trackAOD)continue; + if(trackAOD->Pt()Add(trackAOD); + iCount++; + } } } else if (type == kTrackKineAll||type == kTrackKineCharged){ @@ -1362,7 +1420,7 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ list->Add(part); } else { - if(TMath::Abs(part->Eta())>0.9)continue; + if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue; list->Add(part); } iCount++;