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