8 #include "AliAnalysisTask.h"
9 #include "AliInputEventHandler.h"
10 #include "AliESDtrack.h"
11 #include "AliAODVertex.h"
12 #include "AliAODCluster.h"
17 #include <TInterpreter.h>
26 #include <TLorentzVector.h>
27 #include <TClonesArray.h>
28 #include <TRefArray.h>
30 #include "TDatabasePDG.h"
31 #include "AliAnalysisManager.h"
32 #include "AliJetFinder.h"
33 #include "AliJetHeader.h"
34 #include "AliJetReader.h"
35 #include "AliJetReaderHeader.h"
36 #include "AliUA1JetHeaderV1.h"
37 #include "AliSISConeJetHeader.h"
38 #include "AliESDEvent.h"
39 #include "AliAODEvent.h"
40 #include "AliAODHandler.h"
41 #include "AliAODInputHandler.h"
42 #include "AliAODTrack.h"
43 #include "AliVParticle.h"
44 #include "AliAODJet.h"
45 #include "AliAODJetEventBackground.h"
46 #include "AliMCParticle.h"
47 #include "AliAODMCParticle.h"
48 #include "AliMCEventHandler.h"
49 #include "AliMCEvent.h"
52 #include "AliAODHeader.h"
53 #include "AliAODMCHeader.h"
54 #include "AliJetKineReaderHeader.h"
55 #include "AliGenCocktailEventHeader.h"
56 #include "AliInputEventHandler.h"
57 #include "AliGenEventHeader.h"
58 #include "AliGenDPMjetEventHeader.h"
60 #include "AliAnalysisTaskCheckSingleTrackJetRejection.h"
61 #include "AliAnalysisTaskPhiCorrelations.h"
62 #include "AliAnalysisHelperJetTasks.h"
63 #include "AliPWG4HighPtQAMC.h"
69 ClassImp(AliAnalysisTaskCheckSingleTrackJetRejection)
71 //________________________________________________________________________
72 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection():
91 fJetRecEtaWindow(0.5),
97 fH2jetMCAKT04_Jetpt_maxpt(0x0),
98 fH2jetAKT04_Jetpt_maxpt (0x0)
101 for(int j=0;j<6;j++){
102 fH1jetMCAKT04_pt [j]=0;
103 fH1jetAKT04_pt [j]=0;
105 fH2jetMCAKT04_Eratio [j]=0;
106 fH1jetMCAKT04_match [j]=0;
109 // Default constructor
112 //________________________________________________________________________
113 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection(const char *name):
114 AliAnalysisTaskSE(name),
115 fUseAODInput(kFALSE),
132 fJetRecEtaWindow(0.5),
137 fH2jetMCAKT04_Jetpt_maxpt(0x0),
138 fH2jetAKT04_Jetpt_maxpt (0x0)
142 for(int j=0;j<6;j++){
143 fH1jetMCAKT04_pt [j]=0;
144 fH1jetAKT04_pt [j]=0;
146 fH2jetMCAKT04_Eratio [j]=0;
147 fH1jetMCAKT04_match [j]=0;
151 // Define input and output slots here
152 // Input slot #0 works with a TChain
153 DefineInput(0, TChain::Class());
154 // Output slot #0 id reserved by the base class for AOD
155 // Output slot #1 writes into a TH1 container
156 DefineOutput(1, TList::Class());
160 //________________________________________________________________________
161 void AliAnalysisTaskCheckSingleTrackJetRejection::UserCreateOutputObjects()
166 if (!fHistList){ fHistList = new TList();fHistList->SetOwner(kTRUE); cout<<"TList is created for output "<<endl;}
168 Bool_t oldStatus = TH1::AddDirectoryStatus();
169 TH1::AddDirectory(kFALSE);
173 fH1Xsec = new TProfile("Xsec","Xsec",1,0,1);
174 fH1Trials = new TH1F ("Trials","Trials",1,0,1);
175 histname = Form("fH2jetMCAKT04_Jetpt_maxpt");
176 fH2jetMCAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
177 histname = Form("fH2jetAKT04_Jetpt_maxpt");
178 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
179 for(int j=0;j<6;j++){
180 histname = Form("fH1jetMCAKT04_pt%d",j);
181 fH1jetMCAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
182 histname = Form("fH2jetAKT04_pt%d",j);
183 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
185 histname = Form("fH2jetMCAKT04_Eratio%d",j);
186 fH2jetMCAKT04_Eratio[j] = new TH2F(histname,histname,200,0,200,100,0,30);
187 histname = Form("fH1jetMCAKT04_match%d",j);
188 fH1jetMCAKT04_match[j] = new TH1F(histname,histname,200,0,200);
190 fHistList->Add(fH1Xsec);
191 fHistList->Add(fH1Trials);
192 fHistList->Add(fH2jetMCAKT04_Jetpt_maxpt);
193 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
194 for(int j=0;j<6;j++){
195 fHistList->Add(fH1jetAKT04_pt[j]);
196 fHistList->Add(fH1jetMCAKT04_pt[j]);
197 fHistList->Add(fH2jetMCAKT04_Eratio[j]);
198 fHistList->Add(fH1jetMCAKT04_match[j]);
202 fH1Events = new TH1F("Events","Number of Events",1,0,1);
203 histname = Form("fH2jetAKT04_Jetpt_maxpt");
204 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
205 for(int j=0;j<6;j++){
206 histname = Form("fH2jetAKT04_pt%d",j);
207 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
209 fHistList->Add(fH1Events);
210 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
211 for(int j=0;j<6;j++){
212 fHistList->Add(fH1jetAKT04_pt[j]);
217 // =========== Switch on Sumw2 for all histos ===========
218 for (Int_t i=0; i<fHistList->GetEntries(); ++i)
220 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
226 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
229 TH1::AddDirectory(oldStatus);
232 PostData(1,fHistList);
236 //----------------------------------------------------------------------
237 void AliAnalysisTaskCheckSingleTrackJetRejection::Init()
240 if (fDebug) printf("AnalysisTaskCheckSingleTrackJetRejection::Init() \n");
244 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::Notify()
246 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
247 fAODOut = AODEvent();
248 if(fNonStdFile.Length()!=0){
249 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
250 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
252 if(fDebug>1)Printf("AODExtension found for %s ",fNonStdFile.Data());
256 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
261 TFile *curfile = tree->GetCurrentFile();
263 Error("Notify","No current file");
268 PythiaInfoFromFile(curfile->GetName(),fxsec,ftrial);
269 //cout<<" Xsec "<<fxsec<<" trial "<<ftrial<<endl;
270 fH1Xsec ->Fill(0.,fxsec);
271 fH1Trials->Fill(0.,ftrial);
274 //Float_t totalEvent;
275 //totalEvent = GetTotalEvents(curfile->GetName());
276 //fH1Events->Fill(0.,totalEvent);
281 printf("Reading File %s ",fInputHandler->GetTree()->GetCurrentFile()->GetName());
285 void AliAnalysisTaskCheckSingleTrackJetRejection::FinishTaskOutput()
291 //________________________________________________________________________
292 void AliAnalysisTaskCheckSingleTrackJetRejection::UserExec(Option_t *)
296 // Main loop (called each event)
297 // Execute analysis for current event
300 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
301 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
302 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
303 if(fUseAODInput) fAODIn = fAOD;
304 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
307 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
308 if( handler && handler->InheritsFrom("AliAODHandler") ) {
309 fAOD = ((AliAODHandler*)handler)->GetAOD();
311 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
315 if(!fAODIn && !fUseAODInput){ // case we have AOD in input & output and want jets from output
316 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
317 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
318 fAODIn = ((AliAODHandler*)outHandler)->GetAOD();
319 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
323 //if(fNonStdFile.Length()!=0){
324 // // case we have an AOD extension - fetch the jets from the extended output
325 // AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
326 // fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
327 // if(!fAODExtension){
328 // if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
331 //fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
333 Printf("ERROR %s : fAODIn not available",(char*)__FILE__);
337 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
338 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
339 if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
340 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
343 if(!IsMC)fH1Events->Fill(0);
347 // start jet analysis
350 Double_t Jet_pt [20][10000];
351 Double_t Jet_eta[20][10000];
352 Double_t Jet_phi[20][10000];
354 for(int i=0;i<20;i++){
356 for(int j=0;j<10000;j++){
363 // get AOD event from input/ouput
365 cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
366 cAdd += Form("B%d",(int)BackM);
367 cAdd += Form("_Filter%05d",Filtermask);
368 cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
369 cAdd += Form("_Skip%02d",SkipCone);
370 TString Branchname_gen,Branchname_rec;
371 //Branchname_gen = Form("clustersMCKINE2_%s%s",JFAlg.Data(),cAdd.Data());
372 Branchname_gen = Form("clustersAODMC2_%s%s",JFAlg.Data(),cAdd.Data());
373 Branchname_rec = Form("clustersAOD_%s%s",JFAlg.Data(),cAdd.Data());
377 double maxpt[2][1000];for(int i=0;i<1000;i++){maxpt[0][i]=0;maxpt[1][i]=0;}
378 int nearrecID[1000]; for(int i=0;i<1000;i++){nearrecID[i]=99999;}
380 for(int algorithm=0;algorithm<2;algorithm++){
381 //for LHC11a1 LHC11a2 official
382 if((!IsMC&&algorithm==0))continue;
384 if(algorithm==0)fJetBranch = Branchname_gen.Data();
385 if(algorithm==1)fJetBranch = Branchname_rec.Data();
387 TClonesArray* jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
389 printf(" Tere are no Branch named %s \n",fJetBranch.Data());
392 Int_t nj = jets->GetEntriesFast();
393 if (fDebug) printf("There are %5d jets in the event \n", nj);
396 Jet_n[algorithm] = nj;
397 for(int njet =0;njet<nj;njet++){
398 jetsAOD = (AliAODJet*) (jets->At(njet));
399 Jet_pt [algorithm][njet] = jetsAOD->Pt();
400 Jet_phi [algorithm][njet] = jetsAOD->Phi();
401 Jet_eta [algorithm][njet] = jetsAOD->Eta();
402 double eta_cut_Jet=0.5;
403 TRefArray *reftracks = jetsAOD->GetRefTracks();
404 int ntracks=reftracks->GetEntriesFast();
405 if(TMath::Abs(Jet_eta[algorithm][njet])<eta_cut_Jet){
406 //------------calc max pt in Jet----------------
408 double sumtrackpt=0;//test
409 for(int ntr=0;ntr<ntracks;ntr++){// calc. max pt of track which is in Jet
410 AliAODTrack *AODtrack = dynamic_cast<AliAODTrack*>(reftracks->At(ntr));
414 if(AODtrack->TestFilterMask(Filtermask))bgoodT=true;
417 if(AODtrack->IsHybridGlobalConstrainedGlobal())bgoodT=true; //for hybrid Track cuts
420 trackpt = AODtrack->Pt();
421 sumtrackpt += trackpt;
422 if(trackpt>maxpt[algorithm][njet]){
423 maxpt[algorithm][njet] = trackpt;
429 fH1jetMCAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
430 if(maxpt[algorithm][njet]<40)fH1jetMCAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
431 if(ntracks>1) fH1jetMCAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
432 if(ntracks>2) fH1jetMCAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
433 if(ntracks>3) fH1jetMCAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
434 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetMCAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
435 fH2jetMCAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
438 fH1jetAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
439 if(maxpt[algorithm][njet]<40)fH1jetAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
440 if(ntracks>1) fH1jetAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
441 if(ntracks>2) fH1jetAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
442 if(ntracks>3) fH1jetAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
443 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
444 fH2jetAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
448 if(!(IsMC&&algorithm==1))continue;
449 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
450 double eta_cut_Jet=0.5;
451 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
452 //Find muched jet pare=====================================
453 for(int cut=0;cut<6;cut++){
455 for(int njetAOD=0;njetAOD<Jet_n[1];njetAOD++){
456 jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
458 jetsAOD = (AliAODJet*) (jets->At(njetAOD));
459 TRefArray *reftracks = jetsAOD->GetRefTracks();
460 int ntracks=reftracks->GetEntriesFast();
461 if(cut==1){if(maxpt[1][njetAOD]>=40.)continue;}
462 if(cut==2){if(ntracks==1)continue;}
463 if(cut==3){if(ntracks<=2)continue;}
464 if(cut==4){if(ntracks<=3)continue;}
465 if(cut==5){if(maxpt[1][njetAOD]>=40.)continue;if(ntracks==1)continue;}
466 double DelR = DeltaR(Jet_phi[0][njetMC],Jet_phi[1][njetAOD],Jet_eta[0][njetMC],Jet_eta[1][njetAOD]);
468 nearrecID[njetMC]=njetAOD;
475 for(int njet =0;njet<Jet_n[0];njet++){
476 double DelR = DeltaR(Jet_phi[0][njet],Jet_phi[1][nearrecID[njetMC]],Jet_eta[0][njet],Jet_eta[1][nearrecID[njetMC]]);
482 if((min_R<0.4)&&(neargenID==njetMC))fFIND[cut][njetMC]=true;
485 //======================================================
490 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
491 double eta_cut_Jet=0.5;
492 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
493 for(int cut=0;cut<6;cut++){
494 if(fFIND[cut][njetMC]==true){
495 fH1jetMCAKT04_match[cut]->Fill(Jet_pt[0][njetMC]);
496 fH2jetMCAKT04_Eratio[cut]->Fill(Jet_pt[0][njetMC],Jet_pt[1][nearrecID[njetMC]]/Jet_pt[0][njetMC]);
503 // End of jet analysis ---------
505 PostData(1, fHistList);
509 //________________________________________________________________________
510 void AliAnalysisTaskCheckSingleTrackJetRejection::Terminate(Option_t *)
512 // Terminate analysis
514 if (fDebug) printf("AnalysisTaskPt: Terminate() \n");
519 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::JetSelected(AliAODJet *jet){
520 Bool_t selected = false;
522 if(!jet)return selected;
524 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
532 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaPhi(Double_t phi1,Double_t phi2){
533 Float_t pi=TMath::Pi();
534 Double_t dphi = phi1-phi2;
535 if (dphi<(-1./2*pi))dphi = dphi +2*pi;
536 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
539 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaR(Double_t phi1,Double_t phi2,Double_t eta1,Double_t eta2){
540 Float_t pi=TMath::Pi();
541 Double_t dphi = DeltaPhi(phi1,phi2);
542 if (dphi<(-1./2*pi))dphi = dphi +2*pi;
543 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
544 Double_t Deta = eta1 - eta2;
545 Double_t deltaR = sqrt(pow(dphi,2)+pow(Deta,2));
549 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
551 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
552 // This is to called in Notify and should provide the path to the AOD/ESD file
553 // Copied from AliAnalysisTaskJetSpectrum2
556 TString file(currFile);
560 if(file.Contains("root_archive.zip#")){
561 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
562 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
563 file.Replace(pos+1,20,"");
566 // not an archive take the basename....
567 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
569 // Printf("%s",file.Data());
572 TFile *filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really
573 //test the existance of a file in a archive so we have to lvie with open error message from root
575 // next trial fetch the histgram file
576 filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
578 // not a severe condition but inciate that we have no information
582 // find the tlist we want to be independtent of the name so use the Tkey
583 TKey* key = (TKey*)filexsec->GetListOfKeys()->At(0);
588 TList *list = dynamic_cast<TList*>(key->ReadObj());
593 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
594 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
597 } // no tree pyxsec.root
599 TTree *xtree = (TTree*)filexsec->Get("Xsection");
605 Double_t xsection = 0;
606 xtree->SetBranchAddress("xsection",&xsection);
607 xtree->SetBranchAddress("ntrials",&ntrials);
617 //Float_t AliAnalysisTaskCheckSingleTrackJetRejection::GetTotalEvents(const char* currFile){
618 // Float_t totalevent;
619 // TString file_es(currFile);
620 // if(file_es.Contains("root_archive.zip#")){
621 // Ssiz_t pos1 = file_es.Index("root_archive",12,TString::kExact);
622 // Ssiz_t pos = file_es.Index("#",1,pos1,TString::kExact);
623 // file_es.Replace(pos+1,20,"");
626 // // not an archive take the basename....
627 // file_es.ReplaceAll(gSystem->BaseName(file_es.Data()),"");
630 // TString cAdd = "";
631 // cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
632 // cAdd += Form("B%d",(int)BackM);
633 // cAdd += Form("_Filter%05d",Filtermask);
634 // cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
635 // cAdd += Form("_Skip%02d",SkipCone);
636 // TString Dirname,Listname;
637 // Dirname = Form("PWG4_cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
638 // Listname = Form("pwg4cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
640 // TFile *feventstat = TFile::Open(Form("%s%s",file_es.Data(),"JetTasksOutput.root"));
641 // gROOT->Cd(Dirname.Data());
642 // TList *templist = (TList*)gROOT->FindObject(Listname.Data());
643 // TH1F* temphist = (TH1F*)templist->FindObject("fh1Trials");
644 // totalevent = temphist->Integral();
645 // //cout<<temphist->Integral()<<endl;
646 // delete feventstat;
647 // return totalevent;