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