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    if(checkpileup == kTRUE){
 145         AliRsnCutEventUtils *cutEventUtils = new AliRsnCutEventUtils("cutEventUtils", kFALSE, kTRUE);
 146         if(SPDpileup == kTRUE) {cutEventUtils->SetRemovePileUppA2013(kTRUE);
 147                                 //cutEventUtils->SetRemoveFirstEvtInChunk(kTRUE);
 148                                 //cutEventUtils->SetUseVertexSelection2013pA(kTRUE);
 149                                 cutEventUtils->SetUseMVPlpSelection(kFALSE);
 150                                 cutEventUtils->SetMinPlpContribSPD(5);
 151                                                                           }
 152         if(SPDpileup == kFALSE){cutEventUtils->SetRemovePileUppA2013(kTRUE);
 153                                 //cutEventUtils->SetRemoveFirstEvtInChunk(kTRUE);
 154                                 //cutEventUtils->SetUseVertexSelection2013pA(kTRUE);
 155                                 cutEventUtils->SetUseMVPlpSelection(kTRUE);
 156                                 cutEventUtils->SetMinPlpContribMV(5);
 157                                                                           }
 158    }
 159
 160    ::Info("AddAnalysisTaskD0", Form("Checking Pile up? %s", (checkpileup ? "yes" : "no") ));
 161    ::Info("AddAnalysisTaskD0", Form("Which Method? %s", (checkpileup? (SPDpileup ? "SPD" : "Multi Vertex") : "none" )));
 162
 163    // define and fill cut set for event cut
 164    AliRsnCutSet *eventCuts = new AliRsnCutSet("eventCuts", AliRsnTarget::kEvent);
 165    eventCuts->AddCut(cutVertex);
 166    //if(SPDpileup == kFALSE) eventCuts->AddCut(cutEventUtils);
 167    eventCuts->AddCut(cutEventUtils);
 168    eventCuts->SetCutScheme(Form("%s&%s", cutEventUtils->GetName(), cutVertex->GetName() ));
 169    eventCuts->ShowCuts();
 170    eventCuts->PrintSetInfo();
 171    // set cuts in task
 172    task->SetEventCuts(eventCuts);
 173
 174    //
 175    // -- EVENT-ONLY COMPUTATIONS -------------------------------------------------------------------
 176    //
 177    //vertex
 178    Int_t vtxID = task->CreateValue(AliRsnMiniValue::kVz, kFALSE);
 179    AliRsnMiniOutput *outVtx = task->CreateOutput("eventVtx", "HIST", "EVENT");
 180    outVtx->AddAxis(vtxID, 400, -20.0, 20.0);
 181
 182    //multiplicity or centrality
 183    Int_t multID = task->CreateValue(AliRsnMiniValue::kMult, kFALSE);
 184    AliRsnMiniOutput *outMult = task->CreateOutput("eventMult", "HIST", "EVENT");
 185    if (isPP)
 186      outMult->AddAxis(multID, 400, 0.0, 400.0);
 187    else
 188      outMult->AddAxis(multID, 100, 0.0, 100.0);
 189
 190    //tracklets
 191    Int_t trackletID = task->CreateValue(AliRsnMiniValue::kTracklets, kFALSE);
 192    AliRsnMiniOutput *outTracklets = task->CreateOutput("eventTracklets", "HIST", "EVENT");
 193    outTracklets->AddAxis(trackletID, 400, 0.0, 400.0);
 194
 195
 196    //event plane (only for PbPb)
 197    Int_t planeID = task->CreateValue(AliRsnMiniValue::kPlaneAngle, kFALSE);
 198    AliRsnMiniOutput *outPlane = 0x0;
 199    if (!isPP){
 200      outPlane = task->CreateOutput("eventPlane", "HIST", "EVENT");
 201      outPlane->AddAxis(planeID, 180, 0.0, TMath::Pi());
 202      }
 203
 204    TH2F* hvz=new TH2F("hVzVsCent","",100,0.,100., 220,-11.,11.);
 205    task->SetEventQAHist("vz",hvz);//plugs this histogram into the fHAEventVz data member
 206
 207    TH2F* hmc=new TH2F("MultiVsCent","",100,0.,100.,4000,0.,4000.);
 208    hmc->GetYaxis()->SetTitle("QUALITY");
 209    task->SetEventQAHist("multicent",hmc);//plugs this histogram into the fHAEventMultiCent data member
 210
 211    TH2F* hep=new TH2F("hEventPlaneVsCent","",100,0.,100., 180,0.,TMath::Pi());
 212    task->SetEventQAHist("eventplane",hep);//plugs this histogram into the fHAEventPlane data member
 213
 214    //
 215    // -- PAIR CUTS (common to all resonances) ------------------------------------------------------
 216    //
 217    AliRsnCutMiniPair *cutY = 0x0;
 218
 219    if(fixedYcut == kTRUE){
 220    cutY = new AliRsnCutMiniPair("cutRapidity", AliRsnCutMiniPair::kRapidityRange);
 221    cutY->SetRangeD(minYlab, maxYlab);
 222    }
 223    else {
 224    cutY = new AliRsnCutMiniPair("cutRapidity", AliRsnCutMiniPair::kRapidityFiducialRegion);
 225    cutY->SetRangeD(-0.8, 0.8);
 226    cutY->SetPtDepCut(kTRUE);
 227    cutY->SetMinPt(0.0);
 228    cutY->SetMaxPt(5.0);
 229    cutY->SetPtDepCutMaxFormula("-0.2/15*pt*pt+1.9/15*pt+0.5");
 230    cutY->SetPtDepCutMinFormula("0.2/15*pt*pt-1.9/15*pt-0.5");
 231    }
 232
 233    ::Info("AddAnalysisTaskD0", Form("Fixed Y cut? %s", (fixedYcut ? "yes" : "no") ));
 234
 235    AliRsnCutMiniPair *cutDCAproduct = new AliRsnCutMiniPair("cutDCAproduct", AliRsnCutMiniPair::kDCAproduct);
 236    cutDCAproduct->SetRangeD(-1E20, dcaProduct);
 237
 238    AliRsnCutSet *cutsPairY = new AliRsnCutSet("pairCutsY", AliRsnTarget::kMother);
 239    cutsPairY->AddCut(cutY);
 240    cutsPairY->UseMonitor(kTRUE);
 241    cutsPairY->SetCutScheme("setPairD0_Y");
 242    cutsPairY->ShowCuts();
 243    cutsPairY->PrintSetInfo();
 244
 245    AliRsnCutSet *cutsPair = new AliRsnCutSet("pairCuts", AliRsnTarget::kMother);
 246    cutsPair->AddCut(cutY);
 247    cutsPair->AddCut(cutDCAproduct);
 248    cutsPair->UseMonitor(kTRUE);
 249    cutsPair->SetCutScheme(Form("%s&%s", cutY->GetName(), cutDCAproduct->GetName()));
 250    cutsPair->ShowCuts();
 251    cutsPair->PrintSetInfo();
 252
 253    //
 254    // -- CONFIG ANALYSIS --------------------------------------------------------------------------
 255    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/RESONANCES/macros/mini/ConfigD0.C");
 256
 257    if (isMC) {
 258        Printf("========================== MC analysis - PID cuts used");
 259    } else
 260      Printf("========================== DATA analysis - PID cuts used");
 261    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;
 262
 263    //
 264    // -- CONTAINERS --------------------------------------------------------------------------------
 265    //
 266    TString outputFileName = AliAnalysisManager::GetCommonFileName();
 267    Printf("AddAnalysisTaskD0 - Set OutputFileName : \n %s\n", outputFileName.Data() );
 268
 269    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()),
 270                                                            TList::Class(),
 271                                                            AliAnalysisManager::kOutputContainer,
 272                                                            outputFileName);
 273    mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
 274    mgr->ConnectOutput(task, 1, output);
 275
 276    return task;
 277 }