AliPhysics  608b256 (608b256)
MakeCuts4Charm4Prong.C
Go to the documentation of this file.
1 //gSystem->AddIncludePath("-I/$ALICE_ROOT/PWG3/vertexingHF");
2 
3 #include <Riostream.h>
4 #include <TFile.h>
5 //#include <AliRDHFCutsD0toKpipipi.h>
6 #include <TClonesArray.h>
7 #include <TParameter.h>
8 
9 //Author: Chiara Bianchin, cbianchi@pd.infn.it, Fabio Colamaria, fabio.colamaria@ba.infn.it
10 
11 //macro to make a .root file which contains an AliRDHFCutsD0toKpipipi for AliAnalysisTaskSESelectHF4Prong task
12 
14 
15 gSystem->Load("libTree");
16 gSystem->Load("libGeom");
17 gSystem->Load("libPhysics");
18 gSystem->Load("libVMC");
19 gSystem->Load("libMinuit");
20 gSystem->Load("libSTEERBase");
21 gSystem->Load("libESD");
22 gSystem->Load("libAOD");
23 gSystem->Load("libANALYSIS");
24 gSystem->Load("libANALYSISalice");
25 gSystem->Load("libCORRFW");
26 gSystem->Load("libPWGHFbase");
27 gSystem->Load("libPWGHFvertexingHF");
28 gSystem->Load("libPWGmuon");
29 
30  AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
31  RDHFCharm4Prong->SetName("Charm4ProngCuts");
32  RDHFCharm4Prong->SetTitle("Cuts for D0 analysis");
33 
34  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
35  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
36  //default
37  esdTrackCuts->SetRequireTPCRefit(kTRUE);
38  esdTrackCuts->SetRequireITSRefit(kTRUE);
39  esdTrackCuts->SetMinNClustersITS(4); // default is 5
40  //esdTrackCuts->SetMinNClustersTPC(70);
41  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
42  AliESDtrackCuts::kAny);
43  // default is kBoth, otherwise kAny
44  esdTrackCuts->SetMinDCAToVertexXY(0.);
45  esdTrackCuts->SetPtRange(0.3,1.e10);
46 
47  RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
48 
49 
50 
51  const Int_t nvars=9;
52 
53  const Int_t nptbins=5;
54  Float_t* ptbins;
55  ptbins=new Float_t[nptbins+1];
56  ptbins[0]=0.;
57  ptbins[1]=4.5;
58  ptbins[2]=6.;
59  ptbins[3]=8.;
60  ptbins[4]=12.;
61  ptbins[5]=16.;
62 
63  RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins);
64 
65 
66  Float_t** rdcutsvalmine;
67  rdcutsvalmine=new Float_t*[nvars];
68  for(Int_t iv=0;iv<nvars;iv++){
69  rdcutsvalmine[iv]=new Float_t[nptbins];
70  }
71 
72  //setting my cut values
73  //cuts order
74  // printf(" |M-MD0| [GeV]");
75  // printf(" DCA [cm]");
76  // printf(" Dist 2-trk Vtx to PrimVtx [cm]");
77  // printf(" Dist 3-trk Vtx to PrimVtx [cm]");
78  // printf(" Dist 4-trk Vtx to PrimVtx [cm]");
79  // printf(" cosThetaPoint");
80  // printf(" pt [GeV/c]");
81  // printf(" rho mass");
82  // printf(" PID cut");
83 
84  //setting my cut values
85  rdcutsvalmine[0][0]=0.2;
86  rdcutsvalmine[1][0]=0.025;
87  rdcutsvalmine[2][0]=0.0900;
88  rdcutsvalmine[3][0]=0.0900;
89  rdcutsvalmine[4][0]=0.0600;
90  rdcutsvalmine[5][0]=0.995;
91  rdcutsvalmine[6][0]=0.;
92  rdcutsvalmine[7][0]=0.1;
93  rdcutsvalmine[8][0]=0.;
94 
95  rdcutsvalmine[0][1]=0.2;
96  rdcutsvalmine[1][1]=0.015;
97  rdcutsvalmine[2][1]=0.1200;
98  rdcutsvalmine[3][1]=0.1000;
99  rdcutsvalmine[4][1]=0.0850;
100  rdcutsvalmine[5][1]=0.99;
101  rdcutsvalmine[6][1]=0.;
102  rdcutsvalmine[7][1]=0.1;
103  rdcutsvalmine[8][1]=0.;
104 
105  rdcutsvalmine[0][2]=0.2;
106  rdcutsvalmine[1][2]=0.025;
107  rdcutsvalmine[2][2]=0.0975;
108  rdcutsvalmine[3][2]=0.0900;
109  rdcutsvalmine[4][2]=0.0800;
110  rdcutsvalmine[5][2]=0.99;
111  rdcutsvalmine[6][2]=0.;
112  rdcutsvalmine[7][2]=0.1;
113  rdcutsvalmine[8][2]=0.;
114 
115  rdcutsvalmine[0][3]=0.2;
116  rdcutsvalmine[1][3]=0.015;
117  rdcutsvalmine[2][3]=0.0975;
118  rdcutsvalmine[3][3]=0.0700;
119  rdcutsvalmine[4][3]=0.0750;
120  rdcutsvalmine[5][3]=0.9825;
121  rdcutsvalmine[6][3]=0.;
122  rdcutsvalmine[7][3]=0.1;
123  rdcutsvalmine[8][3]=0.;
124 
125  rdcutsvalmine[0][4]=rdcutsvalmine[0][3];
126  rdcutsvalmine[1][4]=rdcutsvalmine[1][3];
127  rdcutsvalmine[2][4]=rdcutsvalmine[2][3];
128  rdcutsvalmine[3][4]=rdcutsvalmine[3][3];
129  rdcutsvalmine[4][4]=rdcutsvalmine[4][3];
130  rdcutsvalmine[5][4]=rdcutsvalmine[5][3];
131  rdcutsvalmine[6][4]=rdcutsvalmine[6][3];
132  rdcutsvalmine[7][4]=rdcutsvalmine[7][3];
133  rdcutsvalmine[8][4]=rdcutsvalmine[8][3];
134 
135  RDHFCharm4Prong->SetCuts(nvars,nptbins,rdcutsvalmine);
136  cout<<"This is the odject I'm going to save:"<<endl;
137  RDHFCharm4Prong->PrintAll();
138  TFile* fout=new TFile("Charm4ProngCutsDef_16GeV.root","recreate"); //set this!!
139  fout->cd();
140  RDHFCharm4Prong->Write();
141  fout->Close();
142 
143 }
144 
146 
147 gSystem->Clear();
148 gSystem->Load("libTree");
149 gSystem->Load("libGeom");
150 gSystem->Load("libPhysics");
151 gSystem->Load("libVMC");
152 gSystem->Load("libMinuit");
153 gSystem->Load("libSTEERBase");
154 gSystem->Load("libESD");
155 gSystem->Load("libAOD");
156 gSystem->Load("libANALYSIS");
157 gSystem->Load("libANALYSISalice");
158 gSystem->Load("libCORRFW");
159 gSystem->Load("libPWGHFbase");
160 gSystem->Load("libPWGHFvertexingHF");
161 gSystem->Load("libPWGmuon");
162 
163  AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
164  RDHFCharm4Prong->SetName("loosercuts");
165  RDHFCharm4Prong->SetTitle("Cuts for significance maximization");
166 
167  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
168  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
169  //default
170  esdTrackCuts->SetRequireTPCRefit(kTRUE);
171  esdTrackCuts->SetRequireITSRefit(kTRUE);
172  //esdTrackCuts->SetMinNClustersITS(4);
173  //esdTrackCuts->SetMinNClustersTPC(120);
174 
175  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
176  esdTrackCuts->SetMinDCAToVertexXY(0.);
177  esdTrackCuts->SetEtaRange(-0.8,0.8);
178  esdTrackCuts->SetPtRange(0.3,1.e10);
179 
180  RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
181 
182  const Int_t nvars=9;
183 
184  const Int_t nptbins=5; //change this when adding pt bins!
185  Float_t ptbins[nptbins+1];
186  ptbins[0]=0.;
187  ptbins[1]=4.5;
188  ptbins[2]=6.;
189  ptbins[3]=8.;
190  ptbins[4]=12.;
191  ptbins[5]=16.; //o 25!
192 
193  RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins);
194 
195  Float_t** rdcutsvalmine;
196  rdcutsvalmine=new Float_t*[nvars];
197  for(Int_t iv=0;iv<nvars;iv++){
198  rdcutsvalmine[iv]=new Float_t[nptbins];
199  }
200 
201  //setting my cut values
202  //cuts order
203  // printf(" |M-MD0| [GeV]");
204  // printf(" DCA [cm]");
205  // printf(" Dist 2-trk Vtx to PrimVtx [cm]");
206  // printf(" Dist 3-trk Vtx to PrimVtx [cm]");
207  // printf(" Dist 4-trk Vtx to PrimVtx [cm]");
208  // printf(" cosThetaPoint");
209  // printf(" pt [GeV/c]");
210  // printf(" rho mass");
211  // printf(" PID cut");
212 
213  Float_t cutsMatrixD0toKpipipiStand[nptbins][nvars]={{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
214 
215 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
216 
217 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
218 
219 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
220 
221 {0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.}};
222 
223  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
224  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
225  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
226  for (Int_t ibin=0;ibin<nptbins;ibin++){
227  for (Int_t ivar = 0; ivar<nvars; ivar++){
228  cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpipipiStand[ibin][ivar];
229  }
230  }
231  RDHFCharm4Prong->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
232 
233  Int_t nvarsforopt=RDHFCharm4Prong->GetNVarsForOpt();
234  Int_t dim=5; //set this!!
235  Bool_t *boolforopt;
236  boolforopt=new Bool_t[nvars];
237  if(dim>nvarsforopt){
238  cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
239  return;
240  } else {
241  if(dim==nvarsforopt){
242  boolforopt=RDHFCharm4Prong->GetVarsForOpt();
243  }else{
244  TString *names;
245  names=new TString[nvars];
246  TString answer="";
247  Int_t checktrue=0;
248  names=RDHFCharm4Prong->GetVarNames();
249  for(Int_t i=0;i<nvars;i++){
250  cout<<names[i]<<" for opt? (y/n)"<<endl;
251  cin>>answer;
252  if(answer=="y") {
253  boolforopt[i]=kTRUE;
254  checktrue++;
255  }
256  else boolforopt[i]=kFALSE;
257  }
258  if (checktrue!=dim) {
259  cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
260  return;
261  }
262  RDHFCharm4Prong->SetVarsForOpt(dim,boolforopt);
263  }
264  }
265 
266 // Tight values for variables to be maximized
267 // (DCA, Dist12, Dist3, Dist4, CosThetaPoint)
268 // Indexes: variable, ptbin
269 
270  Float_t tighterval[5][5];
271 
272  //number of steps for each variable is set in the AddTask arguments (default=8)
273 
274  tighterval[0][0]=0.025;
275  tighterval[1][0]=0.09;
276  tighterval[2][0]=0.09;
277  tighterval[3][0]=0.06;
278  tighterval[4][0]=0.955;
279 
280  tighterval[0][1]=0.015;
281  tighterval[1][1]=0.12;
282  tighterval[2][1]=0.10;
283  tighterval[3][1]=0.085;
284  tighterval[4][1]=0.99;
285 
286  tighterval[0][2]=0.025;
287  tighterval[1][2]=0.0975;
288  tighterval[2][2]=0.09;
289  tighterval[3][2]=0.08;
290  tighterval[4][2]=0.99;
291 
292  tighterval[0][3]=0.015;
293  tighterval[1][3]=0.0975;
294  tighterval[2][3]=0.07;
295  tighterval[3][3]=0.075;
296  tighterval[4][3]=0.9825;
297 
298  tighterval[0][4]=0.015;
299  tighterval[1][4]=0.0975;
300  tighterval[2][4]=0.07;
301  tighterval[3][4]=0.075;
302  tighterval[4][4]=0.9825;
303 
304  TString name="";
305  Int_t arrdim=dim*nptbins;
306  cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
307  TClonesArray max("TParameter<float>",arrdim);
308  for(Int_t ival=0;ival<dim;ival++){
309  for(Int_t jpt=0;jpt<nptbins;jpt++){
310  name=Form("par%dptbin%d",ival,jpt);
311  cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
312  new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
313  }
314  }
315 
316  Bool_t flagPID=kFALSE;
317  RDHFCharm4Prong->SetUsePID(flagPID);
318 
319  RDHFCharm4Prong->PrintAll();
320  printf("Use PID? %s\n",flagPID ? "yes" : "no");
321 
322 /*
323  //pid settings
324  AliAODPidHF* pidObj=new AliAODPidHF();
325  //pidObj->SetName("pid4D0");
326  Int_t mode=1;
327  const Int_t nlims=2;
328  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
329  Bool_t compat=kTRUE; //effective only for this mode
330  Bool_t asym=kTRUE;
331  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
332  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
333  pidObj->SetMatch(mode);
334  pidObj->SetPLimit(plims,nlims);
335  pidObj->SetSigma(sigmas);
336  pidObj->SetCompat(compat);
337  pidObj->SetTPC(kTRUE);
338  pidObj->SetTOF(kTRUE);
339  RDHFCharm4Prong->SetPidHF(pidObj);
340 
341  RDHFCharm4Prong->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
342 */
343 
344  //activate pileup rejection (for pp)
346 
347  //Do not recalculate the vertex
348  RDHFCharm4Prong->SetRemoveDaughtersFromPrim(kTRUE); //activate for pp
349 
350 /* Per il pp levo la centralit√†!
351 //Metto kCentOff! Per il Pb c'era kCentV0M
352  TString cent="";
353  //centrality selection (Pb-Pb)
354  Float_t minc=20,maxc=80;
355  RDHFCharm4Prong->SetMinCentrality(minc);
356  RDHFCharm4Prong->SetMaxCentrality(maxc);
357  cent=Form("%.0f%.0f",minc,maxc);
358  RDHFCharm4Prong->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
359 */
360 
361  //temporary
362 // RDHFCharm4Prong->SetFixRefs();
363 
364  TFile* fout=new TFile("Charm4ProngCutsForMaxim.root","recreate");
365  fout->cd();
366  RDHFCharm4Prong->Write();
367  max.Write();
368  fout->Close();
369 }
TSystem * gSystem
TString * GetVarNames() const
Definition: AliRDHFCuts.h:265
int Int_t
Definition: External.C:63
void SetCuts(Int_t nVars, Int_t nPtBins, Float_t **cutsRD)
void MakeCuts4Charm4Prong()
float Float_t
Definition: External.C:68
Bool_t * GetVarsForOpt() const
Definition: AliRDHFCuts.h:266
void MakeCuts4Charm4ProngForMaxim()
Int_t GetNVarsForOpt() const
Definition: AliRDHFCuts.h:267
void SetVarsForOpt(Int_t nVars, Bool_t *forOpt)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:219
virtual void PrintAll() const
void SetRemoveDaughtersFromPrim(Bool_t removeDaughtersPrim)
Definition: AliRDHFCuts.h:229
void SetPtBins(Int_t nPtBinLimits, Float_t *ptBinLimits)
void AddTrackCuts(const AliESDtrackCuts *cuts)
Definition: AliRDHFCuts.h:217
bool Bool_t
Definition: External.C:53
TFile * fout
input train file
void SetOptPileup(Int_t opt=0)
Definition: AliRDHFCuts.h:233
Int_t nptbins