JETAN/AliAnalysisTaskJetCluster.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom3.h>
  25 #include <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TRefArray.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetCluster.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliJetFinder.h"
  43 #include "AliJetHeader.h"
  44 #include "AliJetReader.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODTrack.h"
  49 #include "AliAODJet.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenPythiaEventHeader.h"
  55 #include "AliJetKineReaderHeader.h"
  56 #include "AliGenCocktailEventHeader.h"
  57 #include "AliInputEventHandler.h"
  58 #include "AliAODJetEventBackground.h"
  59
  60 #include "fastjet/PseudoJet.hh"
  61 #include "fastjet/ClusterSequenceArea.hh"
  62 #include "fastjet/AreaDefinition.hh"
  63 #include "fastjet/JetDefinition.hh"
  64 // get info on how fastjet was configured
  65 #include "fastjet/config.h"
  66
  67
  68 ClassImp(AliAnalysisTaskJetCluster)
  69
  70 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
  71   delete fRef;
  72   delete fRandom;
  73 }
  74
  75 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
  76   fAOD(0x0),
  77   fAODExtension(0x0),
  78   fRef(new TRefArray),
  79   fUseAODTrackInput(kFALSE),
  80   fUseAODMCInput(kFALSE),
  81   fUseGlobalSelection(kFALSE),
  82   fUseBackgroundCalc(kFALSE),
  83   fFilterMask(0),
  84   fTrackTypeRec(kTrackUndef),
  85   fTrackTypeGen(kTrackUndef),
  86   fNSkipLeadingRan(0),
  87   fNRandomCones(10),
  88   fAvgTrials(1),
  89   fExternalWeight(1),
  90   fRecEtaWindow(0.5),
  91   fTrackPtCut(0.),
  92   fJetOutputMinPt(1),
  93   fJetTriggerPtCut(0),
  94   fCentCutUp(0),
  95   fCentCutLo(0),
  96   fNonStdBranch(""),
  97   fBackgroundBranch(""),
  98   fNonStdFile(""),
  99   fRparam(1.0),
 100   fAlgorithm(fastjet::kt_algorithm),
 101   fStrategy(fastjet::Best),
 102   fRecombScheme(fastjet::BIpt_scheme),
 103   fAreaType(fastjet::active_area),
 104   fGhostArea(0.01),
 105   fActiveAreaRepeats(1),
 106   fGhostEtamax(1.5),
 107   fRandom(0),
 108   fh1Xsec(0x0),
 109   fh1Trials(0x0),
 110   fh1PtHard(0x0),
 111   fh1PtHardNoW(0x0),
 112   fh1PtHardTrials(0x0),
 113   fh1NJetsRec(0x0),
 114   fh1NConstRec(0x0),
 115   fh1NConstLeadingRec(0x0),
 116   fh1PtJetsRecIn(0x0),
 117   fh1PtJetsLeadingRecIn(0x0),
 118   fh1PtJetConstRec(0x0),
 119   fh1PtJetConstLeadingRec(0x0),
 120   fh1PtTracksRecIn(0x0),
 121   fh1PtTracksLeadingRecIn(0x0),
 122   fh1NJetsRecRan(0x0),
 123   fh1NConstRecRan(0x0),
 124   fh1PtJetsLeadingRecInRan(0x0),
 125   fh1NConstLeadingRecRan(0x0),
 126   fh1PtJetsRecInRan(0x0),
 127   fh1PtTracksGenIn(0x0),
 128   fh1Nch(0x0),
 129   fh1CentralityPhySel(0x0),
 130   fh1Centrality(0x0),
 131   fh1CentralitySelect(0x0),
 132   fh1ZPhySel(0x0),
 133   fh1Z(0x0),
 134   fh1ZSelect(0x0),
 135   fh2NRecJetsPt(0x0),
 136   fh2NRecTracksPt(0x0),
 137   fh2NConstPt(0x0),
 138   fh2NConstLeadingPt(0x0),
 139   fh2JetPhiEta(0x0),
 140   fh2LeadingJetPhiEta(0x0),
 141   fh2JetEtaPt(0x0),
 142   fh2LeadingJetEtaPt(0x0),
 143   fh2TrackEtaPt(0x0),
 144   fh2LeadingTrackEtaPt(0x0),
 145   fh2JetsLeadingPhiEta(0x0),
 146   fh2JetsLeadingPhiPt(0x0),
 147   fh2TracksLeadingPhiEta(0x0),
 148   fh2TracksLeadingPhiPt(0x0),
 149   fh2TracksLeadingJetPhiPt(0x0),
 150   fh2JetsLeadingPhiPtW(0x0),
 151   fh2TracksLeadingPhiPtW(0x0),
 152   fh2TracksLeadingJetPhiPtW(0x0),
 153   fh2NRecJetsPtRan(0x0),
 154   fh2NConstPtRan(0x0),
 155   fh2NConstLeadingPtRan(0x0),
 156   fh2PtNch(0x0),
 157   fh2PtNchRan(0x0),
 158   fh2PtNchN(0x0),
 159   fh2PtNchNRan(0x0),
 160   fh2TracksLeadingJetPhiPtRan(0x0),
 161   fh2TracksLeadingJetPhiPtWRan(0x0),
 162   fHistList(0x0)
 163 {
 164   for(int i = 0;i<3;i++){
 165     fh1BiARandomCones[i] = 0;
 166     fh1BiARandomConesRan[i] = 0;
 167   }
 168   for(int i = 0;i<kMaxCent;i++){
 169     fh2JetsLeadingPhiPtC[i] = 0;
 170     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 171     fh2TracksLeadingJetPhiPtC[i] = 0;
 172     fh2TracksLeadingJetPhiPtWC[i] = 0;
 173   }
 174 }
 175
 176 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
 177   AliAnalysisTaskSE(name),
 178   fAOD(0x0),
 179   fAODExtension(0x0),
 180   fRef(new TRefArray),
 181   fUseAODTrackInput(kFALSE),
 182   fUseAODMCInput(kFALSE),
 183   fUseGlobalSelection(kFALSE),
 184   fUseBackgroundCalc(kFALSE),
 185   fFilterMask(0),
 186   fTrackTypeRec(kTrackUndef),
 187   fTrackTypeGen(kTrackUndef),
 188   fNSkipLeadingRan(0),
 189   fNRandomCones(5),
 190   fAvgTrials(1),
 191   fExternalWeight(1),
 192   fRecEtaWindow(0.5),
 193   fTrackPtCut(0.),
 194   fJetOutputMinPt(1),
 195   fJetTriggerPtCut(0),
 196   fCentCutUp(0),
 197   fCentCutLo(0),
 198   fNonStdBranch(""),
 199   fBackgroundBranch(""),
 200   fNonStdFile(""),
 201   fRparam(1.0),
 202   fAlgorithm(fastjet::kt_algorithm),
 203   fStrategy(fastjet::Best),
 204   fRecombScheme(fastjet::BIpt_scheme),
 205   fAreaType(fastjet::active_area),
 206   fGhostArea(0.01),
 207   fActiveAreaRepeats(1),
 208   fGhostEtamax(1.5),
 209   fRandom(0),
 210   fh1Xsec(0x0),
 211   fh1Trials(0x0),
 212   fh1PtHard(0x0),
 213   fh1PtHardNoW(0x0),
 214   fh1PtHardTrials(0x0),
 215   fh1NJetsRec(0x0),
 216   fh1NConstRec(0x0),
 217   fh1NConstLeadingRec(0x0),
 218   fh1PtJetsRecIn(0x0),
 219   fh1PtJetsLeadingRecIn(0x0),
 220   fh1PtJetConstRec(0x0),
 221   fh1PtJetConstLeadingRec(0x0),
 222   fh1PtTracksRecIn(0x0),
 223   fh1PtTracksLeadingRecIn(0x0),
 224   fh1NJetsRecRan(0x0),
 225   fh1NConstRecRan(0x0),
 226   fh1PtJetsLeadingRecInRan(0x0),
 227   fh1NConstLeadingRecRan(0x0),
 228   fh1PtJetsRecInRan(0x0),
 229   fh1PtTracksGenIn(0x0),
 230   fh1Nch(0x0),
 231   fh1CentralityPhySel(0x0),
 232   fh1Centrality(0x0),
 233   fh1CentralitySelect(0x0),
 234   fh1ZPhySel(0x0),
 235   fh1Z(0x0),
 236   fh1ZSelect(0x0),
 237   fh2NRecJetsPt(0x0),
 238   fh2NRecTracksPt(0x0),
 239   fh2NConstPt(0x0),
 240   fh2NConstLeadingPt(0x0),
 241   fh2JetPhiEta(0x0),
 242   fh2LeadingJetPhiEta(0x0),
 243   fh2JetEtaPt(0x0),
 244   fh2LeadingJetEtaPt(0x0),
 245   fh2TrackEtaPt(0x0),
 246   fh2LeadingTrackEtaPt(0x0),
 247   fh2JetsLeadingPhiEta(0x0),
 248   fh2JetsLeadingPhiPt(0x0),
 249   fh2TracksLeadingPhiEta(0x0),
 250   fh2TracksLeadingPhiPt(0x0),
 251   fh2TracksLeadingJetPhiPt(0x0),
 252   fh2JetsLeadingPhiPtW(0x0),
 253   fh2TracksLeadingPhiPtW(0x0),
 254   fh2TracksLeadingJetPhiPtW(0x0),
 255   fh2NRecJetsPtRan(0x0),
 256   fh2NConstPtRan(0x0),
 257   fh2NConstLeadingPtRan(0x0),
 258   fh2PtNch(0x0),
 259   fh2PtNchRan(0x0),
 260   fh2PtNchN(0x0),
 261   fh2PtNchNRan(0x0),
 262   fh2TracksLeadingJetPhiPtRan(0x0),
 263   fh2TracksLeadingJetPhiPtWRan(0x0),
 264   fHistList(0x0)
 265 {
 266   for(int i = 0;i<3;i++){
 267     fh1BiARandomCones[i] = 0;
 268     fh1BiARandomConesRan[i] = 0;
 269   }
 270   for(int i = 0;i<kMaxCent;i++){
 271     fh2JetsLeadingPhiPtC[i] = 0;
 272     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
 273     fh2TracksLeadingJetPhiPtC[i] = 0;
 274     fh2TracksLeadingJetPhiPtWC[i] = 0;
 275   }
 276  DefineOutput(1, TList::Class());
 277 }
 278
 279
 280
 281 Bool_t AliAnalysisTaskJetCluster::Notify()
 282 {
 283   //
 284   // Implemented Notify() to read the cross sections
 285   // and number of trials from pyxsec.root
 286   //
 287   return kTRUE;
 288 }
 289
 290 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 291 {
 292
 293   //
 294   // Create the output container
 295   //
 296
 297   fRandom = new TRandom3(0);
 298
 299
 300   // Connect the AOD
 301
 302
 303   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 304
 305
 306
 307   if(fNonStdBranch.Length()!=0)
 308     {
 309       // only create the output branch if we have a name
 310       // Create a new branch for jets...
 311       //  -> cleared in the UserExec....
 312       // here we can also have the case that the brnaches are written to a separate file
 313
 314       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
 315       tca->SetName(fNonStdBranch.Data());
 316       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
 317
 318
 319       TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
 320       tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
 321       AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
 322
 323
 324       if(fUseBackgroundCalc){
 325         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
 326           AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
 327           evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
 328           AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
 329         }
 330         // create the branch for the random cones with the same R
 331         TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 332         if(!AODEvent()->FindListObject(cName.Data())){
 333           TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
 334           tcaC->SetName(cName.Data());
 335           AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
 336         }
 337
 338         // create the branch with the random for the random cones on the random event
 339         cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 340         if(!AODEvent()->FindListObject(cName.Data())){
 341           TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
 342           tcaCran->SetName(cName.Data());
 343           AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
 344         }
 345       }
 346
 347       if(fNonStdFile.Length()!=0){
 348         //
 349         // case that we have an AOD extension we need to fetch the jets from the extended output
 350         // we identifay the extension aod event by looking for the branchname
 351         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 352         TObjArray* extArray = aodH->GetExtensions();
 353         if (extArray) {
 354           TIter next(extArray);
 355           while ((fAODExtension=(AliAODExtension*)next())){
 356             TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
 357             if(fDebug>10){
 358               Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
 359               fAODExtension->GetAOD()->Dump();
 360             }
 361             if(obj){
 362               if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
 363               break;
 364             }
 365             fAODExtension = 0;
 366           }
 367         }
 368       }
 369     }
 370
 371
 372   if(!fHistList)fHistList = new TList();
 373   fHistList->SetOwner();
 374
 375   Bool_t oldStatus = TH1::AddDirectoryStatus();
 376   TH1::AddDirectory(kFALSE);
 377
 378   //
 379   //  Histogram
 380
 381   const Int_t nBinPt = 200;
 382   Double_t binLimitsPt[nBinPt+1];
 383   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 384     if(iPt == 0){
 385       binLimitsPt[iPt] = 0.0;
 386     }
 387     else {// 1.0
 388       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
 389     }
 390   }
 391
 392   const Int_t nBinPhi = 90;
 393   Double_t binLimitsPhi[nBinPhi+1];
 394   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
 395     if(iPhi==0){
 396       binLimitsPhi[iPhi] = -1.*TMath::Pi();
 397     }
 398     else{
 399       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
 400     }
 401   }
 402
 403
 404
 405   const Int_t nBinEta = 40;
 406   Double_t binLimitsEta[nBinEta+1];
 407   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
 408     if(iEta==0){
 409       binLimitsEta[iEta] = -2.0;
 410     }
 411     else{
 412       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
 413     }
 414   }
 415
 416   const int nChMax = 4000;
 417
 418   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 419   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 420
 421   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 422   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 423
 424
 425   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
 426   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
 427
 428   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
 429   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
 430   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
 431   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
 432
 433
 434   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 435   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 436   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 437
 438   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 439   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 440   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 441   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 442   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 443   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 444   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 445   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 446   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
 447   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
 448
 449   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
 450   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
 451   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
 452
 453   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
 454   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
 455   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
 456
 457   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 458   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
 459   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
 460   //
 461
 462
 463   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 464   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 465   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 466   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
 467
 468   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 469   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 470   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 471   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
 472
 473
 474
 475   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
 476                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 477   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
 478                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
 479
 480   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
 481                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 482   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
 483                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 484
 485   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
 486                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 487   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
 488                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
 489
 490
 491
 492   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
 493                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 494   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 495                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 496   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
 497                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
 498   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 499                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 500   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 501                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 502   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 503                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 504
 505   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
 506                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 507   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
 508                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 509
 510   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 511                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 512   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
 513                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 514
 515
 516   if(fUseBackgroundCalc){
 517     for(int i = 0;i<3;i++){
 518       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 519       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
 520     }
 521   }
 522   for(int i = 0;i < kMaxCent;i++){
 523     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
 524     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
 525     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
 526     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
 527   }
 528
 529   const Int_t saveLevel = 3; // large save level more histos
 530   if(saveLevel>0){
 531     fHistList->Add(fh1Xsec);
 532     fHistList->Add(fh1Trials);
 533
 534     fHistList->Add(fh1NJetsRec);
 535     fHistList->Add(fh1NConstRec);
 536     fHistList->Add(fh1NConstLeadingRec);
 537     fHistList->Add(fh1PtJetsRecIn);
 538     fHistList->Add(fh1PtJetsLeadingRecIn);
 539     fHistList->Add(fh1PtTracksRecIn);
 540     fHistList->Add(fh1PtTracksLeadingRecIn);
 541     fHistList->Add(fh1PtJetConstRec);
 542     fHistList->Add(fh1PtJetConstLeadingRec);
 543     fHistList->Add(fh1NJetsRecRan);
 544     fHistList->Add(fh1NConstRecRan);
 545     fHistList->Add(fh1PtJetsLeadingRecInRan);
 546     fHistList->Add(fh1NConstLeadingRecRan);
 547     fHistList->Add(fh1PtJetsRecInRan);
 548     fHistList->Add(fh1Nch);
 549     fHistList->Add(fh1Centrality);
 550     fHistList->Add(fh1CentralitySelect);
 551     fHistList->Add(fh1CentralityPhySel);
 552     fHistList->Add(fh1Z);
 553     fHistList->Add(fh1ZSelect);
 554     fHistList->Add(fh1ZPhySel);
 555     if(fUseBackgroundCalc){
 556       for(int i = 0;i<3;i++){
 557         fHistList->Add(fh1BiARandomCones[i]);
 558         fHistList->Add(fh1BiARandomConesRan[i]);
 559       }
 560     }
 561   for(int i = 0;i < kMaxCent;i++){
 562     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
 563     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
 564     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
 565     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
 566   }
 567
 568     fHistList->Add(fh2NRecJetsPt);
 569     fHistList->Add(fh2NRecTracksPt);
 570     fHistList->Add(fh2NConstPt);
 571     fHistList->Add(fh2NConstLeadingPt);
 572     fHistList->Add(fh2PtNch);
 573     fHistList->Add(fh2PtNchRan);
 574     fHistList->Add(fh2PtNchN);
 575     fHistList->Add(fh2PtNchNRan);
 576     fHistList->Add(fh2JetPhiEta);
 577     fHistList->Add(fh2LeadingJetPhiEta);
 578     fHistList->Add(fh2JetEtaPt);
 579     fHistList->Add(fh2LeadingJetEtaPt);
 580     fHistList->Add(fh2TrackEtaPt);
 581     fHistList->Add(fh2LeadingTrackEtaPt);
 582     fHistList->Add(fh2JetsLeadingPhiEta );
 583     fHistList->Add(fh2JetsLeadingPhiPt);
 584     fHistList->Add(fh2TracksLeadingPhiEta);
 585     fHistList->Add(fh2TracksLeadingPhiPt);
 586     fHistList->Add(fh2TracksLeadingJetPhiPt);
 587     fHistList->Add(fh2JetsLeadingPhiPtW);
 588     fHistList->Add(fh2TracksLeadingPhiPtW);
 589     fHistList->Add(fh2TracksLeadingJetPhiPtW);
 590     fHistList->Add(fh2NRecJetsPtRan);
 591     fHistList->Add(fh2NConstPtRan);
 592     fHistList->Add(fh2NConstLeadingPtRan);
 593     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
 594     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
 595     }
 596
 597   // =========== Switch on Sumw2 for all histos ===========
 598   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 599     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 600     if (h1){
 601       h1->Sumw2();
 602       continue;
 603     }
 604     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 605     if(hn)hn->Sumw2();
 606   }
 607   TH1::AddDirectory(oldStatus);
 608 }
 609
 610 void AliAnalysisTaskJetCluster::Init()
 611 {
 612   //
 613   // Initialization
 614   //
 615
 616   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
 617
 618 }
 619
 620 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 621 {
 622
 623   if(fUseGlobalSelection){
 624     // no selection by the service task, we continue
 625     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
 626     PostData(1, fHistList);
 627     return;
 628   }
 629
 630
 631
 632   // handle and reset the output jet branch
 633   // only need this once
 634   TClonesArray* jarray = 0;
 635   TClonesArray* jarrayran = 0;
 636   TClonesArray* rConeArray = 0;
 637   TClonesArray* rConeArrayRan = 0;
 638   AliAODJetEventBackground*  evBkg = 0;
 639   if(fNonStdBranch.Length()!=0) {
 640     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
 641     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
 642     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
 643     if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 644     if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
 645     if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
 646
 647     if(fUseBackgroundCalc){
 648       if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 649       if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
 650       if(evBkg)evBkg->Reset();
 651
 652       TString cName = Form("%sRandomCone",fNonStdBranch.Data());
 653       if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
 654       if(!rConeArray)rConeArray =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
 655       if(rConeArray)rConeArray->Delete();
 656
 657       cName = Form("%sRandomCone_random",fNonStdBranch.Data());
 658       if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
 659       if(!rConeArrayRan)rConeArrayRan =  (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
 660       if(rConeArrayRan)rConeArrayRan->Delete();
 661     }
 662   }
 663
 664   static AliAODJetEventBackground* externalBackground = 0;
 665   if(!externalBackground&&fBackgroundBranch.Length()){
 666     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
 667   }
 668   //
 669   // Execute analysis for current event
 670   //
 671   AliESDEvent *fESD = 0;
 672   if(fUseAODTrackInput){
 673     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 674     if(!fAOD){
 675       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
 676       return;
 677     }
 678     // fethc the header
 679   }
 680   else{
 681     //  assume that the AOD is in the general output...
 682     fAOD  = AODEvent();
 683     if(!fAOD){
 684       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 685       return;
 686     }
 687     if(fDebug>0){
 688       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 689     }
 690   }
 691
 692   Bool_t selectEvent =  false;
 693   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 694
 695   Float_t cent = 0;
 696   Float_t zVtx  = 0;
 697   Int_t cenClass = -1;
 698   if(fAOD){
 699     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
 700     TString vtxTitle(vtxAOD->GetTitle());
 701     zVtx = vtxAOD->GetZ();
 702
 703     cent = fAOD->GetHeader()->GetCentrality();
 704     if(cent<10)cenClass = 0;
 705     else if(cent<30)cenClass = 1;
 706     else if(cent<50)cenClass = 2;
 707     else if(cent<80)cenClass = 3;
 708     if(physicsSelection){
 709       fh1CentralityPhySel->Fill(cent);
 710       fh1ZPhySel->Fill(zVtx);
 711     }
 712
 713
 714     if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 715         Float_t yvtx = vtxAOD->GetY();
 716         Float_t xvtx = vtxAOD->GetX();
 717         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 718         if(TMath::Abs(zVtx)<20.&&r2<1.){ // apply vertex cut later on
 719           if(physicsSelection){
 720             selectEvent = true;
 721           }
 722         }
 723     }
 724     if(fCentCutUp>0){
 725       if(cent<fCentCutLo||cent>fCentCutUp){
 726         selectEvent = false;
 727       }
 728     }
 729
 730   }
 731   if(!selectEvent){
 732     PostData(1, fHistList);
 733     return;
 734   }
 735   fh1Centrality->Fill(cent);
 736   fh1Z->Fill(zVtx);
 737   fh1Trials->Fill("#sum{ntrials}",1);
 738
 739
 740   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 741
 742   // ==== General variables needed
 743
 744
 745
 746   // we simply fetch the tracks/mc particles as a list of AliVParticles
 747
 748   TList recParticles;
 749   TList genParticles;
 750
 751   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
 752   Float_t nCh = recParticles.GetEntries();
 753   fh1Nch->Fill(nCh);
 754   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
 755   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
 756   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
 757
 758   // find the jets....
 759
 760   vector<fastjet::PseudoJet> inputParticlesRec;
 761   vector<fastjet::PseudoJet> inputParticlesRecRan;
 762
 763   // Generate the random cones
 764
 765   AliAODJet vTmpRan(1,0,0,1);
 766   for(int i = 0; i < recParticles.GetEntries(); i++){
 767     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 768     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
 769     // we take total momentum here
 770     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
 771     jInp.set_user_index(i);
 772     inputParticlesRec.push_back(jInp);
 773
 774     // the randomized input changes eta and phi, but keeps the p_T
 775     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 776       Double_t pT = vp->Pt();
 777       Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
 778       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
 779
 780       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
 781       Double_t pZ = pT/TMath::Tan(theta);
 782
 783       Double_t pX = pT * TMath::Cos(phi);
 784       Double_t pY = pT * TMath::Sin(phi);
 785       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ);
 786       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
 787
 788       jInpRan.set_user_index(i);
 789       inputParticlesRecRan.push_back(jInpRan);
 790       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
 791     }
 792
 793     // fill the tref array, only needed when we write out jets
 794     if(jarray){
 795       if(i == 0){
 796         fRef->Delete(); // make sure to delete before placement new...
 797         new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
 798       }
 799       fRef->Add(vp);
 800     }
 801   }// recparticles
 802
 803   if(inputParticlesRec.size()==0){
 804     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 805     PostData(1, fHistList);
 806     return;
 807   }
 808
 809   // run fast jet
 810   // employ setters for these...
 811
 812
 813   // now create the object that holds info about ghosts
 814   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
 815     // reduce CPU time...
 816     fGhostArea = 0.5;
 817     fActiveAreaRepeats = 0;
 818   }
 819
 820  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
 821   fastjet::AreaType areaType =   fastjet::active_area;
 822   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 823   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
 824   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 825   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
 826
 827   //range where to compute background
 828   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
 829   phiMin = 0;
 830   phiMax = 2*TMath::Pi();
 831   rapMax = fGhostEtamax - fRparam;
 832   rapMin = - fGhostEtamax + fRparam;
 833   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
 834
 835
 836
 837   inclusiveJets = clustSeq.inclusive_jets();
 838   sortedJets    = sorted_by_pt(inclusiveJets);
 839
 840   fh1NJetsRec->Fill(sortedJets.size());
 841
 842  // loop over all jets an fill information, first on is the leading jet
 843
 844   Int_t nRecOver = inclusiveJets.size();
 845   Int_t nRec     = inclusiveJets.size();
 846   if(inclusiveJets.size()>0){
 847     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
 848     Double_t area = clustSeq.area(sortedJets[0]);
 849     leadingJet.SetEffArea(area,0);
 850     Float_t pt = leadingJet.Pt();
 851     Int_t nAodOutJets = 0;
 852     Int_t nAodOutTracks = 0;
 853     AliAODJet *aodOutJet = 0;
 854
 855     Int_t iCount = 0;
 856     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
 857       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
 858       while(pt<ptCut&&iCount<nRec){
 859         nRecOver--;
 860         iCount++;
 861         if(iCount<nRec){
 862           pt = sortedJets[iCount].perp();
 863         }
 864       }
 865       if(nRecOver<=0)break;
 866       fh2NRecJetsPt->Fill(ptCut,nRecOver);
 867     }
 868     Float_t phi = leadingJet.Phi();
 869     if(phi<0)phi+=TMath::Pi()*2.;
 870     Float_t eta = leadingJet.Eta();
 871     Float_t pTback = 0;
 872     if(externalBackground){
 873       // carefull has to be filled in a task before
 874       // todo, ReArrange to the botom
 875       pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
 876     }
 877     pt = leadingJet.Pt() - pTback;
 878     // correlation of leading jet with tracks
 879     TIterator *recIter = recParticles.MakeIterator();
 880     recIter->Reset();
 881     AliVParticle *tmpRecTrack = 0;
 882     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
 883       Float_t tmpPt = tmpRecTrack->Pt();
 884       // correlation
 885       Float_t tmpPhi =  tmpRecTrack->Phi();
 886       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
 887       Float_t dPhi = phi - tmpPhi;
 888       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 889       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 890       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
 891       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
 892       if(cenClass>=0){
 893         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
 894         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
 895       }
 896
 897     }
 898
 899
 900
 901     for(int j = 0; j < nRec;j++){
 902       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
 903       aodOutJet = 0;
 904       nAodOutTracks = 0;
 905       Float_t tmpPt = tmpRec.Pt();
 906       Float_t tmpPtBack = 0;
 907       if(externalBackground){
 908         // carefull has to be filled in a task before
 909        // todo, ReArrange to the botom
 910         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
 911       }
 912       tmpPt = tmpPt - tmpPtBack;
 913       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
 914
 915       fh1PtJetsRecIn->Fill(tmpPt);
 916       // Fill Spectra with constituents
 917       vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
 918       fh1NConstRec->Fill(constituents.size());
 919       fh2PtNch->Fill(nCh,tmpPt);
 920       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
 921       fh2NConstPt->Fill(tmpPt,constituents.size());
 922       // loop over constiutents and fill spectrum
 923
 924       // Add the jet information and the track references to the output AOD
 925
 926       if(tmpPt>fJetOutputMinPt&&jarray){
 927         aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
 928         Double_t area1 = clustSeq.area(sortedJets[j]);
 929         aodOutJet->SetEffArea(area1,0);
 930       }
 931
 932
 933       if(fUseBackgroundCalc){
 934
 935         // create a random jet within the acceptance
 936         Double_t etaMax = 0.9 - fRparam;
 937         Int_t nCone = 0;
 938         Int_t nConeRan = 0;
 939         Double_t pTC = 1; // small number
 940         for(int ir = 0;ir < fNRandomCones;ir++){
 941           Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
 942           Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
 943           // massless jet
 944           Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
 945           Double_t pZC = pTC/TMath::Tan(thetaC);
 946           Double_t pXC = pTC * TMath::Cos(phiC);
 947           Double_t pYC = pTC * TMath::Sin(phiC);
 948           Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
 949           AliAODJet tmpRecC (pXC,pYC,pZC, pC);
 950           bool skip = false;
 951           for(int jj = 0; jj < TMath::Min(nRec,2);jj++){
 952             AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
 953             if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
 954               skip = true;
 955               break;
 956             }
 957           }
 958           if(skip)continue;
 959           tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
 960           if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
 961           if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
 962         }// random cones
 963
 964
 965         // loop over the reconstructed particles and add up the pT in the random cones
 966         // maybe better to loop over randomized particles not in the real jets...
 967         // but this by definition brings dow average energy in the whole  event
 968         AliAODJet vTmpRanR(1,0,0,1);
 969         for(int i = 0; i < recParticles.GetEntries(); i++){
 970           AliVParticle *vp = (AliVParticle*)recParticles.At(i);
 971           if(rConeArray){
 972             for(int ir = 0;ir < fNRandomCones;ir++){
 973               AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
 974               if(jC&&jC->DeltaR(vp)<fRparam){
 975                 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
 976               }
 977             }
 978           }// add up energy in cone
 979
 980           // the randomized input changes eta and phi, but keeps the p_T
 981           if(i>=fNSkipLeadingRan){// eventually skip the leading particles
 982             Double_t pTR = vp->Pt();
 983             Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
 984             Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
 985
 986             Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
 987             Double_t pZR = pTR/TMath::Tan(thetaR);
 988
 989             Double_t pXR = pTR * TMath::Cos(phiR);
 990             Double_t pYR = pTR * TMath::Sin(phiR);
 991             Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR);
 992             vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
 993             if(rConeArrayRan){
 994               for(int ir = 0;ir < fNRandomCones;ir++){
 995                 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
 996                 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
 997                   jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
 998                 }
 999               }
1000             }
1001           }
1002         }// loop over recparticles
1003
1004         Float_t jetArea = fRparam*fRparam*TMath::Pi();
1005         for(int ir = 0;ir < fNRandomCones;ir++){
1006           // rescale the momntum vectors for the random cones
1007           if(!rConeArray)continue;
1008           AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
1009           if(rC){
1010             Double_t etaC = rC->Eta();
1011             Double_t phiC = rC->Phi();
1012             // massless jet, unit vector
1013             pTC = rC->ChargedBgEnergy();
1014             if(pTC<=0)pTC = 0.1; // for almost empty events
1015             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1016             Double_t pZC = pTC/TMath::Tan(thetaC);
1017             Double_t pXC = pTC * TMath::Cos(phiC);
1018             Double_t pYC = pTC * TMath::Sin(phiC);
1019             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1020             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1021             rC->SetBgEnergy(0,0);
1022             rC->SetEffArea(jetArea,0);
1023           }
1024           rC = (AliAODJet*)rConeArrayRan->At(ir);
1025           // same wit random
1026           if(rC){
1027             Double_t etaC = rC->Eta();
1028             Double_t phiC = rC->Phi();
1029             // massless jet, unit vector
1030             pTC = rC->ChargedBgEnergy();
1031             if(pTC<=0)pTC = 0.1;// for almost empty events
1032             Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1033             Double_t pZC = pTC/TMath::Tan(thetaC);
1034             Double_t pXC = pTC * TMath::Cos(phiC);
1035             Double_t pYC = pTC * TMath::Sin(phiC);
1036             Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC);
1037             rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1038             rC->SetBgEnergy(0,0);
1039             rC->SetEffArea(jetArea,0);
1040           }
1041         }
1042       }// if(fUseBackgroundCalc
1043
1044      for(unsigned int ic = 0; ic < constituents.size();ic++){
1045        AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1046        fh1PtJetConstRec->Fill(part->Pt());
1047        if(aodOutJet){
1048          aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1049        }
1050        if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1051      }
1052
1053      // correlation
1054      Float_t tmpPhi =  tmpRec.Phi();
1055      Float_t tmpEta =  tmpRec.Eta();
1056      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1057
1058
1059
1060      if(j==0){
1061        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1062        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1063        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1064        fh1NConstLeadingRec->Fill(constituents.size());
1065        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1066        continue;
1067      }
1068      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1069      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1070      Float_t dPhi = phi - tmpPhi;
1071      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1072      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1073      Float_t dEta = eta - tmpRec.Eta();
1074      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1075      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1076      if(cenClass>=0){
1077        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1078        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1079      }
1080      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1081     }
1082    delete recIter;
1083
1084    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1085
1086
1087    if(evBkg){
1088      vector<fastjet::PseudoJet> jets2=sortedJets;
1089      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1090      Double_t bkg1=0;
1091      Double_t sigma1=0.;
1092      Double_t meanarea1=0.;
1093      Double_t bkg2=0;
1094      Double_t sigma2=0.;
1095      Double_t meanarea2=0.;
1096
1097      clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
1098      evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
1099
1100      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1101      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1102
1103      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1104      evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
1105      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1106      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1107
1108    }
1109   }
1110
1111
1112
1113
1114
1115   // fill track information
1116   Int_t nTrackOver = recParticles.GetSize();
1117   // do the same for tracks and jets
1118
1119   if(nTrackOver>0){
1120    TIterator *recIter = recParticles.MakeIterator();
1121    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1122    Float_t pt = tmpRec->Pt();
1123
1124    //    Printf("Leading track p_t %3.3E",pt);
1125    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1126      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1127      while(pt<ptCut&&tmpRec){
1128        nTrackOver--;
1129        tmpRec = (AliVParticle*)(recIter->Next());
1130        if(tmpRec){
1131          pt = tmpRec->Pt();
1132        }
1133      }
1134      if(nTrackOver<=0)break;
1135      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1136    }
1137
1138    recIter->Reset();
1139    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1140    Float_t phi = leading->Phi();
1141    if(phi<0)phi+=TMath::Pi()*2.;
1142    Float_t eta = leading->Eta();
1143    pt  = leading->Pt();
1144    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1145      Float_t tmpPt = tmpRec->Pt();
1146      Float_t tmpEta = tmpRec->Eta();
1147      fh1PtTracksRecIn->Fill(tmpPt);
1148      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1149      if(tmpRec==leading){
1150        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1151        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1152        continue;
1153      }
1154       // correlation
1155      Float_t tmpPhi =  tmpRec->Phi();
1156
1157      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1158      Float_t dPhi = phi - tmpPhi;
1159      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1160      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1161      Float_t dEta = eta - tmpRec->Eta();
1162      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1163      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1164      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1165    }
1166    delete recIter;
1167  }
1168
1169  // find the random jets
1170  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1171  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1172
1173  inclusiveJetsRan = clustSeqRan.inclusive_jets();
1174  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1175
1176  // fill the jet information from random track
1177
1178   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1179
1180  // loop over all jets an fill information, first on is the leading jet
1181
1182  Int_t nRecOverRan = inclusiveJetsRan.size();
1183  Int_t nRecRan     = inclusiveJetsRan.size();
1184  if(inclusiveJetsRan.size()>0){
1185    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1186    Float_t pt = leadingJet.Pt();
1187
1188    Int_t iCount = 0;
1189    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1190      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1191       while(pt<ptCut&&iCount<nRecRan){
1192         nRecOverRan--;
1193         iCount++;
1194         if(iCount<nRecRan){
1195           pt = sortedJetsRan[iCount].perp();
1196         }
1197       }
1198       if(nRecOverRan<=0)break;
1199       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1200     }
1201     Float_t phi = leadingJet.Phi();
1202     if(phi<0)phi+=TMath::Pi()*2.;
1203     pt = leadingJet.Pt();
1204
1205     // correlation of leading jet with random tracks
1206
1207     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1208       {
1209         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1210         // correlation
1211         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1212         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1213         Float_t dPhi = phi - tmpPhi;
1214         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1215         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1216         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1217         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1218       }
1219
1220     Int_t nAodOutJetsRan = 0;
1221      AliAODJet *aodOutJetRan = 0;
1222     for(int j = 0; j < nRecRan;j++){
1223       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1224       Float_t tmpPt = tmpRec.Pt();
1225       fh1PtJetsRecInRan->Fill(tmpPt);
1226       // Fill Spectra with constituents
1227       vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1228       fh1NConstRecRan->Fill(constituents.size());
1229       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1230       fh2PtNchRan->Fill(nCh,tmpPt);
1231       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1232
1233
1234      if(tmpPt>fJetOutputMinPt&&jarrayran){
1235        aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1236        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1237
1238        aodOutJetRan->SetEffArea(arearan,0);    }
1239
1240
1241
1242
1243      for(unsigned int ic = 0; ic < constituents.size();ic++){
1244        // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1245        // fh1PtJetConstRec->Fill(part->Pt());
1246        if(aodOutJetRan){
1247          aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
1248        }
1249      }
1250
1251
1252
1253
1254       // correlation
1255       Float_t tmpPhi =  tmpRec.Phi();
1256       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1257
1258       if(j==0){
1259         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1260         fh1NConstLeadingRecRan->Fill(constituents.size());
1261         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1262         continue;
1263       }
1264     }
1265
1266
1267     if(evBkg){
1268      Double_t bkg3=0.;
1269      Double_t sigma3=0.;
1270      Double_t meanarea3=0.;
1271      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1272      evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
1273      //     float areaRandomCone = rRandomCone2 *TMath::Pi();
1274      /*
1275      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1276      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1277      */
1278     }
1279
1280
1281
1282  }
1283
1284
1285  // do the event selection if activated
1286  if(fJetTriggerPtCut>0){
1287    bool select = false;
1288    Float_t minPt = fJetTriggerPtCut;
1289    /*
1290    // hard coded for now ...
1291    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1292    if(cent<10)minPt = 50;
1293    else if(cent<30)minPt = 42;
1294    else if(cent<50)minPt = 28;
1295    else if(cent<80)minPt = 18;
1296    */
1297    float rho = 0;
1298    if(externalBackground)rho = externalBackground->GetBackground(2);
1299    if(jarray){
1300      for(int i = 0;i < jarray->GetEntriesFast();i++){
1301        AliAODJet *jet = (AliAODJet*)jarray->At(i);
1302        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1303        if(ptSub>=minPt){
1304          select = true;
1305          break;
1306        }
1307      }
1308    }
1309
1310    if(select){
1311      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1312      fh1CentralitySelect->Fill(cent);
1313      fh1ZSelect->Fill(zVtx);
1314      aodH->SetFillAOD(kTRUE);
1315    }
1316  }
1317
1318  if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1319  PostData(1, fHistList);
1320 }
1321
1322 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1323 {
1324 // Terminate analysis
1325 //
1326     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1327 }
1328
1329
1330 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1331
1332   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1333
1334   Int_t iCount = 0;
1335   if(type==kTrackAOD){
1336     AliAODEvent *aod = 0;
1337     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1338     else aod = AODEvent();
1339     if(!aod){
1340       return iCount;
1341     }
1342     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1343       AliAODTrack *tr = aod->GetTrack(it);
1344       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1345       if(TMath::Abs(tr->Eta())>0.9)continue;
1346       if(tr->Pt()<fTrackPtCut)continue;
1347       list->Add(tr);
1348       iCount++;
1349     }
1350   }
1351   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1352     AliMCEvent* mcEvent = MCEvent();
1353     if(!mcEvent)return iCount;
1354     // we want to have alivpartilces so use get track
1355     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1356       if(!mcEvent->IsPhysicalPrimary(it))continue;
1357       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1358       if(type == kTrackKineAll){
1359         if(part->Pt()<fTrackPtCut)continue;
1360         list->Add(part);
1361         iCount++;
1362       }
1363       else if(type == kTrackKineCharged){
1364         if(part->Particle()->GetPDG()->Charge()==0)continue;
1365         if(part->Pt()<fTrackPtCut)continue;
1366         list->Add(part);
1367         iCount++;
1368       }
1369     }
1370   }
1371   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1372     AliAODEvent *aod = 0;
1373     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1374     else aod = AODEvent();
1375     if(!aod)return iCount;
1376     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1377     if(!tca)return iCount;
1378     for(int it = 0;it < tca->GetEntriesFast();++it){
1379       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1380       if(!part->IsPhysicalPrimary())continue;
1381       if(type == kTrackAODMCAll){
1382         if(part->Pt()<fTrackPtCut)continue;
1383         list->Add(part);
1384         iCount++;
1385       }
1386       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1387         if(part->Charge()==0)continue;
1388         if(part->Pt()<fTrackPtCut)continue;
1389         if(kTrackAODMCCharged){
1390           list->Add(part);
1391         }
1392         else {
1393           if(TMath::Abs(part->Eta())>0.9)continue;
1394           list->Add(part);
1395         }
1396         iCount++;
1397       }
1398     }
1399   }// AODMCparticle
1400   list->Sort();
1401   return iCount;
1402 }
1403
1404 /*
1405 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1406   for(int i = 0; i < particles.GetEntries(); i++){
1407     AliVParticle *vp = (AliVParticle*)particles.At(i);
1408     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1409     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1410     jInp.set_user_index(i);
1411     inputParticles.push_back(jInp);
1412   }
1413
1414   return 0;
1415
1416 }
1417 */