Update by Chitrasen - Modified histrogram axis
[u/mrichter/AliRoot.git] / 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 }