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