]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/hbtanalysis.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AntiMerging cut implemented
[u/mrichter/AliRoot.git] / HBTAN / hbtanalysis.C
CommitLineData
f1e7ff24 1 #include "AliHBTAnalysis.h"
6e0d7f90 2// #include "alles.h"
f1e7ff24 3 #include "AliHBTReader.h"
4 #include "AliHBTReaderKineTree.h"
5 #include "AliHBTReaderITSv2.h"
6 #include "AliHBTReaderITSv1.h"
7 #include "AliHBTReaderTPC.h"
8 #include "AliHBTReaderInternal.h"
9 #include "AliHBTParticleCut.h"
10 #include "AliHBTEvent.h"
11 #include "AliHBTPairCut.h"
12 #include "AliHBTQResolutionFctns.h"
13 #include "AliHBTQDistributionFctns.h"
14 #include "AliHBTTwoTrackEffFctn.h"
15 #include "AliHBTCorrelFctn.h"
16 #include "TSystem.h"
17 #include "TObjString.h"
18 #include "TString.h"
6e0d7f90 19 #include "TPDGCode.h"
f1e7ff24 20 #include "AliHBTMonDistributionFctns.h"
21 #include "AliHBTMonResolutionFctns.h"
22
e477958c 23void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
24 Int_t first = -1,Int_t last = -1,
25 char *outfile = "hbtanalysis.root")
0a6dde4c 26 {
9616170a 27
28 AliHBTTrackPoints::SetDebug(1);
0a6dde4c 29//HBT Anlysis Macro
30//Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
e477958c 31
32//datatype defines type of data to be read
33// Kine - analyzes Kine Tree: simulated particles
34// TPC - analyzes TPC tracking + particles corresponding to these tracks
35// ITSv1 - analyzes ITSv1 ----------------------//--------------------------
36// ITSv2 - analyzes ITSv2 ----------------------//--------------------------
37
38//processopt defines option passed to AliHBTAnlysis::Process method
39// default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them
40// Tracks - process only recontructed tracks
41// Particles - process only simulated particles
42
0a6dde4c 43//Reads data from diroctories from first to last(including)
44// For examples if first=3 and last=5 it reads from
45// ./3/
46// ./4/
47// ./5/
48//if first or last is negative (or both), it reads from current directory
49//
50//these names I use when analysis is done directly from CASTOR files via RFIO
51// const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
52// const char* serie="standard";
53// const char* field = "0.4";
9616170a 54
0a6dde4c 55 const char* basedir=".";
56 const char* serie="";
57 const char* field = "";
58
f1e7ff24 59 Bool_t threeDcuts = kFALSE;
60
0a6dde4c 61 // Dynamically link some shared libs
f1e7ff24 62// gROOT->LoadMacro("loadlibs.C");
63// loadlibs();
0a6dde4c 64
71094efa 65 gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
0a6dde4c 66
67
68 /***********************************************************/
69 //Create analysis and reader
70 AliHBTAnalysis * analysis = new AliHBTAnalysis();
71
72 AliHBTReader* reader;
e477958c 73 Int_t kine = strcmp(datatype,"Kine");
6e0d7f90 74 Int_t ESD = strcmp(datatype,"ESD");
e477958c 75 Int_t TPC = strcmp(datatype,"TPC");
76 Int_t ITSv1 = strcmp(datatype,"ITSv1");
77 Int_t ITSv2 = strcmp(datatype,"ITSv2");
91923a38 78 Int_t intern = strcmp(datatype,"Intern");
79
e477958c 80 if(!kine)
81 {
82 reader = new AliHBTReaderKineTree();
83 processopt="Particles"; //this reader by definition reads only simulated particles
84 }
6e0d7f90 85 else if(!ESD)
86 {
87 AliHBTReaderESD* esdreader = new AliHBTReaderESD();
88 esdreader->ReadParticles(kTRUE);
89 reader = esdreader;
90 }
e477958c 91 else if(!TPC)
92 {
93 reader = new AliHBTReaderTPC();
9616170a 94 ((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
e477958c 95 }
96 else if(!ITSv1)
97 {
98 reader = new AliHBTReaderITSv1();
99 }
100 else if(!ITSv2)
101 {
102 reader = new AliHBTReaderITSv2();
103 }
91923a38 104 else if(!intern)
105 {
6e0d7f90 106 reader = new AliHBTReaderInternal("data.root");
91923a38 107 }
e477958c 108 else
109 {
110 cerr<<"Option "<<datatype<<" not recognized. Exiting"<<endl;
111 return;
112 }
113
0a6dde4c 114 TObjArray* dirs=0;
e477958c 115 if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
116 {//read from many dirs dirs
117 char buff[50];
118 dirs = new TObjArray(last-first+1);
119 for (Int_t i = first; i<=last; i++)
120 {
121 sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
122 TObjString *odir= new TObjString(buff);
123 dirs->Add(odir);
124 }
0a6dde4c 125 }
e477958c 126 reader->SetDirs(dirs);
0a6dde4c 127
128 /***********************************************************/
f1e7ff24 129 /***** R E A D E R ***************************************/
130 /***********************************************************/
0a6dde4c 131
132 //we want have only low pt pi+ so set a cut to reader
133 AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
6e0d7f90 134 readerpartcut->SetPtRange(0.0,10.0);
0a6dde4c 135 readerpartcut->SetPID(kPiPlus);
136 reader->AddParticleCut(readerpartcut);//read this particle type with this cut
137
6e0d7f90 138// readerpartcut->SetPtRange(0.0,10.0);
139// readerpartcut->SetPID(kPiMinus);
140// reader->AddParticleCut(readerpartcut);//read this particle type with this cut
141
0a6dde4c 142 analysis->SetReader(reader);
f1e7ff24 143
9616170a 144 AliHBTPairCut *paircut = new AliHBTPairCut();
145 Float_t qinvmin = 0.0;
146 Float_t qinvmax = 0.05;//50MeV
147 paircut->SetQInvRange(qinvmin,qinvmax);
148
149 // paircut->SetAvSeparationRange(10.,10000.);
150
151// AliHBTParticleCut* partcut= new AliHBTParticleCut();
152// partcut->SetPID(kPiPlus);
153// paircut->SetFirstPartCut(partcut);
154// partcut->SetPID(kPiMinus);
155// paircut->SetSecondPartCut(partcut);
156
157 analysis->SetGlobalPairCut(paircut);
158 analysis->SetAntiMergingCut(5.);//in cm
159
16f9289f 160 /***********************************************************/
161 /***** W E I G H T S ******************************/
162 /***********************************************************/
163
164 AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
165 AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
166 AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
167
168 AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
169//AliHBTLLWeights::Instance()->SetR1dw(8.0);
170
171 AliHBTLLWeights::Instance()->Init();
172 AliHBTLLWeights::Instance()->Set();
173
174 //example function
9616170a 175 AliHBTWeightTheorQInvFctn *wqinvcfP = new AliHBTWeightTheorQInvFctn(100,qinvmax);
176 wqinvcfP->SetNumberOfBinsToScale(30);
177 wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function");
16f9289f 178 analysis->AddParticleFunction(wqinvcfP);
179
f1e7ff24 180 /************************************************************/
181 /**** Q INV Correlation Function ************************/
0a6dde4c 182 /************************************************************/
6e0d7f90 183
0a6dde4c 184
f1e7ff24 185 AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
186 AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
0a6dde4c 187
188 analysis->AddTrackFunction(qinvcfT);
189 analysis->AddParticleFunction(qinvcfP);
190
0a6dde4c 191
192 qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
193 qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
91923a38 194
f1e7ff24 195 /************************************************************/
196 /**** Q OUT Correlation Function ************************/
197 /************************************************************/
198
91923a38 199 AliHBTQOutCMSLCCorrelFctn* qoutP = new AliHBTQOutCMSLCCorrelFctn();
200 qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
201 AliHBTQOutCMSLCCorrelFctn* qoutT = new AliHBTQOutCMSLCCorrelFctn();
202 qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
203
204 AliHBTPairCut *outPairCut = new AliHBTPairCut();
205 outPairCut->SetQOutCMSLRange(0.0,0.15);
206 outPairCut->SetQSideCMSLRange(0.0,0.02);
207 outPairCut->SetQLongCMSLRange(0.0,0.02);
208 qoutP->SetPairCut(outPairCut);
209 qoutT->SetPairCut(outPairCut);
f1e7ff24 210
6e0d7f90 211 //analysis->AddTrackFunction(qoutP);
212 //analysis->AddParticleFunction(qoutT);
f1e7ff24 213
214 /************************************************************/
215 /**** Q SIDE Correlation Function ***********************/
216 /************************************************************/
91923a38 217
f1e7ff24 218 AliHBTQSideCMSLCCorrelFctn* qsideP = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax);
91923a38 219 qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
f1e7ff24 220 AliHBTQSideCMSLCCorrelFctn* qsideT = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax);
91923a38 221 qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
222
223 AliHBTPairCut *sidePairCut = new AliHBTPairCut();
224 sidePairCut->SetQOutCMSLRange(0.0,0.02);
225 sidePairCut->SetQSideCMSLRange(0.0,0.15);
226 sidePairCut->SetQLongCMSLRange(0.0,0.02);
227 qsideP->SetPairCut(sidePairCut);
228 qsideT->SetPairCut(sidePairCut);
0a6dde4c 229
6e0d7f90 230 //analysis->AddTrackFunction(qsideP);
231 //analysis->AddParticleFunction(qsideT);
f1e7ff24 232
233 /************************************************************/
234 /**** Q LONG Correlation Function ***********************/
235 /************************************************************/
91923a38 236
f1e7ff24 237 AliHBTQLongCMSLCCorrelFctn* qlongP = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax);
91923a38 238 qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
f1e7ff24 239 AliHBTQLongCMSLCCorrelFctn* qlongT = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax);
91923a38 240 qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
241
242 AliHBTPairCut *longPairCut = new AliHBTPairCut();
243 longPairCut->SetQOutCMSLRange(0.0,0.02);
244 longPairCut->SetQSideCMSLRange(0.0,0.02);
245 longPairCut->SetQLongCMSLRange(0.0,0.15);
246 qlongP->SetPairCut(longPairCut);
247 qlongT->SetPairCut(longPairCut);
248
6e0d7f90 249 //analysis->AddTrackFunction(qlongT);
250 //analysis->AddParticleFunction(qlongT);
f1e7ff24 251
252
253 /************************************************************/
254 /************************************************************/
255 /************************************************************/
256 /**** R E S O L U T I O N S ***********************/
257 /************************************************************/
258 /************************************************************/
259 /************************************************************/
260
261
262
263 AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015); //QInvCMSLC Res Vs Kt
264 if (threeDcuts) qinvVsktF.SetPairCut(paircut);
6e0d7f90 265 //analysis->AddResolutionFunction(&qinvVsktF);
f1e7ff24 266
267 AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05); //QOutCMSLC Res Vs Kt
268 if (threeDcuts) qoutVsktF.SetPairCut(outPairCut);
6e0d7f90 269 //analysis->AddResolutionFunction(&qoutVsktF);
f1e7ff24 270
271 AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015); //QSideCMSLC Res Vs Kt
272 if (threeDcuts) qsideVsktF.SetPairCut(sidePairCut);
6e0d7f90 273 //analysis->AddResolutionFunction(&qsideVsktF);
f1e7ff24 274
275 AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015); //QLongCMSLC Res Vs Kt
276 if (threeDcuts) qlongVsktF.SetPairCut(longPairCut);
6e0d7f90 277 //analysis->AddResolutionFunction(&qlongVsktF);
f1e7ff24 278
279 AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvCMSLC Res Vs QInvCMSLC
6e0d7f90 280 //analysis->AddResolutionFunction(&qInvVsqinvF);
f1e7ff24 281
282 AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05); //QOutCMSLC Res Vs QInvCMSLC
283 if (threeDcuts) qoutVsqinvF.SetPairCut(outPairCut);
6e0d7f90 284 //analysis->AddResolutionFunction(&qoutVsqinvF);
f1e7ff24 285
286 AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideCMSLC Res Vs QInvCMSLC
287 if (threeDcuts) qsideVsqinvF.SetPairCut(sidePairCut);
6e0d7f90 288 //analysis->AddResolutionFunction(&qsideVsqinvF);
f1e7ff24 289
290 AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongCMSLC Res Vs QInvCMSLC
291 if (threeDcuts) qlongVsqinvF.SetPairCut(longPairCut);
6e0d7f90 292 //analysis->AddResolutionFunction(&qlongVsqinvF);
f1e7ff24 293
294
295 AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
296 if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
6e0d7f90 297 //analysis->AddResolutionFunction(&pairThetaVsqinv);
f1e7ff24 298
299 AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
300 if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
6e0d7f90 301 //analysis->AddResolutionFunction(&pairThetaVsKt);
f1e7ff24 302
303 AliHBTPairCut phipc;
304 phipc.SetQLongCMSLRange(0.0,0.02);
305
306 AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
307 if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
6e0d7f90 308 //analysis->AddResolutionFunction(&pairPhiVsqinv);
f1e7ff24 309
310 AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
311 if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
6e0d7f90 312 //analysis->AddResolutionFunction(&pairPhiVsKt);
f1e7ff24 313
314
315 AliHBTQOutResolVsQOutFctn qoutVsqoutF; //QOutCMSLC Res Vs QOut
316 if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
6e0d7f90 317 //analysis->AddResolutionFunction(&qoutVsqoutF);
f1e7ff24 318
319 AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideCMSLC Res Vs QSide
320 if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
6e0d7f90 321 //analysis->AddResolutionFunction(&qsideVsqsideF);
f1e7ff24 322
91923a38 323
f1e7ff24 324 AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongCMSLC Res Vs QLong
325 if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
6e0d7f90 326 //analysis->AddResolutionFunction(&qlongVsqlongF);
f1e7ff24 327
328
329 /************************************************************/
330 /************************************************************/
331 /************************************************************/
332 /**** D I S T R B U T I O N S ***********************/
333 /************************************************************/
334 /************************************************************/
335 /************************************************************/
91923a38 336
f1e7ff24 337 /************************************************************/
338 /**** OUT VS KT ***********************/
339 /************************************************************/
340 AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
341 qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
342 AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
343 qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
344 if (threeDcuts) qoutVsKtDistP.SetPairCut(outPairCut);
345 if (threeDcuts) qoutVsKtDistT.SetPairCut(outPairCut);
6e0d7f90 346 //analysis->AddParticleFunction(&qoutVsKtDistP);
347 //analysis->AddTrackFunction(&qoutVsKtDistT);
f1e7ff24 348 /************************************************************/
349 /**** Side VS KT ***********************/
350 /************************************************************/
351 AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
352 qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
353 AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
354 qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
355 if (threeDcuts) qSideVsKtDistP.SetPairCut(sidePairCut);
356 if (threeDcuts) qSideVsKtDistT.SetPairCut(sidePairCut);
6e0d7f90 357 //analysis->AddParticleFunction(&qSideVsKtDistP);
358 //analysis->AddTrackFunction(&qSideVsKtDistT);
f1e7ff24 359 /************************************************************/
360 /**** Long VS KT ***********************/
361 /************************************************************/
362 AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
363 qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
364 AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
365 qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
366 if (threeDcuts) qLongVsKtDistP.SetPairCut(longPairCut);
367 if (threeDcuts) qLongVsKtDistT.SetPairCut(longPairCut);
6e0d7f90 368 //analysis->AddParticleFunction(&qLongVsKtDistP);
369 //analysis->AddTrackFunction(&qLongVsKtDistT);
f1e7ff24 370
371 /************************************************************/
372 /************************************************************/
373 /************************************************************/
374 /**** M O N I T O R R E S O L U T I O N S *************/
375 /************************************************************/
376 /************************************************************/
377 /************************************************************/
378
379 AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90 380 //analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
f1e7ff24 381
382 AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90 383 //analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
f1e7ff24 384
385 AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90 386 //analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);
f1e7ff24 387
388 AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
6e0d7f90 389 //analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);
f1e7ff24 390
6e0d7f90 391 AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
f1e7ff24 392 analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);
393
394 AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
395 analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);
396
397 AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
398 analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
399
400
401
402
403 /************************************************************/
404 /**** P R O C E S S ***********************/
405 /************************************************************/
6e0d7f90 406 analysis->SetBufferSize(2);
91923a38 407 analysis->Process(processopt);
408
0a6dde4c 409 TFile histoOutput(outfile,"recreate");
410 analysis->WriteFunctions();
411 histoOutput.Close();
412
413 delete qinvcfP;
414 delete qinvcfT;
d4ebbc5f 415 delete paircut;
0a6dde4c 416 delete readerpartcut;
417 if (dirs)
418 {
419 dirs->SetOwner();
420 delete dirs;
421 }
422 delete reader;
423 delete analysis;
424
425 }
426