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