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 "AliESDEvent.h"
  37 #include "AliAODEvent.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODInputHandler.h"
  40 #include "AliAODTrack.h"
  41 #include "AliVParticle.h"
  42 #include "AliAODJet.h"
  43 #include "AliAODJetEventBackground.h"
  44 #include "AliMCParticle.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliMCEventHandler.h"
  47 #include "AliMCEvent.h"
  48 #include "AliStack.h"
  49
  50 #include "AliAODHeader.h"
  51 #include "AliAODMCHeader.h"
  52 #include "AliJetKineReaderHeader.h"
  53 #include "AliGenCocktailEventHeader.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliGenEventHeader.h"
  56 #include "AliGenDPMjetEventHeader.h"
  57
  58
  59
  60 #include "AliAnalysisTaskCheckSingleTrackJetRejection.h"
  61 #include "AliAnalysisTaskPhiCorrelations.h"
  62 #include "AliAnalysisHelperJetTasks.h"
  63 #include "AliPWG4HighPtQAMC.h"
  64
  65
  66 using std::cout;
  67 using std::endl;
  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_ptmatch        [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_ptmatch        [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("fH1jetAKT04_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_ptmatch%d",j);
 189                                                                                                 fH1jetMCAKT04_ptmatch[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_ptmatch[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("fH1jetAKT04_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                                 //                              AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 327                                 //                              fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 328                                 //                              if(!fAODExtension){
 329                                 //                                                              if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
 330                                 //                              }
 331                                 //}
 332                                 //fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 333                                 if (!fAODIn) {
 334                                                                 Printf("ERROR %s : fAODIn not available",(char*)__FILE__);
 335                                                                 return;
 336                                 }
 337
 338                                 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 339                                                                 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 340                                 if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
 341                                                                 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
 342                                                                 return;
 343                                 }
 344                                 if(!IsMC)fH1Events->Fill(0);
 345
 346
 347
 348                                 // start jet analysis
 349
 350                                 Double_t Jet_n  [20];
 351                                 Double_t Jet_pt [20][10000];
 352                                 Double_t Jet_eta[20][10000];
 353                                 Double_t Jet_phi[20][10000];
 354
 355                                 for(int i=0;i<20;i++){
 356                                                                 Jet_n[i]=0;
 357                                                                 for(int j=0;j<10000;j++){
 358                                                                                                 Jet_pt[i][j]=0.;
 359                                                                                                 Jet_phi[i][j]=999.;
 360                                                                                                 Jet_eta[i][j]=999.;
 361                                                                 }
 362                                 }
 363
 364                                 // get AOD event from input/ouput
 365                                 TString cAdd = "";
 366                                 cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
 367                                 cAdd += Form("B%d",(int)BackM);
 368                                 cAdd += Form("_Filter%05d",Filtermask);
 369                                 cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
 370                                 cAdd += Form("_Skip%02d",SkipCone);
 371                                 TString Branchname_gen,Branchname_rec;
 372                                 //Branchname_gen = Form("clustersMCKINE2_%s%s",JFAlg.Data(),cAdd.Data());
 373                                 Branchname_gen = Form("clustersAODMC2_%s%s",JFAlg.Data(),cAdd.Data());
 374                                 Branchname_rec = Form("clustersAOD_%s%s",JFAlg.Data(),cAdd.Data());
 375
 376
 377                                 bool fFIND[6][1000];
 378                                 double maxpt[2][1000];for(int i=0;i<1000;i++){maxpt[0][i]=0;maxpt[1][i]=0;}
 379                                 int nearrecID[1000];  for(int i=0;i<1000;i++){nearrecID[i]=99999;}
 380                                 AliAODJet* jetsAOD;
 381                                 for(int algorithm=0;algorithm<2;algorithm++){
 382                                                                 //for LHC11a1  LHC11a2 official
 383                                                                 if((!IsMC&&algorithm==0))continue;
 384
 385                                                                 if(algorithm==0)fJetBranch   = Branchname_gen.Data();
 386                                                                 if(algorithm==1)fJetBranch   = Branchname_rec.Data();
 387
 388                                                                 TClonesArray* jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
 389                                                                 if(!jets){
 390                                                                                                 printf(" Tere are no Branch named %s \n",fJetBranch.Data());
 391                                                                                                 continue;
 392                                                                 }
 393                                                                 Int_t nj = jets->GetEntriesFast();
 394                                                                 if (fDebug) printf("There are %5d jets in the event \n", nj);
 395
 396
 397                                                                 Jet_n[algorithm] = nj;
 398                                                                 for(int njet =0;njet<nj;njet++){
 399                                                                                                 jetsAOD = (AliAODJet*) (jets->At(njet));
 400                                                                                                 Jet_pt   [algorithm][njet] = jetsAOD->Pt();
 401                                                                                                 Jet_phi  [algorithm][njet] = jetsAOD->Phi();
 402                                                                                                 Jet_eta  [algorithm][njet] = jetsAOD->Eta();
 403                                                                                                 double eta_cut_Jet=0.5;
 404                                                                                                 TRefArray *reftracks = jetsAOD->GetRefTracks();
 405                                                                                                 int ntracks=reftracks->GetEntriesFast();
 406                                                                                                 if(TMath::Abs(Jet_eta[algorithm][njet])<eta_cut_Jet){
 407                                                                                                                                 //------------calc max pt in Jet----------------
 408                                                                                                                                 double trackpt;
 409                                                                                                                                 double sumtrackpt=0;//test
 410                                                                                                                                 for(int ntr=0;ntr<ntracks;ntr++){// calc. max pt of track which is in Jet
 411                                                                                                                                                                 AliAODTrack *AODtrack = dynamic_cast<AliAODTrack*>(reftracks->At(ntr));
 412                                                                                                                                                                 if(AODtrack){
 413                                                                                                                                                                                                 Bool_t bgoodT=false;
 414                                                                                                                                                                                                 if(Filtermask!=272){
 415                                                                                                                                                                                                                                 if(AODtrack->TestFilterMask(Filtermask))bgoodT=true;
 416                                                                                                                                                                                                 }
 417                                                                                                                                                                                                 else{
 418                                                                                                                                                                                                                                 if(AODtrack->IsHybridGlobalConstrainedGlobal())bgoodT=true; //for hybrid Track cuts
 419                                                                                                                                                                                                 }
 420                                                                                                                                                                                                 if(!bgoodT)continue;
 421                                                                                                                                                                                                 trackpt = AODtrack->Pt();
 422                                                                                                                                                                                                 sumtrackpt += trackpt;
 423                                                                                                                                                                                                 if(trackpt>maxpt[algorithm][njet]){
 424                                                                                                                                                                                                                                 maxpt[algorithm][njet] = trackpt;
 425                                                                                                                                                                                                 }
 426                                                                                                                                                                 }
 427                                                                                                                                 }// track Loop
 428
 429                                                                                                                                 if(algorithm==0){
 430                                                                                                                                                                 fH1jetMCAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
 431                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetMCAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
 432                                                                                                                                                                 if(ntracks>1)     fH1jetMCAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
 433                                                                                                                                                                 if(ntracks>2)     fH1jetMCAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
 434                                                                                                                                                                 if(ntracks>3)     fH1jetMCAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
 435                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetMCAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
 436                                                                                                                                                                 fH2jetMCAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
 437                                                                                                                                 }
 438                                                                                                                                 if(algorithm==1){
 439                                                                                                                                                                 fH1jetAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
 440                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
 441                                                                                                                                                                 if(ntracks>1)     fH1jetAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
 442                                                                                                                                                                 if(ntracks>2)     fH1jetAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
 443                                                                                                                                                                 if(ntracks>3)     fH1jetAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
 444                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
 445                                                                                                                                                                 fH2jetAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
 446                                                                                                                                 }
 447                                                                                                 }
 448                                                                 }
 449                                                                 if(!(IsMC&&algorithm==1))continue;
 450                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
 451                                                                                                 double eta_cut_Jet=0.5;
 452                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
 453                                                                                                                                 //Find muched jet pare=====================================
 454                                                                                                                                 for(int cut=0;cut<6;cut++){
 455                                                                                                                                                                 double min_R=10.;
 456                                                                                                                                                                 for(int njetAOD=0;njetAOD<Jet_n[1];njetAOD++){
 457                                                                                                                                                                   jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
 458                                                                                                                                                                   if(!jets)continue;
 459                                                                                                                                                                   jetsAOD = (AliAODJet*) (jets->At(njetAOD));
 460                                                                                                                                                                   TRefArray *reftracks = jetsAOD->GetRefTracks();
 461                                                                                                                                                                   int ntracks=reftracks->GetEntriesFast();
 462                                                                                                                                                                                                 if(cut==1){if(maxpt[1][njetAOD]>=40.)continue;}
 463                                                                                                                                                                                                 if(cut==2){if(ntracks==1)continue;}
 464                                                                                                                                                                                                 if(cut==3){if(ntracks<=2)continue;}
 465                                                                                                                                                                                                 if(cut==4){if(ntracks<=3)continue;}
 466                                                                                                                                                                                                 if(cut==5){if(maxpt[1][njetAOD]>=40.)continue;if(ntracks==1)continue;}
 467                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njetMC],Jet_phi[1][njetAOD],Jet_eta[0][njetMC],Jet_eta[1][njetAOD]);
 468                                                                                                                                                                                                 if(DelR<min_R){
 469                                                                                                                                                                                                                                 nearrecID[njetMC]=njetAOD;
 470                                                                                                                                                                                                                                 min_R=DelR;
 471                                                                                                                                                                                                 }
 472                                                                                                                                                                 }
 473                                                                                                                                                                 if(min_R<0.4){
 474                                                                                                                                                                                                 min_R=10.;
 475                                                                                                                                                                                                 int neargenID=99999;
 476                                                                                                                                                                                                 for(int njet =0;njet<Jet_n[0];njet++){
 477                                                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njet],Jet_phi[1][nearrecID[njetMC]],Jet_eta[0][njet],Jet_eta[1][nearrecID[njetMC]]);
 478                                                                                                                                                                                                                                 if(DelR<min_R){
 479                                                                                                                                                                                                                                                                 neargenID=njet;
 480                                                                                                                                                                                                                                                                 min_R=DelR;
 481                                                                                                                                                                                                                                 }
 482                                                                                                                                                                                                 }
 483                                                                                                                                                                                                 if((min_R<0.4)&&(neargenID==njetMC))fFIND[cut][njetMC]=true;
 484                                                                                                                                                                 }
 485                                                                                                                                 }
 486                                                                                                                                 //======================================================
 487                                                                                                 }
 488                                                                 }
 489                                 }//algorithm
 490                                 if(IsMC){
 491                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
 492                                                                                                 double eta_cut_Jet=0.5;
 493                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
 494                                                                                                                                 for(int cut=0;cut<6;cut++){
 495                                                                                                                                                                 if(fFIND[cut][njetMC]==true){
 496                                                                                                                                                                                                 fH1jetMCAKT04_ptmatch[cut]->Fill(Jet_pt[0][njetMC]);
 497                                                                                                                                                                                                 fH2jetMCAKT04_Eratio[cut]->Fill(Jet_pt[0][njetMC],Jet_pt[1][nearrecID[njetMC]]/Jet_pt[0][njetMC]);
 498                                                                                                                                                                 }
 499                                                                                                                                 }
 500                                                                                                 }//eta window
 501                                                                 }//njet loop
 502                                 }
 503
 504                                 // End of jet analysis ---------
 505                                 // Post output data.
 506                                 PostData(1, fHistList);
 507                                 return;
 508 }
 509
 510 //________________________________________________________________________
 511 void AliAnalysisTaskCheckSingleTrackJetRejection::Terminate(Option_t *)
 512 {
 513                                 // Terminate analysis
 514                                 //
 515                                 if (fDebug) printf("AnalysisTaskPt: Terminate() \n");
 516
 517 }
 518
 519
 520 Bool_t  AliAnalysisTaskCheckSingleTrackJetRejection::JetSelected(AliAODJet *jet){
 521                                 Bool_t selected = false;
 522
 523                                 if(!jet)return selected;
 524
 525                                 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
 526                                                                 selected = kTRUE;
 527                                 }
 528                                 return selected;
 529
 530 }
 531
 532
 533 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaPhi(Double_t phi1,Double_t phi2){
 534                                 Float_t pi=TMath::Pi();
 535                                 Double_t dphi = phi1-phi2;
 536                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
 537                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
 538                                 return dphi;
 539 }
 540 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaR(Double_t phi1,Double_t phi2,Double_t eta1,Double_t eta2){
 541                                 Float_t pi=TMath::Pi();
 542                                 Double_t dphi = DeltaPhi(phi1,phi2);
 543                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
 544                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
 545                                 Double_t Deta = eta1 - eta2;
 546                                 Double_t deltaR = sqrt(pow(dphi,2)+pow(Deta,2));
 547                                 return deltaR;
 548 }
 549
 550 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 551                                 //
 552                                 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 553                                 // This is to called in Notify and should provide the path to the AOD/ESD file
 554                                 // Copied from AliAnalysisTaskJetSpectrum2
 555                                 //
 556
 557                                 TString file(currFile);
 558                                 fXsec = 0;
 559                                 fTrials = 1;
 560
 561                                 if(file.Contains("root_archive.zip#")){
 562                                                                 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 563                                                                 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 564                                                                 file.Replace(pos+1,20,"");
 565                                 }
 566                                 else {
 567                                                                 // not an archive take the basename....
 568                                                                 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 569                                 }
 570                                 //  Printf("%s",file.Data());
 571
 572
 573                                 TFile *filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really
 574                                 //test the existance of a file in a archive so we have to lvie with open error message from root
 575                                                                 if(!filexsec){
 576                                                                                                 // next trial fetch the histgram file
 577                                                                                                 filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 578                                                                                                 if(!filexsec){
 579                                                                                                                                 // not a severe condition but inciate that we have no information
 580                                                                                                                                 return kFALSE;
 581                                                                                                 }
 582                                                                                                 else{
 583                                                                                                                                 // find the tlist we want to be independtent of the name so use the Tkey
 584                                                                                                                                 TKey* key = (TKey*)filexsec->GetListOfKeys()->At(0);
 585                                                                                                                                 if(!key){
 586                                                                                                                                                                 filexsec->Close();
 587                                                                                                                                                                 return kFALSE;
 588                                                                                                                                 }
 589                                                                                                                                 TList *list = dynamic_cast<TList*>(key->ReadObj());
 590                                                                                                                                 if(!list){
 591                                                                                                                                                                 filexsec->Close();
 592                                                                                                                                                                 return kFALSE;
 593                                                                                                                                 }
 594                                                                                                                                 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 595                                                                                                                                 fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 596                                                                                                                                 filexsec->Close();
 597                                                                                                 }
 598                                                                 } // no tree pyxsec.root
 599                                                                 else {
 600                                                                                                 TTree *xtree = (TTree*)filexsec->Get("Xsection");
 601                                                                                                 if(!xtree){
 602                                                                                                                                 filexsec->Close();
 603                                                                                                                                 return kFALSE;
 604                                                                                                 }
 605                                                                                                 UInt_t   ntrials  = 0;
 606                                                                                                 Double_t  xsection  = 0;
 607                                                                                                 xtree->SetBranchAddress("xsection",&xsection);
 608                                                                                                 xtree->SetBranchAddress("ntrials",&ntrials);
 609                                                                                                 xtree->GetEntry(0);
 610                                                                                                 fTrials = ntrials;
 611                                                                                                 fXsec = xsection;
 612                                                                                                 filexsec->Close();
 613                                                                 }
 614                                 return kTRUE;
 615 }
 616
 617
 618 //Float_t AliAnalysisTaskCheckSingleTrackJetRejection::GetTotalEvents(const char* currFile){
 619 //                              Float_t totalevent;
 620 //                              TString file_es(currFile);
 621 //                              if(file_es.Contains("root_archive.zip#")){
 622 //                                                              Ssiz_t pos1 = file_es.Index("root_archive",12,TString::kExact);
 623 //                                                              Ssiz_t pos = file_es.Index("#",1,pos1,TString::kExact);
 624 //                                                              file_es.Replace(pos+1,20,"");
 625 //                              }
 626 //                              else {
 627 //                                                              // not an archive take the basename....
 628 //                                                              file_es.ReplaceAll(gSystem->BaseName(file_es.Data()),"");
 629 //                              }
 630 //
 631 //                              TString cAdd = "";
 632 //                              cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
 633 //                              cAdd += Form("B%d",(int)BackM);
 634 //                              cAdd += Form("_Filter%05d",Filtermask);
 635 //                              cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
 636 //                              cAdd += Form("_Skip%02d",SkipCone);
 637 //                              TString Dirname,Listname;
 638 //                              Dirname  = Form("PWG4_cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
 639 //                              Listname = Form("pwg4cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
 640 //
 641 //                              TFile *feventstat = TFile::Open(Form("%s%s",file_es.Data(),"JetTasksOutput.root"));
 642 //                              gROOT->Cd(Dirname.Data());
 643 //                              TList *templist     = (TList*)gROOT->FindObject(Listname.Data());
 644 //                              TH1F* temphist = (TH1F*)templist->FindObject("fh1Trials");
 645 //                              totalevent = temphist->Integral();
 646 //                              //cout<<temphist->Integral()<<endl;
 647 //                              delete feventstat;
 648 //                              return totalevent;
 649 //
 650 //}
 651
 652