PWGJE/UserTasks/AliAnalysisTaskCheckSingleTrackJetRejection.cxx

   1 //#include <string.h>
   2 //#include <TStyle.h>
   3 #include <list>
   4 #include <string>
   5
   6 #include "TTree.h"
   7 #include "TCanvas.h"
   8 #include "AliAnalysisTask.h"
   9 #include "AliInputEventHandler.h"
  10 #include "AliESDtrack.h"
  11 #include "AliAODVertex.h"
  12 #include "AliAODCluster.h"
  13
  14 #include <TROOT.h>
  15 #include <TRandom.h>
  16 #include <TSystem.h>
  17 #include <TInterpreter.h>
  18 #include <TChain.h>
  19 #include <TFile.h>
  20 #include <TKey.h>
  21 #include <TH1F.h>
  22 #include <TH2F.h>
  23 #include <TH3F.h>
  24 #include <TProfile.h>
  25 #include <TList.h>
  26 #include <TLorentzVector.h>
  27 #include <TClonesArray.h>
  28 #include <TRefArray.h>
  29
  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"
  50 #include "AliStack.h"
  51
  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"
  59
  60
  61
  62 #include "AliAnalysisTaskCheckSingleTrackJetRejection.h"
  63 #include "AliAnalysisTaskPhiCorrelations.h"
  64 #include "AliAnalysisHelperJetTasks.h"
  65 #include "AliPWG4HighPtQAMC.h"
  66
  67
  68
  69
  70 using std::cout;
  71 using std::endl;
  72
  73 ClassImp(AliAnalysisTaskCheckSingleTrackJetRejection)
  74
  75                                 //________________________________________________________________________
  76                                 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection():
  77                                                                 AliAnalysisTaskSE(),
  78                                                                 fUseAODInput(kFALSE),
  79                                                                 fFillAOD(kFALSE),
  80                                                                 fNonStdFile(""),
  81                                                                 fJetBranch("jets"),
  82                                                                 fAODIn(0x0),
  83                                                                 fAODOut(0x0),
  84                                                                 fAODExtension(0x0),
  85                                                                 JFAlg("ANTIKT"),
  86                                                                 Radius(0.4),
  87                                                                 Filtermask(256),
  88                                                                 BackM(0),
  89                                                                 TrackPtcut(0.15),
  90                                                                 SkipCone(0),
  91                                                                 IsMC(kTRUE),
  92                                                                 fHistList(0x0),
  93                                                                 fxsec(0.),
  94                                                                 ftrial(1.),
  95                                                                 fJetRecEtaWindow(0.5),
  96                                                                 fMinJetPt(0.15),
  97                                                                 fH1Xsec(0x0),
  98                                                                 fH1Trials(0x0),
  99                                                                 fH1Events(0x0),
 100
 101                                                                 fH2jetMCAKT04_Jetpt_maxpt(0x0),
 102                                                                 fH2jetAKT04_Jetpt_maxpt  (0x0)
 103
 104 {
 105                                 for(int j=0;j<6;j++){
 106                                                                 fH1jetMCAKT04_pt         [j]=0;
 107                                                                 fH1jetAKT04_pt           [j]=0;
 108
 109                                                                 fH2jetMCAKT04_Eratio       [j]=0;
 110                                                                 fH1jetMCAKT04_match        [j]=0;
 111                                 }
 112
 113                                 // Default constructor
 114 }
 115
 116 //________________________________________________________________________
 117 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection(const char *name):
 118                                 AliAnalysisTaskSE(name),
 119                                 fUseAODInput(kFALSE),
 120                                 fFillAOD(kFALSE),
 121                                 fNonStdFile(""),
 122                                 fJetBranch("jets"),
 123                                 fAODIn(0x0),
 124                                 fAODOut(0x0),
 125                                 fAODExtension(0x0),
 126                                 JFAlg("ANTIKT"),
 127                                 Radius(0.4),
 128                                 Filtermask(256),
 129                                 BackM(0),
 130                                 TrackPtcut(0.15),
 131                                 SkipCone(0),
 132                                 IsMC(kTRUE),
 133                                 fHistList(0x0),
 134                                 fxsec(0.),
 135                                 ftrial(1.),
 136                                 fJetRecEtaWindow(0.5),
 137                                 fMinJetPt(0.15),
 138                                 fH1Xsec(0x0),
 139                                 fH1Trials(0x0),
 140                                 fH1Events(0x0),
 141                                 fH2jetMCAKT04_Jetpt_maxpt(0x0),
 142                                 fH2jetAKT04_Jetpt_maxpt  (0x0)
 143
 144 {
 145
 146                                 for(int j=0;j<6;j++){
 147                                                                 fH1jetMCAKT04_pt         [j]=0;
 148                                                                 fH1jetAKT04_pt           [j]=0;
 149
 150                                                                 fH2jetMCAKT04_Eratio       [j]=0;
 151                                                                 fH1jetMCAKT04_match        [j]=0;
 152                                 }
 153
 154                                 // Constructor
 155                                 // Define input and output slots here
 156                                 // Input slot #0 works with a TChain
 157                                 DefineInput(0, TChain::Class());
 158                                 // Output slot #0 id reserved by the base class for AOD
 159                                 // Output slot #1 writes into a TH1 container
 160                                 DefineOutput(1, TList::Class());
 161
 162 }
 163
 164 //________________________________________________________________________
 165 void AliAnalysisTaskCheckSingleTrackJetRejection::UserCreateOutputObjects()
 166 {
 167                                 // Create histograms
 168                                 // Called once
 169
 170                                 if (!fHistList){ fHistList = new TList();fHistList->SetOwner(kTRUE); cout<<"TList is created for output "<<endl;}
 171
 172                                 Bool_t oldStatus = TH1::AddDirectoryStatus();
 173                                 TH1::AddDirectory(kFALSE);
 174
 175                                 char *histname;
 176                                 if(IsMC){
 177                                                                 fH1Xsec           = new TProfile("Xsec","Xsec",1,0,1);
 178                                                                 fH1Trials         = new TH1F    ("Trials","Trials",1,0,1);
 179                                                                 histname = Form("fH2jetMCAKT04_Jetpt_maxpt");
 180                                                                 fH2jetMCAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
 181                                                                 histname = Form("fH2jetAKT04_Jetpt_maxpt");
 182                                                                 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
 183                                                                 for(int j=0;j<6;j++){
 184                                                                                                 histname = Form("fH1jetMCAKT04_pt%d",j);
 185                                                                                                 fH1jetMCAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
 186                                                                                                 histname = Form("fH2jetAKT04_pt%d",j);
 187                                                                                                 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
 188
 189                                                                                                 histname = Form("fH2jetMCAKT04_Eratio%d",j);
 190                                                                                                 fH2jetMCAKT04_Eratio[j] = new TH2F(histname,histname,200,0,200,100,0,30);
 191                                                                                                 histname = Form("fH1jetMCAKT04_match%d",j);
 192                                                                                                 fH1jetMCAKT04_match[j] = new TH1F(histname,histname,200,0,200);
 193                                                                 }
 194                                                                 fHistList->Add(fH1Xsec);
 195                                                                 fHistList->Add(fH1Trials);
 196                                                                 fHistList->Add(fH2jetMCAKT04_Jetpt_maxpt);
 197                                                                 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
 198                                                                 for(int j=0;j<6;j++){
 199                                                                                                 fHistList->Add(fH1jetAKT04_pt[j]);
 200                                                                                                 fHistList->Add(fH1jetMCAKT04_pt[j]);
 201                                                                                                 fHistList->Add(fH2jetMCAKT04_Eratio[j]);
 202                                                                                                 fHistList->Add(fH1jetMCAKT04_match[j]);
 203                                                                 }
 204                                 }
 205                                 else    {
 206                                                                 fH1Events         = new TH1F("Events","Number of Events",1,0,1);
 207                                                                 histname = Form("fH2jetAKT04_Jetpt_maxpt");
 208                                                                 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
 209                                                                 for(int j=0;j<6;j++){
 210                                                                                                 histname = Form("fH2jetAKT04_pt%d",j);
 211                                                                                                 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
 212                                                                 }
 213                                                                 fHistList->Add(fH1Events);
 214                                                                 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
 215                                                                 for(int j=0;j<6;j++){
 216                                                                                                 fHistList->Add(fH1jetAKT04_pt[j]);
 217                                                                 }
 218                                 }
 219
 220
 221                                 // =========== Switch on Sumw2 for all histos ===========
 222                                 for (Int_t i=0; i<fHistList->GetEntries(); ++i)
 223                                 {
 224                                                                 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 225                                                                 if (h1)
 226                                                                 {
 227                                                                                                 h1->Sumw2();
 228                                                                                                 continue;
 229                                                                 }
 230                                                                 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 231                                                                 if(hn)hn->Sumw2();
 232                                 }
 233                                 TH1::AddDirectory(oldStatus);
 234
 235
 236                                 PostData(1,fHistList);
 237
 238 }
 239
 240 //----------------------------------------------------------------------
 241 void AliAnalysisTaskCheckSingleTrackJetRejection::Init()
 242 {
 243                                 // Initialization
 244                                 if (fDebug) printf("AnalysisTaskCheckSingleTrackJetRejection::Init() \n");
 245
 246 }
 247
 248 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::Notify()
 249 {
 250                                 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 251                                 fAODOut = AODEvent();
 252                                 if(fNonStdFile.Length()!=0){
 253                                                                 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 254                                                                 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 255                                                                 if(fAODExtension){
 256                                                                                                 if(fDebug>1)Printf("AODExtension found for %s ",fNonStdFile.Data());
 257                                                                 }
 258                                 }
 259
 260                                 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 261                                 fxsec=0.;
 262                                 ftrial=1.;
 263
 264                                 if(tree){
 265                                                                 TFile *curfile = tree->GetCurrentFile();
 266                                                                 if(!curfile){
 267                                                                                                 Error("Notify","No current file");
 268                                                                                                 return kFALSE;
 269                                                                 }
 270
 271                                                                 if(IsMC){
 272                                                                                                 PythiaInfoFromFile(curfile->GetName(),fxsec,ftrial);
 273                                                                                                 //cout<<" Xsec "<<fxsec<<" trial "<<ftrial<<endl;
 274                                                                                                 fH1Xsec  ->Fill(0.,fxsec);
 275                                                                                                 fH1Trials->Fill(0.,ftrial);
 276                                                                 }
 277                                                                 else{
 278                                                                                                 //Float_t totalEvent;
 279                                                                                                 //totalEvent = GetTotalEvents(curfile->GetName());
 280                                                                                                 //fH1Events->Fill(0.,totalEvent);
 281                                                                 }
 282
 283                                 }
 284
 285                                 printf("Reading File %s ",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 286                                 return kTRUE;
 287 }
 288
 289 void AliAnalysisTaskCheckSingleTrackJetRejection::FinishTaskOutput()
 290 {
 291 }
 292
 293
 294
 295 //________________________________________________________________________
 296 void AliAnalysisTaskCheckSingleTrackJetRejection::UserExec(Option_t *)
 297 {
 298
 299
 300                                 // Main loop (called each event)
 301                                 // Execute analysis for current event
 302
 303                                 AliAODEvent *fAOD;
 304                                 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 305                                 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
 306                                                                 fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
 307                                                                 if(fUseAODInput) fAODIn = fAOD;
 308                                                                 if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
 309                                 }
 310                                 else {
 311                                                                 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 312                                                                 if( handler && handler->InheritsFrom("AliAODHandler") ) {
 313                                                                                                 fAOD = ((AliAODHandler*)handler)->GetAOD();
 314                                                                                                 fAODIn = fAOD;
 315                                                                                                 if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
 316                                                                 }
 317                                 }
 318
 319                                 if(!fAODIn && !fUseAODInput){ // case we have AOD in input & output and want jets from output
 320                                                                 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 321                                                                 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
 322                                                                                                 fAODIn = ((AliAODHandler*)outHandler)->GetAOD();
 323                                                                                                 if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
 324                                                                 }
 325                                 }
 326
 327                                 //if(fNonStdFile.Length()!=0){
 328                                 //                              // case we have an AOD extension - fetch the jets from the extended output
 329                                 //                              AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 330                                 //                              fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 331                                 //                              if(!fAODExtension){
 332                                 //                                                              if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
 333                                 //                              }
 334                                 //}
 335                                 //fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 336                                 if (!fAODIn) {
 337                                                                 Printf("ERROR %s : fAODIn not available",(char*)__FILE__);
 338                                                                 return;
 339                                 }
 340
 341                                 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 342                                                                 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 343                                 if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
 344                                                                 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
 345                                                                 return;
 346                                 }
 347                                 if(!IsMC)fH1Events->Fill(0);
 348
 349
 350
 351                                 // start jet analysis
 352
 353                                 Double_t Jet_n  [20];
 354                                 Double_t Jet_pt [20][10000];
 355                                 Double_t Jet_eta[20][10000];
 356                                 Double_t Jet_phi[20][10000];
 357
 358                                 for(int i=0;i<20;i++){
 359                                                                 Jet_n[i]=0;
 360                                                                 for(int j=0;j<10000;j++){
 361                                                                                                 Jet_pt[i][j]=0.;
 362                                                                                                 Jet_phi[i][j]=999.;
 363                                                                                                 Jet_eta[i][j]=999.;
 364                                                                 }
 365                                 }
 366
 367                                 // get AOD event from input/ouput
 368                                 TString cAdd = "";
 369                                 cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
 370                                 cAdd += Form("B%d",(int)BackM);
 371                                 cAdd += Form("_Filter%05d",Filtermask);
 372                                 cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
 373                                 cAdd += Form("_Skip%02d",SkipCone);
 374                                 TString Branchname_gen,Branchname_rec;
 375                                 //Branchname_gen = Form("clustersMCKINE2_%s%s",JFAlg.Data(),cAdd.Data());
 376                                 Branchname_gen = Form("clustersAODMC2_%s%s",JFAlg.Data(),cAdd.Data());
 377                                 Branchname_rec = Form("clustersAOD_%s%s",JFAlg.Data(),cAdd.Data());
 378
 379
 380                                 bool fFIND[6][1000];
 381                                 double maxpt[2][1000];for(int i=0;i<1000;i++){maxpt[0][i]=0;maxpt[1][i]=0;}
 382                                 int nearrecID[1000];  for(int i=0;i<1000;i++){nearrecID[i]=99999;}
 383                                 AliAODJet* jetsAOD;
 384                                 for(int algorithm=0;algorithm<2;algorithm++){
 385                                                                 //for LHC11a1  LHC11a2 official
 386                                                                 if((!IsMC&&algorithm==0))continue;
 387
 388                                                                 if(algorithm==0)fJetBranch   = Branchname_gen.Data();
 389                                                                 if(algorithm==1)fJetBranch   = Branchname_rec.Data();
 390
 391                                                                 TClonesArray* jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
 392                                                                 if(!jets){
 393                                                                                                 printf(" Tere are no Branch named %s \n",fJetBranch.Data());
 394                                                                                                 continue;
 395                                                                 }
 396                                                                 Int_t nj = jets->GetEntriesFast();
 397                                                                 if (fDebug) printf("There are %5d jets in the event \n", nj);
 398
 399
 400                                                                 Jet_n[algorithm] = nj;
 401                                                                 for(int njet =0;njet<nj;njet++){
 402                                                                                                 jetsAOD = (AliAODJet*) (jets->At(njet));
 403                                                                                                 Jet_pt   [algorithm][njet] = jetsAOD->Pt();
 404                                                                                                 Jet_phi  [algorithm][njet] = jetsAOD->Phi();
 405                                                                                                 Jet_eta  [algorithm][njet] = jetsAOD->Eta();
 406                                                                                                 double eta_cut_Jet=0.5;
 407                                                                                                 TRefArray *reftracks = jetsAOD->GetRefTracks();
 408                                                                                                 int ntracks=reftracks->GetEntriesFast();
 409                                                                                                 if(TMath::Abs(Jet_eta[algorithm][njet])<eta_cut_Jet){
 410                                                                                                                                 //------------calc max pt in Jet----------------
 411                                                                                                                                 double trackpt;
 412                                                                                                                                 double sumtrackpt=0;//test
 413                                                                                                                                 for(int ntr=0;ntr<ntracks;ntr++){// calc. max pt of track which is in Jet
 414                                                                                                                                                                 AliAODTrack *AODtrack = dynamic_cast<AliAODTrack*>(reftracks->At(ntr));
 415                                                                                                                                                                 if(AODtrack){
 416                                                                                                                                                                                                 Bool_t bgoodT=false;
 417                                                                                                                                                                                                 if(Filtermask!=768){
 418                                                                                                                                                                                                                                 if(AODtrack->TestFilterMask(Filtermask))bgoodT=true;
 419                                                                                                                                                                                                 }
 420                                                                                                                                                                                                 else{
 421                                                                                                                                                                                                                                 if(AODtrack->IsHybridGlobalConstrainedGlobal())bgoodT=true; //for hybrid Track cuts
 422                                                                                                                                                                                                 }
 423                                                                                                                                                                                                 if(!bgoodT)continue;
 424                                                                                                                                                                                                 trackpt = AODtrack->Pt();
 425                                                                                                                                                                                                 sumtrackpt += trackpt;
 426                                                                                                                                                                                                 if(trackpt>maxpt[algorithm][njet]){
 427                                                                                                                                                                                                                                 maxpt[algorithm][njet] = trackpt;
 428                                                                                                                                                                                                 }
 429                                                                                                                                                                 }
 430                                                                                                                                 }// track Loop
 431
 432                                                                                                                                 if(algorithm==0){
 433                                                                                                                                                                 fH1jetMCAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
 434                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetMCAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
 435                                                                                                                                                                 if(ntracks>1)     fH1jetMCAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
 436                                                                                                                                                                 if(ntracks>2)     fH1jetMCAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
 437                                                                                                                                                                 if(ntracks>3)     fH1jetMCAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
 438                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetMCAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
 439                                                                                                                                                                 fH2jetMCAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
 440                                                                                                                                 }
 441                                                                                                                                 if(algorithm==1){
 442                                                                                                                                                                 fH1jetAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
 443                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
 444                                                                                                                                                                 if(ntracks>1)     fH1jetAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
 445                                                                                                                                                                 if(ntracks>2)     fH1jetAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
 446                                                                                                                                                                 if(ntracks>3)     fH1jetAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
 447                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
 448                                                                                                                                                                 fH2jetAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
 449                                                                                                                                 }
 450                                                                                                 }
 451                                                                 }
 452                                                                 if(!(IsMC&&algorithm==1))continue;
 453                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
 454                                                                                                 double eta_cut_Jet=0.5;
 455                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
 456                                                                                                                                 //Find muched jet pare=====================================
 457                                                                                                                                 for(int cut=0;cut<6;cut++){
 458                                                                                                                                                                 double min_R=10.;
 459                                                                                                                                                                 for(int njetAOD=0;njetAOD<Jet_n[1];njetAOD++){
 460                                                                                                                                                                   jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
 461                                                                                                                                                                   if(!jets)continue;
 462                                                                                                                                                                   jetsAOD = (AliAODJet*) (jets->At(njetAOD));
 463                                                                                                                                                                   TRefArray *reftracks = jetsAOD->GetRefTracks();
 464                                                                                                                                                                   int ntracks=reftracks->GetEntriesFast();
 465                                                                                                                                                                                                 if(cut==1){if(maxpt[1][njetAOD]>=40.)continue;}
 466                                                                                                                                                                                                 if(cut==2){if(ntracks==1)continue;}
 467                                                                                                                                                                                                 if(cut==3){if(ntracks<=2)continue;}
 468                                                                                                                                                                                                 if(cut==4){if(ntracks<=3)continue;}
 469                                                                                                                                                                                                 if(cut==5){if(maxpt[1][njetAOD]>=40.)continue;if(ntracks==1)continue;}
 470                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njetMC],Jet_phi[1][njetAOD],Jet_eta[0][njetMC],Jet_eta[1][njetAOD]);
 471                                                                                                                                                                                                 if(DelR<min_R){
 472                                                                                                                                                                                                                                 nearrecID[njetMC]=njetAOD;
 473                                                                                                                                                                                                                                 min_R=DelR;
 474                                                                                                                                                                                                 }
 475                                                                                                                                                                 }
 476                                                                                                                                                                 if(min_R<0.4){
 477                                                                                                                                                                                                 min_R=10.;
 478                                                                                                                                                                                                 int neargenID=99999;
 479                                                                                                                                                                                                 for(int njet =0;njet<Jet_n[0];njet++){
 480                                                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njet],Jet_phi[1][nearrecID[njetMC]],Jet_eta[0][njet],Jet_eta[1][nearrecID[njetMC]]);
 481                                                                                                                                                                                                                                 if(DelR<min_R){
 482                                                                                                                                                                                                                                                                 neargenID=njet;
 483                                                                                                                                                                                                                                                                 min_R=DelR;
 484                                                                                                                                                                                                                                 }
 485                                                                                                                                                                                                 }
 486                                                                                                                                                                                                 if((min_R<0.4)&&(neargenID==njetMC))fFIND[cut][njetMC]=true;
 487                                                                                                                                                                 }
 488                                                                                                                                 }
 489                                                                                                                                 //======================================================
 490                                                                                                 }
 491                                                                 }
 492                                 }//algorithm
 493                                 if(IsMC){
 494                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
 495                                                                                                 double eta_cut_Jet=0.5;
 496                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
 497                                                                                                                                 for(int cut=0;cut<6;cut++){
 498                                                                                                                                                                 if(fFIND[cut][njetMC]==true){
 499                                                                                                                                                                                                 fH1jetMCAKT04_match[cut]->Fill(Jet_pt[0][njetMC]);
 500                                                                                                                                                                                                 fH2jetMCAKT04_Eratio[cut]->Fill(Jet_pt[0][njetMC],Jet_pt[1][nearrecID[njetMC]]/Jet_pt[0][njetMC]);
 501                                                                                                                                                                 }
 502                                                                                                                                 }
 503                                                                                                 }//eta window
 504                                                                 }//njet loop
 505                                 }
 506
 507                                 // End of jet analysis ---------
 508                                 // Post output data.
 509                                 PostData(1, fHistList);
 510                                 return;
 511 }
 512
 513 //________________________________________________________________________
 514 void AliAnalysisTaskCheckSingleTrackJetRejection::Terminate(Option_t *)
 515 {
 516                                 // Terminate analysis
 517                                 //
 518                                 if (fDebug) printf("AnalysisTaskPt: Terminate() \n");
 519
 520 }
 521
 522
 523 Bool_t  AliAnalysisTaskCheckSingleTrackJetRejection::JetSelected(AliAODJet *jet){
 524                                 Bool_t selected = false;
 525
 526                                 if(!jet)return selected;
 527
 528                                 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
 529                                                                 selected = kTRUE;
 530                                 }
 531                                 return selected;
 532
 533 }
 534
 535
 536 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaPhi(Double_t phi1,Double_t phi2){
 537                                 Float_t pi=TMath::Pi();
 538                                 Double_t dphi = phi1-phi2;
 539                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
 540                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
 541                                 return dphi;
 542 }
 543 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaR(Double_t phi1,Double_t phi2,Double_t eta1,Double_t eta2){
 544                                 Float_t pi=TMath::Pi();
 545                                 Double_t dphi = DeltaPhi(phi1,phi2);
 546                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
 547                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
 548                                 Double_t Deta = eta1 - eta2;
 549                                 Double_t deltaR = sqrt(pow(dphi,2)+pow(Deta,2));
 550                                 return deltaR;
 551 }
 552
 553 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 554                                 //
 555                                 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 556                                 // This is to called in Notify and should provide the path to the AOD/ESD file
 557                                 // Copied from AliAnalysisTaskJetSpectrum2
 558                                 //
 559
 560                                 TString file(currFile);
 561                                 fXsec = 0;
 562                                 fTrials = 1;
 563
 564                                 if(file.Contains("root_archive.zip#")){
 565                                                                 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 566                                                                 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 567                                                                 file.Replace(pos+1,20,"");
 568                                 }
 569                                 else {
 570                                                                 // not an archive take the basename....
 571                                                                 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 572                                 }
 573                                 //  Printf("%s",file.Data());
 574
 575
 576                                 TFile *filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really
 577                                 //test the existance of a file in a archive so we have to lvie with open error message from root
 578                                                                 if(!filexsec){
 579                                                                                                 // next trial fetch the histgram file
 580                                                                                                 filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 581                                                                                                 if(!filexsec){
 582                                                                                                                                 // not a severe condition but inciate that we have no information
 583                                                                                                                                 return kFALSE;
 584                                                                                                 }
 585                                                                                                 else{
 586                                                                                                                                 // find the tlist we want to be independtent of the name so use the Tkey
 587                                                                                                                                 TKey* key = (TKey*)filexsec->GetListOfKeys()->At(0);
 588                                                                                                                                 if(!key){
 589                                                                                                                                                                 filexsec->Close();
 590                                                                                                                                                                 return kFALSE;
 591                                                                                                                                 }
 592                                                                                                                                 TList *list = dynamic_cast<TList*>(key->ReadObj());
 593                                                                                                                                 if(!list){
 594                                                                                                                                                                 filexsec->Close();
 595                                                                                                                                                                 return kFALSE;
 596                                                                                                                                 }
 597                                                                                                                                 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 598                                                                                                                                 fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 599                                                                                                                                 filexsec->Close();
 600                                                                                                 }
 601                                                                 } // no tree pyxsec.root
 602                                                                 else {
 603                                                                                                 TTree *xtree = (TTree*)filexsec->Get("Xsection");
 604                                                                                                 if(!xtree){
 605                                                                                                                                 filexsec->Close();
 606                                                                                                                                 return kFALSE;
 607                                                                                                 }
 608                                                                                                 UInt_t   ntrials  = 0;
 609                                                                                                 Double_t  xsection  = 0;
 610                                                                                                 xtree->SetBranchAddress("xsection",&xsection);
 611                                                                                                 xtree->SetBranchAddress("ntrials",&ntrials);
 612                                                                                                 xtree->GetEntry(0);
 613                                                                                                 fTrials = ntrials;
 614                                                                                                 fXsec = xsection;
 615                                                                                                 filexsec->Close();
 616                                                                 }
 617                                 return kTRUE;
 618 }
 619
 620
 621 //Float_t AliAnalysisTaskCheckSingleTrackJetRejection::GetTotalEvents(const char* currFile){
 622 //                              Float_t totalevent;
 623 //                              TString file_es(currFile);
 624 //                              if(file_es.Contains("root_archive.zip#")){
 625 //                                                              Ssiz_t pos1 = file_es.Index("root_archive",12,TString::kExact);
 626 //                                                              Ssiz_t pos = file_es.Index("#",1,pos1,TString::kExact);
 627 //                                                              file_es.Replace(pos+1,20,"");
 628 //                              }
 629 //                              else {
 630 //                                                              // not an archive take the basename....
 631 //                                                              file_es.ReplaceAll(gSystem->BaseName(file_es.Data()),"");
 632 //                              }
 633 //
 634 //                              TString cAdd = "";
 635 //                              cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
 636 //                              cAdd += Form("B%d",(int)BackM);
 637 //                              cAdd += Form("_Filter%05d",Filtermask);
 638 //                              cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
 639 //                              cAdd += Form("_Skip%02d",SkipCone);
 640 //                              TString Dirname,Listname;
 641 //                              Dirname  = Form("PWG4_cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
 642 //                              Listname = Form("pwg4cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
 643 //
 644 //                              TFile *feventstat = TFile::Open(Form("%s%s",file_es.Data(),"JetTasksOutput.root"));
 645 //                              gROOT->Cd(Dirname.Data());
 646 //                              TList *templist     = (TList*)gROOT->FindObject(Listname.Data());
 647 //                              TH1F* temphist = (TH1F*)templist->FindObject("fh1Trials");
 648 //                              totalevent = temphist->Integral();
 649 //                              //cout<<temphist->Integral()<<endl;
 650 //                              delete feventstat;
 651 //                              return totalevent;
 652 //
 653 //}
 654
 655