b67dc968f132cd4a64dab5340e5855ca0ce98d77
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / AddTaskPIDFlowSP.C
1 //AddTask macro for higher harmonic flow analysis with SP method only.
2
3 class AliAnalysisDataContainer;
4 class AliFlowTrackCuts;
5 class AliFlowEventCuts;
6
7
8 void AddTaskFlowSP(Int_t triggerSelectionString=AliVEvent::kMB,
9                                    Float_t etamin=-0.8,
10                                    Float_t etamax=0.8,
11                                    Float_t EtaGap=0.2,
12                                    TString fileNameBase="AnalysisResults",
13                                    TString uniqueStr="Pion_02",
14                                    TString Qvector ="Qa",
15                                    Int_t AODfilterBitRP = 272,
16                                    Int_t AODfilterBitPOI = 272,
17                                    Int_t charge=0,
18                                    Int_t MinTPCdedx = 10,
19                                    Int_t ncentrality = 6,
20                                    Bool_t doQA=kTRUE,
21                                    Bool_t isPID = kTRUE,
22                                    AliPID::EParticleType particleType=AliPID::kPion,
23                                    AliFlowTrackCuts::PIDsource sourcePID=AliFlowTrackCuts::kTOFbayesian) {
24
25 // Define a range of the detector to exclude
26 Bool_t ExcludeRegion = kFALSE;
27 Double_t excludeEtaMin = -0.;
28 Double_t excludeEtaMax = 0.;
29 Double_t excludePhiMin = 0.;
30 Double_t excludePhiMax = 0.;
31    
32 //Define the range for eta subevents (for SP method)
33 Double_t minA = -0.8;//
34 Double_t maxA = -0.5*EtaGap;//
35 Double_t minB = +0.5*EtaGap;//
36 Double_t maxB = +0.8;//
37
38     
39 int centrMin[9] = {0,5,10,20,30,40,50,60,70};
40 int centrMax[9] = {5,10,20,30,40,50,60,70,80};
41 const int ncentr = ncentrality;
42     
43     
44 //---------Data selection----------
45 //kMC, kGlobal, kESD_TPConly, kESD_SPDtracklet
46 AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kTPCstandalone;
47 AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kTPCstandalone;
48     
49 //---------Parameter mixing--------
50 //kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
51 AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
52 AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
53     
54 const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype);
55 const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype);
56
57     
58 //===========================================================================
59 // EVENTS CUTS:
60 AliFlowEventCuts* cutsEvent[ncentr];
61 AliFlowTrackCuts* cutsRP[ncentr];
62 AliFlowTrackCuts* cutsPOI[ncentr];
63 TString outputSlotName[ncentr][4];
64 TString suffixName[ncentr];
65     
66 for(int icentr=0;icentr<ncentr;icentr++){
67     cutsEvent[icentr] = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
68     //cutsEvent[icentr]->SetUsedDataset(is2011);
69     cutsEvent[icentr]->SetCentralityPercentileRange(centrMin[icentr],centrMax[icentr]);
70     cutsEvent[icentr]->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
71     //  cutsEvent->SetRefMultMethod(AliFlowEventCuts::kVZERO);
72     //cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
73     //cutsEvent->SetNContributorsRange(2);
74     cutsEvent[icentr]->SetPrimaryVertexZrange(-10.,10.);
75     cutsEvent[icentr]->SetQA(doQA);
76     cutsEvent[icentr]->SetCutTPCmultiplicityOutliers();
77     
78     
79     // RP TRACK CUTS:
80     cutsRP[icentr] = new AliFlowTrackCuts(Form("TPConlyRP_%d",icentr));
81     cutsRP[icentr]->SetParamType(rptype);
82     cutsRP[icentr]->SetParamMix(rpmix);
83     cutsRP[icentr]->SetPtRange(0.2,5.);
84     cutsRP[icentr]->SetEtaRange(etamin,etamax);
85     cutsRP[icentr]->SetMinNClustersTPC(70);
86     //  cutsRP->SetMinChi2PerClusterTPC(0.1);//
87     // cutsRP->SetMaxChi2PerClusterTPC(4.0);//
88     cutsRP[icentr]->SetMaxDCAToVertexXY(3.0);
89     cutsRP[icentr]->SetMaxDCAToVertexZ(3.0);
90     cutsRP[icentr]->SetAcceptKinkDaughters(kFALSE);
91     cutsRP[icentr]->SetMinimalTPCdedx(MinTPCdedx);
92     cutsRP[icentr]->SetAODfilterBit(AODfilterBitRP);
93     cutsRP[icentr]->SetQA(doQA);
94
95     
96     //POIs for SP and QC method
97     //===========================================================================
98     AliFlowTrackCuts  *SP_POI[ncentr];
99     //half window for POIs
100     //=======================SP POI Cuts
101     SP_POI[icentr] = DefinePOIcuts(icentr);
102
103     SP_POI[icentr]->GetBayesianResponse()->ForceOldDedx(); // for 2010 data to use old TPC PID Response instead of the official one
104     SP_POI[icentr]->SetParamType(poitype);
105     SP_POI[icentr]->SetParamMix(poimix);
106     SP_POI[icentr]->SetPtRange(0.2,5.);//
107     SP_POI[icentr]->SetMinNClustersTPC(70);
108     if(Qvector=="Qa"){
109         SP_POI[icentr]->SetEtaRange( +0.5*EtaGap, etamax );
110     }
111     if(Qvector=="Qb"){
112         SP_POI[icentr]->SetEtaRange( etamin,-0.5*EtaGap );
113     }
114
115     // SP_POI->SetMinChi2PerClusterTPC(0.1); //
116     // SP_POI->SetMaxChi2PerClusterTPC(4.0); //
117     //  SP_POI->SetRequireITSRefit(kTRUE);
118     //  SP_POI->SetRequireTPCRefit(kTRUE);
119     //  SP_POI->SetMinNClustersITS(2);
120     //  SP_POI->SetMaxChi2PerClusterITS(1.e+09);
121     SP_POI[icentr]->SetMaxDCAToVertexXY(3.0);
122     SP_POI[icentr]->SetMaxDCAToVertexZ(3.0);
123     //SP_POI->SetDCAToVertex2D(kTRUE);
124     //SP_POI->SetMaxNsigmaToVertex(1.e+10);
125     //SP_POI->SetRequireSigmaToVertex(kFALSE);
126     SP_POI[icentr]->SetAcceptKinkDaughters(kFALSE);
127     if(isPID) SP_POI[icentr]->SetPID(particleType, sourcePID);//particleType, sourcePID
128     if (charge!=0) SP_POI[icentr]->SetCharge(charge);
129     //SP_POI->SetAllowTOFmismatch(kFALSE);
130     SP_POI[icentr]->SetRequireStrictTOFTPCagreement(kTRUE);
131     SP_POI[icentr]->SetMinimalTPCdedx(MinTPCdedx);
132     SP_POI[icentr]->SetAODfilterBit(AODfilterBitPOI);
133     SP_POI[icentr]->SetQA(doQA);
134     SP_POI[icentr]->SetPriors((centrMin[icentr]+centrMax[icentr])*0.5);
135
136
137
138
139     //=====================================================================
140  
141     if(Qvector=="Qa") suffixName[icentr] = "Qa";
142     if(Qvector=="Qb") suffixName[icentr] = "Qb";
143     suffixName[icentr] += "-highharmflow";
144     suffixName[icentr] += Form("%i_", centrMin[icentr]);
145     suffixName[icentr] += Form("%i_", centrMax[icentr]);
146     suffixName[icentr] += Form("%.f_", EtaGap*10);
147
148     if(isPID){
149         suffixName[icentr]+=AliFlowTrackCuts::PIDsourceName(sourcePID);
150         suffixName[icentr]+="_";
151         suffixName[icentr]+=AliPID::ParticleName(particleType);//particleType
152     }
153     else{
154         suffixName[icentr]+="AllCharged";
155     }
156     if (charge<0) suffixName[icentr]+="-";
157     if (charge>0) suffixName[icentr]+="+";
158
159
160     for(int harmonic=2;harmonic<6;harmonic++){  //for v2,v3,v4 and v5
161         outputSlotName[icentr][harmonic-2] = "";
162         outputSlotName[icentr][harmonic-2]+=uniqueStr;
163         outputSlotName[icentr][harmonic-2]+=Form("_v%i_",harmonic);
164         outputSlotName[icentr][harmonic-2]+=cutsRP[icentr]->GetName();
165         outputSlotName[icentr][harmonic-2]+="_";
166         outputSlotName[icentr][harmonic-2]+=SP_POI[icentr]->GetName();
167         outputSlotName[icentr][harmonic-2]+=Form("_%i-",centrMin[icentr]);
168         outputSlotName[icentr][harmonic-2]+=Form("%i_",centrMax[icentr]);
169         
170         
171         if(isPID){
172             outputSlotName[icentr][harmonic-2]+=AliFlowTrackCuts::PIDsourceName(sourcePID);//sourcePID
173             outputSlotName[icentr][harmonic-2]+="_";
174             outputSlotName[icentr][harmonic-2]+=AliPID::ParticleName(particleType);//particleType
175         }
176         else{
177             outputSlotName[icentr][harmonic-2]+="AllCharged";
178         }
179         if (charge<0) outputSlotName[icentr][harmonic-2]+="-";
180         if (charge>0) outputSlotName[icentr][harmonic-2]+="+";
181     }
182 }
183
184
185 TString fileName(fileNameBase);
186 fileName.Append(".root");
187
188
189     
190 //====================================FLOWPACKAGE TASKS=========================//
191 AliAnalysisDataContainer *cinput1[ncentr];
192 AliAnalysisDataContainer *coutputFE[ncentr];
193 AliAnalysisDataContainer* coutputFEQA[ncentr];
194 AliAnalysisTaskFlowEvent *taskFE[ncentr];
195
196 AliAnalysisDataContainer *flowEvent[ncentr][4];
197 AliAnalysisTaskFilterFE *tskFilter[ncentr][4];
198
199 AliAnalysisDataContainer *coutputSP[ncentr][4];
200 AliAnalysisTaskScalarProduct *taskSP[ncentr][4];
201
202 TString outputQA[ncentr];
203 TString myNameSP[ncentr][4];
204 TString slot[ncentr][4];
205
206 for (int icentr=0; icentr<ncentr; icentr++) {
207         
208     // Get the pointer to the existing analysis manager via the static access method.
209         //==============================================================================
210     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
211     if (!mgr) {
212         Error("AddTaskFlowEvent", "No analysis manager to connect to.");
213         return NULL;
214     }
215         
216     // Check the analysis type using the event handlers connected to the analysis
217     // manager. The availability of MC handler can also be checked here.
218     //==============================================================================
219     if (!mgr->GetInputEventHandler()) {
220         ::Error("AddTaskFlowEvent", "This task requires an input event handler");
221         return NULL;
222     }
223
224     taskFE[icentr] = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent_%s",suffixName[icentr].Data()),"",doQA);
225     taskFE[icentr]->SelectCollisionCandidates(triggerSelectionString);
226     taskFE[icentr]->SetSubeventEtaRange(minA, maxA, minB, maxB);
227     mgr->AddTask(taskFE[icentr]);
228
229     // Pass cuts for RPs and POIs to the task:
230     taskFE[icentr]->SetCutsEvent(cutsEvent[icentr]);
231     taskFE[icentr]->SetCutsRP(cutsRP[icentr]);
232     taskFE[icentr]->SetCutsPOI(SP_POI[icentr]);
233     if (cutsRP[icentr]->GetParamType()==AliFlowTrackCuts::kVZERO)
234     {
235         //TODO: since this is set in a static object all analyses in an analysis train
236         //will be affected.
237         taskFE[icentr]->SetHistWeightvsPhiMin(0.);
238         taskFE[icentr]->SetHistWeightvsPhiMax(200.);
239     }
240     
241     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
242     
243     cinput1[icentr] = mgr->GetCommonInputContainer();
244     
245     coutputFE[icentr] = mgr->CreateContainer(Form("FlowEvent_%s",suffixName[icentr].Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
246
247     mgr->ConnectInput(taskFE[icentr],0,cinput1[icentr]);
248     mgr->ConnectOutput(taskFE[icentr],1,coutputFE[icentr]);
249     //==========================================================
250     
251     TString Species = "";
252     if(isPID) Species += AliPID::ParticleName(particleType);
253     else      Species += "Allcharged";
254         
255
256     for(int harm=2;harm<6;harm++){
257         myNameSP[icentr][harm-2] = "SP_";
258         myNameSP[icentr][harm-2] += Qvector;
259         myNameSP[icentr][harm-2] += Form("_v%i_%s_%.f",harm,outputSlotName[icentr][harm-2].Data(),EtaGap*10);
260         
261         flowEvent[icentr][harm-2] = mgr->CreateContainer( Form("Filter_%s", myNameSP[icentr][harm-2].Data()),
262                                                                    AliFlowEventSimple::Class(),
263                                                                    AliAnalysisManager::kExchangeContainer );
264         
265         tskFilter[icentr][harm-2] = new AliAnalysisTaskFilterFE( Form("TaskFilter_%s",myNameSP[icentr][harm-2].Data()),cutsRP[icentr], NULL);
266         tskFilter[icentr][harm-2]->SetSubeventEtaRange(etamin, -.5*EtaGap, +.5*EtaGap, etamax);
267         
268         mgr->AddTask(tskFilter[icentr][harm-2]);
269         mgr->ConnectInput( tskFilter[icentr][harm-2],0,coutputFE[icentr]);
270         mgr->ConnectOutput(tskFilter[icentr][harm-2],1,flowEvent[icentr][harm-2]);
271         
272         
273         taskSP[icentr][harm-2] = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s",outputSlotName[icentr][harm-2].Data()),kFALSE);
274         taskSP[icentr][harm-2]->SetHarmonic(harm);
275         taskSP[icentr][harm-2]->SelectCollisionCandidates(triggerSelectionString);
276         taskSP[icentr][harm-2]->SetRelDiffMsub(1.0);
277         taskSP[icentr][harm-2]->SetTotalQvector(Qvector);
278         taskSP[icentr][harm-2]->SetApplyCorrectionForNUA(kTRUE);
279         
280         TString outputSP = fileName;
281         outputSP += ":outputSPanalysis";
282         outputSP+= rptypestr;
283         slot[icentr][harm-2] = "SP_";
284         slot[icentr][harm-2] += outputSlotName[icentr][harm-2];
285         slot[icentr][harm-2] += "_";
286         slot[icentr][harm-2] += Qvector;
287         coutputSP[icentr][harm-2] = mgr->CreateContainer(Form("%s_%.f",slot[icentr][harm-2].Data(),EtaGap*10),
288                                          TList::Class(),AliAnalysisManager::kOutputContainer,outputSP);
289         mgr->AddTask(taskSP[icentr][harm-2]);
290         mgr->ConnectInput(taskSP[icentr][harm-2],0,flowEvent[icentr][harm-2]);
291         mgr->ConnectOutput(taskSP[icentr][harm-2],1,coutputSP[icentr][harm-2]);
292     }
293
294
295     if (taskFE[icentr]->GetQAOn()) {
296         outputQA[icentr] = fileName;
297         outputQA[icentr] += ":QA";
298         coutputFEQA[icentr] = mgr->CreateContainer(Form("QA_%s",suffixName[icentr].Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputQA[icentr]);
299         mgr->ConnectOutput(taskFE[icentr],2,coutputFEQA[icentr]);
300     }
301
302
303 }
304 }
305
306 //===========================================================================
307
308 AliFlowEventCuts* DefinecutsEvent(Int_t icentr){
309     AliFlowEventCuts* cutsEvent = new AliFlowEventCuts(Form("eventcuts_%d",icentr));
310     return cutsEvent;
311 }
312 AliFlowTrackCuts* DefineRPcuts(Int_t icentr){
313     AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts(Form("TPConlyRP_%d",icentr));
314     return cutsRP;
315 }
316 AliFlowTrackCuts* DefinePOIcuts(Int_t icentr){
317     AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts(Form("TPConlyPOI_%d",icentr));
318     return cutsPOI;
319 }
320
321
322
323