PWGLF/RESONANCES/macros/mini/AddAnalysisTaskD0.C

   1 /***************************************************************************
   2               fbellini@cern.ch - last modified on 06/08/2012
   3 //
   4 // General macro to configure the RSN analysis task.
   5 // It calls all configs desired by the user, by means
   6 // of the boolean switches defined in the first lines.
   7 // ---
   8 // Inputs:
   9 //  1) flag to know if running on MC or data
  10 //  2) path where all configs are stored
  11 // ---
  12 // Returns:
  13 //  kTRUE  --> initialization successful
  14 //  kFALSE --> initialization failed (some config gave errors)
  15 //
  16 ****************************************************************************/
  17
  18 AliRsnMiniAnalysisTask * AddAnalysisTaskD0
  19 (
  20    Bool_t      isMC,
  21    Bool_t      isPP,
  22    Bool_t      ispPb,
  23    Bool_t      monitor = kTRUE,
  24    Bool_t      centortracklets = kTRUE,
  25    Bool_t      sanityhistos = kTRUE,
  26    TString     centrality = "V0M",
  27    Int_t       aodFilterBit = 5,
  28    Float_t     nsigmaTPCPi = 3.0,
  29    Float_t     nsigmaTPCKa = 3.0,
  30    Float_t     nsigmaTOFPi = 2.0,
  31    Float_t     nsigmaTOFKa = 2.0,
  32    Float_t     trackDCAcutMax = 7.0,
  33    Float_t     trackDCAZcutMax = 2.0,
  34    Int_t       NTPCcluster = 70,
  35    Double_t    NTPCcrratio = 0.8,
  36    Int_t       minSPDclt = 0,
  37    Double_t    minpt = 0.15,
  38    TString     triggerMask = AliVEvent::kMB,
  39    Bool_t      useNTPCclt = kTRUE,
  40    Bool_t      maxDCAcutFixed = kFALSE,
  41    Bool_t      ptdepPIDcut = kFALSE,
  42    Bool_t      fixedYcut = kTRUE,
  43    Bool_t      checkpileup = kFALSE,
  44    Bool_t      SPDpileup = kTRUE,
  45    Bool_t      doCalculationInMC = kTRUE,
  46    UShort_t    originDselection = 0,
  47    Int_t       nmix = 5,
  48    Double_t    minYlab = -0.5,
  49    Double_t    maxYlab = 0.5,
  50    Float_t     mineta = -0.8,
  51    Float_t     maxeta = 0.8,
  52    Float_t     min_inv_mass = 0.6,
  53    Float_t     max_inv_mass = 2.2,
  54    Int_t       bins = 320,
  55    Float_t     maxDiffVzMix = 1.0,
  56    Float_t     maxDiffMultMix = 10.0,
  57    Float_t     maxDiffAngleMixDeg = 20.0,
  58    TString     outNameSuffix = "D0"
  59 )
  60 {
  61   //
  62   // -- INITIALIZATION ----------------------------------------------------------------------------
  63   // retrieve analysis manager
  64   //
  65   Float_t     trackDCAcutMin = 0.0;
  66   Bool_t      minDCAcutFixed = kTRUE;
  67   Double_t    dcaProduct = 100.0;
  68   Float_t     cutV = 10.0;
  69   Short_t     maxSisters = 2;
  70   Bool_t      checkP = kTRUE;
  71   Bool_t      checkFeedDown = kTRUE;
  72   Bool_t      checkQuark = kTRUE;
  73   Int_t       aodN = 0;
  74
  75
  76   UInt_t trigger = 0;
  77
  78
  79
  80   TString eventType = "";
  81    if(triggerMask=="AliVEvent::kMB | AliVEvent::kCentral") {trigger = AliVEvent::kMB | AliVEvent::kCentral; eventType+="Central";}
  82    else if(triggerMask=="AliVEvent::kMB | AliVEvent::kSemiCentral") {trigger = AliVEvent::kMB | AliVEvent::kSemiCentral; eventType+="SemiCentral";}
  83    else if(triggerMask=="AliVEvent::kINT7") {trigger = AliVEvent::kINT7; eventType+="kINT7";}
  84    else if(triggerMask=="AliVEvent::kMB") {trigger = AliVEvent::kMB; eventType+="MinimumBias";}
  85
  86   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  87    if (!mgr) {
  88       ::Error("AddAnalysisTaskD0", "No analysis manager to connect to.");
  89       return NULL;
  90    }
  91
  92    // create the task and configure
  93    TString taskName = Form("D0%s%s_%.1f_%d_%.2f_%d_%.1f_%.1f_%.1f_%.1f_%.1f_%.1f_%.2f_%d_%s", (isPP? "pp" : ispPb? "pPB": "PbPb"), (isMC ? "MC" : "Data"), cutV, NTPCcluster, NTPCcrratio, minSPDclt, nsigmaTPCPi, nsigmaTPCKa, nsigmaTOFPi, nsigmaTOFKa, trackDCAcutMax, trackDCAZcutMax, minpt, originDselection, eventType.Data());
  94    AliRsnMiniAnalysisTask *task = new AliRsnMiniAnalysisTask(taskName.Data(), isMC);
  95    if (!isMC && !isPP){
  96      Printf(Form("========== SETTING USE CENTRALITY PATCH AOD049 : %s", (aodN==49)? "yes" : "no"));
  97      task->SetUseCentralityPatch(aodN==49);
  98    }
  99
 100    ::Info("AddAnalysisTaskD0", Form("TriggerMask: %i",trigger));
 101
 102    task->SelectCollisionCandidates(trigger);
 103    task->SetMaxNDaughters(maxSisters);
 104    task->SetCheckMomentumConservation(checkP);
 105    task->SetCheckFeedDown(checkFeedDown);
 106    task->SetRejectCandidateIfNotFromQuark(checkQuark);
 107    task->SetDselection(originDselection);
 108    task->KeepMotherInAcceptance(kTRUE);
 109    task->SetMotherAcceptanceCutMinPt(minpt);
 110    task->SetMotherAcceptanceCutMaxEta(maxeta);
 111
 112
 113    ::Info("AddAnalysisTaskD0", Form("Maximum numbers of daughters allowed (-1 means cut not applied): %i",maxSisters));
 114    ::Info("AddAnalysisTaskD0", Form("Are we checking the momentum conservation? %s", checkP? "yes" : "no"));
 115    ::Info("AddAnalysisTaskD0", Form("Are we checking the feeddown? %s", checkFeedDown? "yes" : "no"));
 116    ::Info("AddAnalysisTaskD0", Form("Are we rejecting the Hijing generated? %s", checkQuark? "yes" : "no"));
 117    ::Info("AddAnalysisTaskD0", Form("Which D0 are we keeping? %s", (originDselection==0? "only from c quark" : originDselection==1? "only from b quark" : "both from c and b quark") ));
 118    ::Info("AddAnalysisTaskD0", Form("Selecting Mother in Acceptance: Min pT %.1f, Eta Range %.1f - %.1f", minpt, mineta, maxeta));
 119
 120
 121    if (isPP)
 122      task->UseMultiplicity("QUALITY");
 123    else
 124      task->UseCentrality(centrality);
 125    // set event mixing options
 126    task->UseContinuousMix();
 127    //task->UseBinnedMix();
 128    task->SetNMix(nmix);
 129    task->SetMaxDiffVz(maxDiffVzMix);
 130    task->SetMaxDiffMult(maxDiffMultMix);
 131    if (!isPP && !ispPb) task->SetMaxDiffAngle(maxDiffAngleMixDeg*TMath::DegToRad()); //set angle diff in rad
 132    ::Info("AddAnalysisTaskD0", Form("Event mixing configuration: \n events to mix = %i \n max diff. vtxZ = cm %5.3f \n max diff multi = %5.3f \n max diff EP angle = %5.3f deg", nmix, maxDiffVzMix, maxDiffMultMix, (isPP ? 0.0 : ispPb ? 0.0 : maxDiffAngleMixDeg)));
 133
 134    mgr->AddTask(task);
 135
 136    //
 137    // -- EVENT CUTS (same for all configs) ---------------------------------------------------------
 138    //
 139    // cut on primary vertex:
 140    // - 2nd argument --> |Vz| range
 141    // - 3rd argument --> minimum required number of contributors
 142    // - 4th argument --> tells if TPC stand-alone vertexes must be accepted
 143    AliRsnCutPrimaryVertex *cutVertex = new AliRsnCutPrimaryVertex("cutVertex", cutV, 1, kFALSE, kFALSE);
 144    AliRsnCutEventUtils *cutEventUtils = new AliRsnCutEventUtils("cutEventUtils", kFALSE, kTRUE);
 145
 146    if(checkpileup == kTRUE){
 147
 148         if(SPDpileup == kTRUE) {cutEventUtils->SetRemovePileUppA2013(kTRUE);
 149                                 //cutEventUtils->SetRemoveFirstEvtInChunk(kTRUE);
 150                                 //cutEventUtils->SetUseVertexSelection2013pA(kTRUE);
 151                                 cutEventUtils->SetUseMVPlpSelection(kFALSE);
 152                                 cutEventUtils->SetMinPlpContribSPD(5);
 153                                                                           }
 154         if(SPDpileup == kFALSE){cutEventUtils->SetRemovePileUppA2013(kTRUE);
 155                                 //cutEventUtils->SetRemoveFirstEvtInChunk(kTRUE);
 156                                 //cutEventUtils->SetUseVertexSelection2013pA(kTRUE);
 157                                 cutEventUtils->SetUseMVPlpSelection(kTRUE);
 158                                 cutEventUtils->SetMinPlpContribMV(5);
 159                                                                           }
 160    }
 161
 162    ::Info("AddAnalysisTaskD0", Form("Checking Pile up? %s", (checkpileup ? "yes" : "no") ));
 163    ::Info("AddAnalysisTaskD0", Form("Which Method? %s", (checkpileup? (SPDpileup ? "SPD" : "Multi Vertex") : "none" )));
 164
 165    // define and fill cut set for event cut
 166    AliRsnCutSet *eventCuts = new AliRsnCutSet("eventCuts", AliRsnTarget::kEvent);
 167    eventCuts->AddCut(cutVertex);
 168    //if(SPDpileup == kFALSE) eventCuts->AddCut(cutEventUtils);
 169    eventCuts->AddCut(cutEventUtils);
 170    eventCuts->SetCutScheme(Form("%s&%s", cutEventUtils->GetName(), cutVertex->GetName() ));
 171    eventCuts->ShowCuts();
 172    eventCuts->PrintSetInfo();
 173    // set cuts in task
 174    task->SetEventCuts(eventCuts);
 175
 176    //
 177    // -- EVENT-ONLY COMPUTATIONS -------------------------------------------------------------------
 178    //
 179    //vertex
 180    Int_t vtxID = task->CreateValue(AliRsnMiniValue::kVz, kFALSE);
 181    AliRsnMiniOutput *outVtx = task->CreateOutput("eventVtx", "HIST", "EVENT");
 182    outVtx->AddAxis(vtxID, 400, -20.0, 20.0);
 183
 184    //multiplicity or centrality
 185    Int_t multID = task->CreateValue(AliRsnMiniValue::kMult, kFALSE);
 186    AliRsnMiniOutput *outMult = task->CreateOutput("eventMult", "HIST", "EVENT");
 187    if (isPP)
 188      outMult->AddAxis(multID, 400, 0.0, 400.0);
 189    else
 190      outMult->AddAxis(multID, 100, 0.0, 100.0);
 191
 192    //tracklets
 193    Int_t trackletID = task->CreateValue(AliRsnMiniValue::kTracklets, kFALSE);
 194    AliRsnMiniOutput *outTracklets = task->CreateOutput("eventTracklets", "HIST", "EVENT");
 195    outTracklets->AddAxis(trackletID, 250, -0.5, 249.5);
 196
 197
 198    //event plane (only for PbPb)
 199    Int_t planeID = task->CreateValue(AliRsnMiniValue::kPlaneAngle, kFALSE);
 200    AliRsnMiniOutput *outPlane = 0x0;
 201    if (!isPP){
 202      outPlane = task->CreateOutput("eventPlane", "HIST", "EVENT");
 203      outPlane->AddAxis(planeID, 180, 0.0, TMath::Pi());
 204      }
 205
 206    TH2F* hvz=new TH2F("hVzVsCent","",100,0.,100., 220,-11.,11.);
 207    task->SetEventQAHist("vz",hvz);//plugs this histogram into the fHAEventVz data member
 208
 209    TH2F* hmc=new TH2F("MultiVsCent","",100,0.,100.,1000,0.,1000.);
 210    hmc->GetYaxis()->SetTitle("QUALITY");
 211    task->SetEventQAHist("multicent",hmc);//plugs this histogram into the fHAEventMultiCent data member
 212
 213    TH2F* hep=new TH2F("hEventPlaneVsCent","",100,0.,100., 180,0.,TMath::Pi());
 214    task->SetEventQAHist("eventplane",hep);//plugs this histogram into the fHAEventPlane data member
 215
 216    //
 217    // -- PAIR CUTS (common to all resonances) ------------------------------------------------------
 218    //
 219    AliRsnCutMiniPair *cutY = 0x0;
 220
 221    if(fixedYcut == kTRUE){
 222    cutY = new AliRsnCutMiniPair("cutRapidity", AliRsnCutMiniPair::kRapidityRange);
 223    cutY->SetRangeD(minYlab, maxYlab);
 224    }
 225    else {
 226    cutY = new AliRsnCutMiniPair("cutRapidity", AliRsnCutMiniPair::kRapidityFiducialRegion);
 227    cutY->SetRangeD(-0.8, 0.8);
 228    cutY->SetPtDepCut(kTRUE);
 229    cutY->SetMinPt(0.0);
 230    cutY->SetMaxPt(5.0);
 231    cutY->SetPtDepCutMaxFormula("-0.2/15*pt*pt+1.9/15*pt+0.5");
 232    cutY->SetPtDepCutMinFormula("0.2/15*pt*pt-1.9/15*pt-0.5");
 233    }
 234
 235    ::Info("AddAnalysisTaskD0", Form("Fixed Y cut? %s", (fixedYcut ? "yes" : "no") ));
 236
 237    AliRsnCutMiniPair *cutDCAproduct = new AliRsnCutMiniPair("cutDCAproduct", AliRsnCutMiniPair::kDCAproduct);
 238    cutDCAproduct->SetRangeD(-1E20, dcaProduct);
 239
 240    AliRsnCutSet *cutsPairY = new AliRsnCutSet("pairCutsY", AliRsnTarget::kMother);
 241    cutsPairY->AddCut(cutY);
 242    cutsPairY->UseMonitor(kTRUE);
 243    cutsPairY->SetCutScheme("setPairD0_Y");
 244    cutsPairY->ShowCuts();
 245    cutsPairY->PrintSetInfo();
 246
 247    AliRsnCutSet *cutsPair = new AliRsnCutSet("pairCuts", AliRsnTarget::kMother);
 248    cutsPair->AddCut(cutY);
 249    cutsPair->AddCut(cutDCAproduct);
 250    cutsPair->UseMonitor(kTRUE);
 251    cutsPair->SetCutScheme(Form("%s&%s", cutY->GetName(), cutDCAproduct->GetName()));
 252    cutsPair->ShowCuts();
 253    cutsPair->PrintSetInfo();
 254
 255    //
 256    // -- CONFIG ANALYSIS --------------------------------------------------------------------------
 257    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/RESONANCES/macros/mini/ConfigD0.C");
 258
 259    if (isMC) {
 260        Printf("========================== MC analysis - PID cuts used");
 261    } else
 262      Printf("========================== DATA analysis - PID cuts used");
 263    if (!ConfigD0(task, isPP, isMC, monitor, centortracklets, sanityhistos, nsigmaTPCPi, nsigmaTPCKa, nsigmaTOFPi, nsigmaTOFKa, aodFilterBit, trackDCAcutMax, trackDCAcutMin, trackDCAZcutMax, NTPCcluster, NTPCcrratio, minSPDclt, minpt, maxSisters, checkP, useNTPCclt, minDCAcutFixed, maxDCAcutFixed, ptdepPIDcut, checkFeedDown, checkQuark, doCalculationInMC, originDselection, mineta, maxeta, min_inv_mass, max_inv_mass, bins, "", cutsPairY, cutsPair)) return 0x0;
 264
 265    //
 266    // -- CONTAINERS --------------------------------------------------------------------------------
 267    //
 268    TString outputFileName = AliAnalysisManager::GetCommonFileName();
 269    Printf("AddAnalysisTaskD0 - Set OutputFileName : \n %s\n", outputFileName.Data() );
 270
 271    AliAnalysisDataContainer *output = mgr->CreateContainer(Form("%s_%.1f_%d_%.2f_%d_%.1f_%.1f_%.1f_%.1f_%.1f_%.1f_%.2f_%d_%s",outNameSuffix.Data(),cutV,NTPCcluster,NTPCcrratio,minSPDclt,nsigmaTPCPi,nsigmaTPCKa,nsigmaTOFPi,nsigmaTOFKa,trackDCAcutMax,trackDCAZcutMax,minpt,originDselection,eventType.Data()),
 272                                                            TList::Class(),
 273                                                            AliAnalysisManager::kOutputContainer,
 274                                                            outputFileName);
 275    mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
 276    mgr->ConnectOutput(task, 1, output);
 277
 278    return task;
 279 }