updated analysis (Jason Ulery)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / macros / 3particlecorrelations / AddTask_DiHadron.C
1 //#include "exception.h"
2 //For running on PbPb data 0-50% most central
3 AliAnalysisTask *AddTask_DiHadron(Int_t IncludeLowPtBins=0){
4
5   //get the current analysis manager
6   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
7   if (!mgr) {
8     Error("AddTask_sma_PriVtx", "No analysis manager found.");
9     return 0;
10   }
11
12  
13   //========= SetInitial Parameters =====
14   //Int_t IncludeLowPtBins=0;//Set to 1 to include low pt triggers
15
16   //Track Quality Cuts
17   Int_t MinimumClustersTPC=70;
18   Float_t MinClusterRatio=0.51;//Must have at least this ratio not shared
19   Float_t MaxTPCchi2=4;
20   Int_t MinimumClustersITS=0;
21   Float_t EtaCut=0.9;//Tracks in +/- Eta are used
22   Float_t TriggerEtaCut=0.5;//Trigger particle restriction so flat area in acceptanc can be created
23   Float_t NearPhiCut=1.5;//Cut used to seperate near and away side for delta eta plots
24   Float_t XECut=NearPhiCut;//For XE distribution near and away seperation
25   Float_t MaxDCA=3;//Total DCA Cut
26   Float_t MaxDCAXY=2.4;
27   Float_t MaxDCAZ=3.2;
28   Int_t DCAMethod=1;//0 MaxDCA used, 1 MaxDCAXY and MaxDCAZ used 2 pT dependent DCA cut
29   Int_t TPCRefit=1;
30   Int_t ITSRefit=0;//1 for all particles, 2 for particles above 5 GeV/c
31   Int_t SPDCut=0;//check for a point in 1 of first 2 layers of the its
32  Float_t MinimumPt=0.25;//Minimum Pt considered by the code
33   Float_t MaximumPt=50;
34   Float_t ZVertexCut=10;//in cm
35
36   //Options
37   Int_t RunOnAOD=1;
38   Int_t EfficiencyCorrection=1;//do efficiency corrections in this code
39   Int_t MakeMCHistos=0;//if 0 MC histograms are not made (will be empty if 1 and ran on real data)
40   Int_t DEBUG=0;//for debugging
41  
42   
43  
44   //Binning
45   Int_t nBinPhi=60;//Number of bins for #Delta#phi histograms
46   Int_t nBinEta=54;//Number of bins for #Delta#eta histograms
47   Int_t nBinsPhiEtaPhi=20;//Number of bins for #Delta#phi-#Delta#eta in #Delta#phi
48   Int_t nBinsPhiEtaEta=18;//Number of bins for #Delta#phi-#Delta#eta in #Delta#phi
49   Int_t nBinsPhiPhi=30;//Number of bins for #Delta#phi-#Delta#phi
50   Int_t nBinsEtaEta=27;//Number of bins for #Delta#eta-#Delta#eta
51   const Float_t fPi=3.1415926535898;
52   Float_t PhiPlotMin=-fPi/3;//Min bin edge in #Delta#phi
53   Float_t PhiPlotMax=2*fPi+PhiPlotMin;//Max bin edge
54   
55   //Size of some arrays change array contents below
56   if(IncludeLowPtBins)const Int_t NTriggerPtBins=11;
57   else const Int_t NTriggerPtBins=8;//max=20
58   const Int_t NEventsToMix=10;//max=100
59   const Int_t NCentralityBins=7;//max=10  //6
60   const Int_t PercentageCentralityBins=1;//0 or 1
61   const Int_t NAssociatedPtBins=25;//max=50
62   const Int_t N3ParticleAssociatedPtBins=4;//max=50
63   const Int_t NZVertexBinsForMixing=1;//max=20
64   const Int_t NXEBins=1;//max=20
65   const Int_t NumberOfTriggerIDs=1;
66   Float_t EffFitPtCut=3;
67
68   TF1 *EfficiencyFitLow=new TF1("EfficiencyFitLow","[0]/[1]*exp(-0.5*pow(x/[1],2))+[2]+[3]*x+[4]*x**2+[5]*x**3",MinimumPt,EffFitPtCut);
69   TF1 *EfficiencyFitHigh=new TF1("EfficiencyFitHigh","[0]+[1]*(x-3)",EffFitPtCut,MaximumPt);
70   const Int_t NParamFitLow=6;
71   const Int_t NParamFitHigh=2;
72   
73   //Not high enough occupancy to worry about the centrality in pp
74   //For overlapping centrality bins efficiencies from first bin are used
75   //For 1% bins from pass2 
76 Float_t FitLowParam[NCentralityBins*NParamFitLow]={
77   -0.0227546, 0.21379, 0.916013, 0.0586031, -0.0429047, 0.0064962,
78   -0.0224257, 0.215105, 0.916696, 0.0593096, -0.0429726, 0.00650312,
79   -0.0220935, 0.208194, 0.915207, 0.0613371, -0.0416018, 0.00576639,
80   -0.0229694, 0.208645, 0.919995, 0.0527738, -0.0344433, 0.00430404,
81   -0.0235676, 0.206099, 0.922875, 0.0476261, -0.0298983, 0.0032692,
82   -0.0238139, 0.206366, 0.925183, 0.0430703, -0.0264516, 0.00250388,
83   -0.024815, 0.204409, 0.928448, 0.03423, -0.019894, 0.00121746
84 };
85  Float_t FitHighParam[NCentralityBins*NParamFitHigh]={
86    0.881078, 0.00158217,
87    0.883456, 0.00144851,
88    0.880495, 0.00369636,
89    0.884535, 0.00297087,
90    0.884937, 0.00166201,
91    0.883934, 0.00290625,
92    0.884964, -0.000784607
93  };
94  
95   Float_t V2FitPtCut=2.75;
96   Float_t V3FitPtCut=2.75;
97   Float_t V4FitPtCut=2.75;
98   TF1 *V2FitLow=new TF1("V2FitLow","[0]*(x+[1]*x**2+[2]*x**3+[3]*x**4)",MinimumPt,V2FitPtCut);
99   TF1 *V2FitHigh=new TF1("V2FitHigh","[0]*exp(-0.5*((x-[1])/[2])**2)+[3]",V2FitPtCut,MaximumPt);
100   TF1 *V3FitLow=new TF1("V3FitLow","[0]*(x+[1]*x**2+[2]*x**3+[3]*x**4)",MinimumPt,V2FitPtCut);
101   TF1 *V3FitHigh=new TF1("V3FitHigh","[0]*exp(-0.5*((x-[1])/[2])**2)+[3]",V2FitPtCut,MaximumPt);
102  TF1 *V4FitLow=new TF1("V4FitLow","[0]*(x+[1]*x**2+[2]*x**3+[3]*x**4)",MinimumPt,V2FitPtCut);
103   TF1 *V4FitHigh=new TF1("V4FitHigh","[0]*exp(-0.5*((x-[1])/[2])**2)+[3]",V2FitPtCut,MaximumPt);
104
105   const Int_t NParamV2FitLow=4;
106   const Int_t NParamV2FitHigh=4;
107   const Int_t NParamV3FitLow=4;
108   const Int_t NParamV3FitHigh=4;
109   const Int_t NParamV4FitLow=4;
110   const Int_t NParamV4FitHigh=4;
111
112   //From FlowScale2 macro
113 Float_t FitLowParamV2[NCentralityBins*NParamV2FitLow]={
114 0.046931, -0.504895, 0.186345, -0.030504,
115 0.060953, -0.503882, 0.197361, -0.030849,
116 0.088939, -0.396433, 0.127467, -0.016969,
117 0.122602, -0.347615, 0.098626, -0.012576,
118 0.155403, -0.296297, 0.065681, -0.007492,
119 0.177546, -0.271958, 0.049388, -0.005644,
120 0.187634, -0.223967, 0.011281, 0.001462};
121 Float_t FitHighParamV2[NCentralityBins*NParamV2FitHigh]={
122 0.000000, 4.478095, 0.500000, 0.049007,
123 0.012659, 3.603740, 0.656049, 0.071827,
124 0.043836, 3.623140, 1.414293, 0.089569,
125 0.071277, 3.490653, 1.454361, 0.112755,
126 0.082163, 3.324474, 1.429585, 0.145809,
127 0.085348, 3.194132, 1.562318, 0.163847,
128 0.061056, 2.994950, 1.442946, 0.194451};
129 Float_t FitLowParamV3[NCentralityBins*NParamV3FitLow]={
130 0.019078, 0.793599, -0.436079, 0.075946,
131 0.025818, 0.088555, 0.070194, -0.026702,
132 0.028562, 0.297268, -0.112696, 0.013257,
133 0.034723, 0.137282, -0.008578, -0.009524,
134 0.041567, -0.002427, 0.093642, -0.032998,
135 0.045574, 0.059023, 0.029342, -0.019022,
136 0.052432, -0.067332, 0.087218, -0.028239};
137 Float_t FitHighParamV3[NCentralityBins*NParamV3FitHigh]={
138 0.079111, 3.802577, 3.182298, 0.000000,
139 0.105446, 4.216029, 2.377002, 0.000000,
140 0.106553, 4.500000, 2.825647, 0.009652,
141 0.099684, 3.875928, 1.968600, 0.022106,
142 0.103187, 3.597519, 1.980556, 0.022403,
143 0.123473, 3.773586, 2.099909, 0.011130,
144 0.140945, 3.618698, 1.852780, 0.000000};
145 Float_t FitLowParamV4[NCentralityBins*NParamV4FitLow]={
146 0.005289, 1.000000, 0.398581, -0.095251,
147 0.018611, -1.000000, 0.837275, -0.163058,
148 0.022107, -1.000000, 0.847701, -0.168469,
149 0.025434, -1.000000, 0.854238, -0.176690,
150 0.030768, -1.000000, 0.811683, -0.164335,
151 0.037000, -1.000000, 0.781250, -0.158122,
152 0.043484, -1.000000, 0.752546, -0.153147};
153   Float_t FitHighParamV4[NCentralityBins*NParamV4FitHigh]={
154 0.061365, 3.669890, 1.724126, 0.014456,
155 0.086516, 4.194769, 1.723464, 0.000000,
156 0.030263, 4.008413, 0.563914, 0.067831,
157 0.073456, 4.255665, 1.551659, 0.025494,
158 0.058771, 4.102898, 1.191926, 0.048796,
159 0.054682, 4.221963, 1.119594, 0.064515,
160 0.094564, 4.500000, 1.568750, 0.038027};
161
162
163
164   if(IncludeLowPtBins) Float_t TriggerPtBins[(NTriggerPtBins+1)]={0.75,1,2,2.5,3,4,6,8,10,15,20,25};
165   else  Float_t TriggerPtBins[(NTriggerPtBins+1)]={2.5,3,4,6,8,10,15,20,25};
166
167   Float_t AssociatedPtBins[(NAssociatedPtBins+1)]={0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12,15,20,25,30,40,50,70,100};
168  
169   Float_t AssociatedPtBins31[N3ParticleAssociatedPtBins]={0.5,0.75,1,2};
170   Float_t AssociatedPtBins32[N3ParticleAssociatedPtBins]={0.75,1,2,3};
171  
172   Int_t CentralityBins1[NCentralityBins]={0,0,5,10,20,30,40};
173   Int_t CentralityBins2[NCentralityBins]={2,5,10,20,30,40,50};
174
175   Float_t XEBins[(NXEBins+1)]={0,0.01};
176   //char *TriggerIDArray="CINT1B";//seperate multiple with ,
177   //char *TriggerIDArray="CSH1-B";//seperate multiple with ,
178   char *TriggerIDArray="C";//PbPb test, using SelectCollisionsCanidates instead
179
180   ////////////////////////
181   //Add the task
182   ////////////////////////
183   AliAnalysisTaskDiHadron *task = new AliAnalysisTaskDiHadron("julery_DiHadron");
184   task->SetCuts(MinimumClustersTPC,MinClusterRatio,MaxTPCchi2,MinimumClustersITS, EtaCut,TriggerEtaCut,NearPhiCut,XECut,MaxDCA,MaxDCAXY,MaxDCAZ, DCAMethod, TPCRefit,ITSRefit,SPDCut,MinimumPt,MaximumPt,ZVertexCut,NumberOfTriggerIDs,TriggerIDArray);
185   task->SetOptions(RunOnAOD,EfficiencyCorrection,DEBUG,MakeMCHistos);
186   task->SetBins(nBinPhi,nBinEta,nBinsPhiEtaPhi,nBinsPhiEtaEta,nBinsPhiPhi,nBinsEtaEta,PhiPlotMin,PhiPlotMax,NTriggerPtBins,NEventsToMix,NCentralityBins,PercentageCentralityBins,NAssociatedPtBins,N3ParticleAssociatedPtBins,NZVertexBinsForMixing,NXEBins,TriggerPtBins,AssociatedPtBins,AssociatedPtBins31,AssociatedPtBins32,CentralityBins1,CentralityBins2,XEBins);
187   task->SetEfficiencies(EffFitPtCut,EfficiencyFitLow,EfficiencyFitHigh,NParamFitLow,NParamFitHigh,FitLowParam,FitHighParam);
188   task->SetFlow(V2FitPtCut,V3FitPtCut,V4FitPtCut,V2FitLow,V2FitHigh,V3FitLow,V3FitHigh,V4FitLow,V4FitHigh,NParamV2FitLow,NParamV2FitHigh,NParamV3FitLow,NParamV3FitHigh,NParamV4FitLow,NParamV4FitHigh,FitLowParamV2,FitHighParamV2,FitLowParamV3,FitHighParamV3,FitLowParamV4,FitHighParamV4);
189
190 // physics selection
191 Int_t isMC=0;//1 for MC 0 for DATA
192 //gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
193 //AliPhysicsSelectionTask *PhysicsTask=AddTaskPhysicsSelection(isMC, 0); //isMC is true when processing monte carlo, the second 0 disables the cluster vs tracklets
194   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
195   //AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isMC,0);
196   task->SelectCollisionCandidates(AliVEvent::kMB);
197   //AliCentralitySelectionTask *centSelTask = AliCentralitySelectionTask("CentralitySelection");
198  
199   mgr->AddTask(task);
200  
201
202   //================================================
203   //              data containers
204   //================================================
205   //            find input container
206   //below the trunk version
207   //   AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
208   //this is the old way!!!
209   AliAnalysisDataContainer *cinput  = (AliAnalysisDataContainer*)mgr->GetContainers()->FindObject("cAUTO_INPUT");
210   
211   //            define output containers, please use 'username'_'somename'
212   AliAnalysisDataContainer *coutput1 = 
213     mgr->CreateContainer("julery_DiHadron", TList::Class(),
214                              AliAnalysisManager::kOutputContainer,"julery_DiHadron.root");
215
216   //           connect containers
217   mgr->ConnectInput  (task,  0, cinput );
218   mgr->ConnectOutput (task,  0, coutput1);
219
220   return task;
221   
222 }