PWG4/macros/AddTaskJets.C

   1 AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask); // Common config
   2 AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius = -1);
   3
   4 AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf,Float_t radius = -1,UInt_t filterMask = 0,Float_t ptTrackMin = 0.15,Int_t iBack = 0); // for the new AF
   5 Int_t AddTaskJetsDelta(char *nonStdFile = "",UInt_t filterMask = 0,Bool_t kUseAODMC = kTRUE,UInt_t runFlag = 1|4|32|128|256);
   6 AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask = 0);
   7
   8 Float_t kPtTrackMin = 0.15;
   9 Int_t kBackgroundMode = 0;
  10
  11
  12 AliAnalysisTaskJets *AddTaskJets(UInt_t filterMask ){
  13   // fills the standard "jets" branch in the AOD
  14   // need the ESDFilter to run before, to access the AODtracks
  15   // Tracks selected by the first Filter (1<<0)
  16   // needs to be adapted for different cuts
  17
  18   // UA1 as standard chosen, since it is the most robust and simple JF
  19   // R = 0.4 suffficient to provide accurate jet axis for correlation studies
  20   // energy resolution suffers a little
  21   // Acceptance of jets not limited by the Jet Finder but should be done
  22   // by user to abs(eta) < 0.5
  23
  24   return AddTaskJets("AOD","UA1",0.4,filterMask);
  25
  26 }
  27
  28
  29
  30 Int_t AddTaskJetsDelta(char *nonStdFile,UInt_t filterMask,Bool_t kUseAODMC,UInt_t runFlag){
  31
  32   // Adds a whole set of jet finders  all to be written
  33   // to a delta AOD
  34
  35   // this can in principle be done also with on the fly
  36   // if we have an ESD input jet fidner task does automatically fetch the ouput aod
  37
  38    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  39    if (!mgr) {
  40       ::Error("AddTaskJetsDelta", "No analysis manager to connect to.");
  41       return 0;
  42    }
  43
  44    // Check the analysis type using the event handlers connected to the analysis manager.
  45    //==============================================================================
  46    if (!mgr->GetInputEventHandler()) {
  47       ::Error("AddTaskJetsDelta", "This task requires an input event handler");
  48       return 0;
  49    }
  50
  51   AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();
  52   if (!aodh) {
  53     ::Error("AddTaskJetsDelta", "This task needs an output event handler");
  54     return 0;
  55   }
  56
  57   TString outFile(nonStdFile);
  58
  59   AliAnalysisTaskJets *jetana = 0;
  60   Int_t iCount = 0;
  61
  62   // Jet Fidners Selected by run flag first bit 2^0 second by 2^1 etc
  63   const char *cJF[14]        = {"UA1","UA1","UA1", "CDF", "CDF","DA","DA","SISCONE","SISCONE","FASTJET","FASTJET","FASTKT","FASTKT","UA1LO"};
  64   const Float_t radius[14]   = {  0.4,  0.7,  1.0,  0.4,    0.7, 0.4, 0.7,      0.4,     0.7,      0.4,       0.7,     0.4,     0.7,    0.4};
  65   UInt_t  flag[14]           = {    6,    7,    7,    7,      7,   7,   7,        7,       7,        7,         7,       7,       7,      7};
  66   // flag[5] = 0; // set siscone to 0 for proof mode...
  67   // flag first bit AOD, second bit AODMC2 third bit AODMC2 third (8) bit AOODMC2b (limited acceptance)
  68   // i.e. 7 all, 6 only MC2 and MC
  69   // this stay at three
  70   const char *cReader[4] = {"AOD","AODMC","AODMC2","AODMC2b"};
  71
  72
  73
  74   for(int i = 0; i< 14;i++){
  75     if(!(runFlag&(1<<i)))continue;
  76     if(!kUseAODMC)flag[i]&=1; // switch OFF MC if we do not have it
  77     for(int ib = 0;ib<4;ib++){
  78       if(flag[i]&(1<<ib)){
  79         jetana = AddTaskJets(cReader[ib],cJF[i],radius[i],filterMask);
  80         if(jetana){
  81           char *cRadius = "";
  82           if(radius[i]>0)cRadius = Form("%02d",(int)((radius[i]+0.01)*10.)); // add an offset beacuse of precision
  83           //      jetana->SetNonStdBranch(Form("jets%s_%s%s",cReader[ib],cJF[i],cRadius)); // done in addtask jets
  84           if(outFile.Length()>0)jetana->SetNonStdOutputFile(outFile.Data());
  85           iCount++;
  86         }
  87       }
  88     }
  89   }
  90
  91   Printf("Added %d JetFinders",iCount);
  92 }
  93
  94
  95
  96
  97
  98 AliAnalysisTaskJets *AddTaskJets(Char_t *jr, Char_t *jf, Float_t radius,UInt_t filterMask,float ptTrackMin,int iBack)
  99 {
 100   // Creates a jet finder task, configures it and adds it to the analysis manager.
 101
 102   kPtTrackMin = ptTrackMin;
 103   kBackgroundMode = iBack;
 104   // Get the pointer to the existing analysis manager via the static access method.
 105   //==============================================================================
 106   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 107   if (!mgr) {
 108     ::Error("AddTaskJets", "No analysis manager to connect to.");
 109     return NULL;
 110   }
 111
 112   // Check the analysis type using the event handlers connected to the analysis manager.
 113   //==============================================================================
 114   if (!mgr->GetInputEventHandler()) {
 115     ::Error("AddTaskJets", "This task requires an input event handler");
 116     return NULL;
 117    }
 118
 119   AliAODHandler *aodh = (AliAODHandler*)mgr->GetOutputEventHandler();
 120   if (!aodh) {
 121     ::Error("AddTaskJets", "This task needs an output event handler");
 122     return NULL;
 123   }
 124
 125
 126   // Create the task and configure it.
 127   //===========================================================================
 128   AliAnalysisTaskJets *jetana;
 129    AliJetReader *er = CreateJetReader(jr,filterMask);
 130    // Define jet header and jet finder
 131    AliJetFinder *jetFinder = CreateJetFinder(jf,radius);
 132
 133    if (jetFinder){
 134        if (er) jetFinder->SetJetReader(er);
 135    }
 136
 137
 138
 139    TString cAdd = "";
 140    cAdd += Form("%02d",(int)((radius+0.01)*10.));
 141    cAdd += Form("_B%d",(int)((kBackgroundMode)));
 142    cAdd += Form("_Filter%05d",filterMask);
 143    cAdd += Form("_Cut%05d",(int)((1000.*kPtTrackMin)));
 144    Printf("%s",cAdd.Data());
 145    jetana = new AliAnalysisTaskJets(Form("JetAnalysis%s_%s%s",jr,jf,cAdd.Data()));
 146
 147
 148    TString c_jr(jr);
 149    c_jr.ToLower();
 150    TString c_jf(jf);
 151    c_jf.ToLower();
 152
 153    /*
 154    if(c_jf.CompareTo("ua1")==0 && TMath::Abs(radius-0.4) < 0.01 && c_jr.CompareTo("aod") == 0){
 155      // do nothing, this is the standard jet finder R = 0.4, UA1 on AOD
 156    }
 157    else{
 158
 159
 160  }
 161    */
 162    TString bName =  Form("jets%s_%s%s",jr,jf,cAdd.Data());
 163    jetana->SetNonStdBranch(bName.Data());
 164    Printf("Set jet branchname \"%s\"",bName.Data());
 165
 166
 167    AliAnalysisDataContainer *cout_jet = mgr->CreateContainer(
 168                                                              Form("jethist_%s_%s%s",
 169                                                                   c_jr.Data(),
 170                                                                   c_jf.Data(),
 171                                                                   cAdd.Data()),
 172                                                              TList::Class(),
 173                                                              AliAnalysisManager::kOutputContainer,
 174                                                              Form("%s:PWG4_jethist_%s_%s%s",AliAnalysisManager::GetCommonFileName(),
 175                                                                   c_jr.Data(),
 176                                                                   c_jf.Data(),
 177                                                                   cAdd.Data()));
 178
 179    // Connect jet finder to task.
 180    jetana->SetJetFinder(jetFinder);
 181    jetana->SetConfigFile("");
 182    jetana->SetDebugLevel(2);
 183    if(TMath::Abs((radius-0.4))< 0.02&&c_jf.Contains("fastjet")){
 184      //     jetana->SetFilterPt(20.);
 185    }
 186
 187
 188    mgr->AddTask(jetana);
 189
 190
 191    // Create ONLY the output containers for the data produced by the task.
 192    // Get and connect other common input/output containers via the manager as below
 193    //==============================================================================
 194    mgr->ConnectInput  (jetana, 0, mgr->GetCommonInputContainer());
 195 // AOD output slot will be used in a different way in future
 196    mgr->ConnectOutput (jetana, 0, mgr->GetCommonOutputContainer());
 197    mgr->ConnectOutput (jetana, 1, cout_jet);
 198
 199    return jetana;
 200 }
 201
 202 AliJetFinder *CreateJetFinder(Char_t *jf,Float_t radius){
 203   AliJetFinder *jetFinder = 0;
 204
 205   switch (jf) {
 206   case "CDF":
 207     AliCdfJetHeader *jh = new AliCdfJetHeader();
 208     jh->SetRadius(0.7);
 209     if(radius>0)jh->SetRadius(radius);
 210     jh->SetAODwrite(kTRUE);
 211     jh->SetAODtracksWrite(kTRUE);
 212     //    jh->SetDebugCDF(1);
 213     jetFinder = new AliCdfJetFinder();
 214     if (jh) jetFinder->SetJetHeader(jh);
 215     break;
 216
 217   case "DA":
 218     AliDAJetHeader *jh=new AliDAJetHeader();
 219     jh->SetComment("DA jet code with default parameters");
 220     jh->SelectJets(kTRUE);
 221 //    jh->SetNeff(200);
 222 //    jh->SetEtaEff(2.2);
 223     if(radius>0)jh->SetRadius(radius);
 224     jh->SetEtMin(5.);
 225     jetFinder = new AliDAJetFinder();
 226     if (jh) jetFinder->SetJetHeader(jh);
 227     break;
 228
 229   case "FASTJET":
 230     // DEFAULT is ANTI KT
 231     AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
 232     jh->SetRparam(0.4); // setup parameters
 233     Double_t rBkg = 0.6;
 234     if(radius>0)jh->SetRparam(radius);
 235     jh->SetAlgorithm(2); // antikt from fastjet/JetDefinition.hh
 236     /*
 237       $FASTJET/include/fastjet/JetDefinition.hh
 238       enum JetAlgorithm {kt_algorithm, cambridge_algorithm,
 239       antikt_algorithm, genkt_algorithm,
 240       ee_kt_algorithm, ee_genkt_algorithm, ...};
 241
 242      */
 243     jetFinder = new AliFastJetFinder();
 244     jh->SetPtMin(1);
 245     jh->SetRparamBkg(rBkg);
 246     jh->SetGhostEtaMax(0.9+rBkg);
 247     if (jh) jetFinder->SetJetHeader(jh);
 248     break;
 249
 250   case "FASTKT":
 251     AliFastJetHeaderV1 *jh = new AliFastJetHeaderV1();
 252     jh->SetRparam(0.4); // setup parameters
 253     Double_t rBkg = 0.6;
 254     if(radius>0)jh->SetRparam(radius);
 255     jh->SetAlgorithm(0); // kt from fastjet/JetDefinition.hh
 256     jh->SetPtMin(1);
 257     jh->SetRparamBkg(rBkg);
 258     jh->SetGhostEtaMax(0.9+rBkg);
 259     jetFinder = new AliFastJetFinder();
 260     if (jh) jetFinder->SetJetHeader(jh);
 261     break;
 262
 263   case "SISCONE":
 264     AliSISConeJetHeader * jh = new AliSISConeJetHeader();
 265
 266     jh->SetJetEtaMax(1.5);
 267     jh->SetJetEtaMin(-1.5);
 268
 269     //siscone parameters
 270     jh->SetConeRadius(0.4);                   // default cone radius
 271     if(radius>0)jh->SetConeRadius(radius);   // cone radius
 272
 273     jh->SetOverlapThreshold(0.75);            // overlap parameter, between 0 and 1 excluded!! 0.75 value is advised
 274     jh->SetPtProtojetMin(0);                  // pt min of protojets
 275     jh->SetMinJetPt(5);                      // Ptmin of jets (GeV)
 276
 277     //do you want to subtract BG (0 = no, 1 = yes)
 278     jh->SetBGMode(0); // if 1 set also the radius for the background determination..
 279
 280     //for background
 281     jh->SetRapRange( -0.9, 0.9);              // rapidity range for subtracting background must be < ghostmaxrap-0.95*R
 282     jh->SetPhiRange(0 , 2*TMath::Pi());       // phi range for subtracting background
 283
 284     //to determine jets area
 285     jh->SetBGAlgorithm(1);                    // algorithm for rho calculation : 0 = kT, 1 = Cambridge
 286     jh->SetGhostEtaMax(4);                    // eta max where ghosts are generated
 287     jh->SetGhostArea(0.05);                   // area of a ghost
 288     jh->SetMeanGhostKt(1e-100);               // average transverse momentum of the ghosts.
 289     jh->SetAreaTypeNumber(4);                 // from 1 to 5 : 1 = active_area, 2 = active_area_explicit_ghosts, 3 = one_ghost_passive_area, 4 = passive_area, 5 = voronoi_area
 290     jetFinder = new AliSISConeJetFinder();
 291     if (jh) jetFinder->SetJetHeader(jh);
 292     break;
 293
 294   case "UA1":
 295     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
 296     jh->SetComment("UA1 jet code with default parameters");
 297     jh->BackgMode(0);
 298     jh->SetRadius(0.4);
 299     if(radius>0)jh->SetRadius(radius);
 300     jh->SetEtSeed(4.);
 301     jh->SetNAcceptJets(6);
 302     jh->SetLegoNbinPhi(432);
 303     jh->SetLegoNbinEta(274);
 304     jh->SetLegoEtaMin(-2);
 305     jh->SetLegoEtaMax(+2);
 306     jh->SetMinJetEt(5.);
 307     jh->SetJetEtaMax(1.5);
 308     jh->SetJetEtaMin(-1.5);
 309     jh->BackgMode(kBackgroundMode);
 310
 311     jetFinder = new AliUA1JetFinderV1();
 312     if (jh) jetFinder->SetJetHeader(jh);
 313     break;
 314   case "UA1LO":
 315     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
 316     jh->SetComment("UA1 jet code with Lo Pt settings parameters");
 317     jh->BackgMode(0);
 318     jh->SetRadius(0.4);
 319     if(radius>0)jh->SetRadius(radius);
 320     jh->SetEtSeed(1.);
 321     jh->SetNAcceptJets(6);
 322     jh->SetLegoNbinPhi(432);
 323     jh->SetLegoNbinEta(274);
 324     jh->SetLegoEtaMin(-2);
 325     jh->SetLegoEtaMax(+2);
 326     jh->SetMinJetEt(1.);
 327     jh->SetJetEtaMax(1.5);
 328     jh->SetJetEtaMin(-1.5);
 329     jh->BackgMode(kBackgroundMode);
 330
 331     jetFinder = new AliUA1JetFinderV1();
 332     if (jh) jetFinder->SetJetHeader(jh);
 333     break;
 334
 335
 336
 337   case "UA1MC":
 338     AliUA1JetHeaderV1 *jh=new AliUA1JetHeaderV1();
 339     jh->SetComment("UA1 jet code with default MC parameters");
 340     jh->BackgMode(0);
 341     jh->SetRadius(1.0);
 342     if(radius>0)jh->SetRadius(radius);
 343     jh->SetEtSeed(4.);
 344     jh->SetNAcceptJets(6);
 345     jh->SetLegoNbinPhi(432);
 346     jh->SetLegoNbinEta(274);
 347     jh->SetLegoEtaMin(-2);
 348     jh->SetLegoEtaMax(+2);
 349     jh->SetMinJetEt(5.);
 350     jh->SetJetEtaMax(1.5);
 351     jh->SetJetEtaMin(-1.5);
 352     jh->BackgMode(kBackgroundMode);
 353     jetFinder = new AliUA1JetFinderV1();
 354     if (jh) jetFinder->SetJetHeader(jh);
 355     break;
 356   default:
 357     Printf("\n >>>>>>> AddTaskJets Error Wrong jet finder selected\n");
 358     break;
 359   }
 360
 361   return jetFinder;
 362
 363 }
 364
 365 AliJetReader *CreateJetReader(Char_t *jr,UInt_t filterMask){
 366   AliJetReader *er = 0;
 367   Float_t ptCut  = kPtTrackMin ; // cut on track p_T
 368   switch (jr) {
 369   case "MC":
 370     AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
 371     jrh->SetComment("MC full Kinematics");
 372     jrh->SetFastSimTPC(kFALSE);
 373     jrh->SetFastSimEMCAL(kFALSE);
 374     jrh->SetPtCut(ptCut);
 375     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9
 376     // Define reader and set its header
 377     er = new AliJetKineReader();
 378     er->SetReaderHeader(jrh);
 379     break;
 380   case "MC2":
 381     AliJetKineReaderHeader *jrh = new AliJetKineReaderHeader();
 382     jrh->SetComment("MC full Kinematics spearate config charged only");
 383     jrh->SetFastSimTPC(kFALSE);
 384     jrh->SetFastSimEMCAL(kFALSE);
 385     jrh->SetChargedOnly(kTRUE);
 386     jrh->SetPtCut(ptCut);
 387     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0 .9
 388     // Define reader and set its header
 389     er = new AliJetKineReader();
 390     er->SetReaderHeader(jrh);
 391     break;
 392   case "ESD":
 393     AliJetESDReaderHeader *jrh = new AliJetESDReaderHeader();
 394     jrh->SetComment("Testing");
 395     jrh->SetFirstEvent(0.);
 396     jrh->SetLastEvent(1000);
 397     jrh->SetPtCut(ptCut);
 398     jrh->SetReadSignalOnly(kFALSE);
 399     // Define reader and set its header
 400     er = new AliJetESDReader();
 401     er->SetReaderHeader(jrh);
 402     break;
 403
 404   case "AOD":
 405     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 406     jrh->SetComment("AOD Reader");
 407     jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
 408     jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
 409     if(filterMask>0)jrh->SetTestFilterMask(filterMask);
 410     // Define reader and set its header
 411     er = new AliJetAODReader();
 412     er->SetReaderHeader(jrh);
 413     break;
 414   case "AODMC":
 415     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 416     jrh->SetComment("AOD MC Reader");
 417     jrh->SetPtCut(ptCut);
 418     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
 419     jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
 420     // Define reader and set its header
 421     er = new AliJetAODReader();
 422     er->SetReaderHeader(jrh);
 423     break;
 424   case "AODMC2":
 425     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 426     jrh->SetComment("AOD MC Reader");
 427     jrh->SetPtCut(ptCut);
 428     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
 429     jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
 430     // Define reader and set its header
 431     er = new AliJetAODReader();
 432     er->SetReaderHeader(jrh);
 433     break;
 434   case "AODMC2b":
 435     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 436     jrh->SetComment("AOD MC Reader");
 437     jrh->SetPtCut(ptCut);
 438     jrh->SetFiducialEta(-0.9,0.9); // to take all MC particles default is 0.9
 439     jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
 440     // Define reader and set its header
 441     er = new AliJetAODReader();
 442     er->SetReaderHeader(jrh);
 443     break;
 444   case "AODextra":
 445     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 446     jrh->SetComment("AOD Reader with extra branch");
 447     jrh->SetPtCut(ptCut); // set low p_T cut of to ptCut
 448     jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
 449     jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
 450     jrh->SetNonStdBranch("aodExtraTracks");
 451     if(filterMask>0)jrh->SetTestFilterMask(filterMask);
 452     // Define reader and set its header
 453     er = new AliJetAODReader();
 454     er->SetReaderHeader(jrh);
 455     break;
 456   case "AODextraonly":
 457     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 458     jrh->SetComment("AOD Reader with extra branch");
 459     jrh->SetPtCut(ptCut); // set low p_T cut of to 150 MeV
 460     jrh->SetTestFilterMask(32); // Change this one for a different set of cuts
 461     jrh->SetReadBranches(AliJetAODReaderHeader::kReadNonStdBranch);
 462     jrh->SetNonStdBranch("aodExtraTracks");
 463     if(filterMask>0)jrh->SetTestFilterMask(filterMask);
 464     // Define reader and set its header
 465     er = new AliJetAODReader();
 466     er->SetReaderHeader(jrh);
 467     break;
 468   case "AODMCextra":
 469     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 470     jrh->SetComment("AOD MC Reader with extra branch");
 471     jrh->SetPtCut(ptCut);
 472     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
 473     jrh->SetReadAODMC(1);// 1 all primary MC , 2 all primary charged
 474     jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
 475     jrh->SetNonStdBranch("aodExtraMCparticles");
 476     // Define reader and set its header
 477     er = new AliJetAODReader();
 478     er->SetReaderHeader(jrh);
 479     break;
 480   case "AODMC2extra":
 481     AliJetAODReaderHeader *jrh = new AliJetAODReaderHeader();
 482     jrh->SetComment("AOD MC Reader with extra branch");
 483     jrh->SetPtCut(ptCut);
 484     jrh->SetFiducialEta(-2.1,2.1); // to take all MC particles default is 0.9
 485     jrh->SetReadAODMC(2);// 1 all primary MC , 2 all primary charged
 486     jrh->SetReadBranches(AliJetAODReaderHeader::kReadStdBranch | AliJetAODReaderHeader::kReadNonStdBranch);
 487     jrh->SetNonStdBranch("aodExtraMCparticles");
 488     // Define reader and set its header
 489     er = new AliJetAODReader();
 490     er->SetReaderHeader(jrh);
 491     break;
 492   default:
 493     ::Error("AddTaskJets", "Wrong jet reader selected\n");
 494     return 0;
 495   }
 496
 497   return er;
 498
 499 }