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