1 #include "AliHBTAnalysis.h"
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"
17 #include "TObjString.h"
20 #include "AliHBTMonDistributionFctns.h"
21 #include "AliHBTMonResolutionFctns.h"
23 void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
24 Int_t first = -1,Int_t last = -1,
25 char *outfile = "hbtanalysis.root")
28 AliHBTTrackPoints::SetDebug(0);
29 AliHBTParticle::SetDebug(0);
30 AliLoader::SetDebug(0);
33 AliHBTParticleCut c2 = c1;
36 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
38 //datatype defines type of data to be read
39 // Kine - analyzes Kine Tree: simulated particles
40 // TPC - analyzes TPC tracking + particles corresponding to these tracks
41 // ITSv1 - analyzes ITSv1 ----------------------//--------------------------
42 // ITSv2 - analyzes ITSv2 ----------------------//--------------------------
44 //processopt defines option passed to AliHBTAnlysis::Process method
45 // default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them
46 // Tracks - process only recontructed tracks
47 // Particles - process only simulated particles
49 //Reads data from diroctories from first to last(including)
50 // For examples if first=3 and last=5 it reads from
54 //if first or last is negative (or both), it reads from current directory
56 //these names I use when analysis is done directly from CASTOR files via RFIO
57 // const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
58 // const char* serie="standard";
59 // const char* field = "0.4";
61 const char* basedir=".";
63 const char* field = "";
65 Bool_t threeDcuts = kFALSE;
67 // Dynamically link some shared libs
68 // gROOT->LoadMacro("loadlibs.C");
71 gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
74 /***********************************************************/
75 //Create analysis and reader
76 AliHBTAnalysis * analysis = new AliHBTAnalysis();
77 analysis->SetCutsOnTracks();
80 Int_t kine = strcmp(datatype,"Kine");
81 Int_t ESD = strcmp(datatype,"ESD");
82 Int_t TPC = strcmp(datatype,"TPC");
83 Int_t ITSv1 = strcmp(datatype,"ITSv1");
84 Int_t ITSv2 = strcmp(datatype,"ITSv2");
85 Int_t intern = strcmp(datatype,"Intern");
89 reader = new AliHBTReaderKineTree();
90 processopt="Particles"; //this reader by definition reads only simulated particles
94 AliHBTReaderESD* esdreader = new AliHBTReaderESD();
95 esdreader->ReadParticles(kTRUE);
96 esdreader->SetNumberOfTrackPoints(5,30.);
97 esdreader->SetClusterMap();
102 reader = new AliHBTReaderTPC();
103 //((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
104 ((AliHBTReaderTPC*)reader)->SetClusterMap();
108 reader = new AliHBTReaderITSv1();
112 reader = new AliHBTReaderITSv2();
116 reader = new AliHBTReaderInternal("data.root");
120 cerr<<"Option "<<datatype<<" not recognized. Exiting"<<endl;
125 if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
126 {//read from many dirs dirs
128 dirs = new TObjArray(last-first+1);
129 for (Int_t i = first; i<=last; i++)
131 sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
132 TObjString *odir= new TObjString(buff);
136 reader->SetDirs(dirs);
138 /***********************************************************/
139 /***** R E A D E R ***************************************/
140 /***********************************************************/
142 //we want have only low pt pi+ so set a cut to reader
143 AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
144 readerpartcut->SetPtRange(0.0,10.0);
145 readerpartcut->SetPID(kPiPlus);
146 reader->AddParticleCut(readerpartcut);//read this particle type with this cut
148 // readerpartcut->SetPtRange(0.0,10.0);
149 // readerpartcut->SetPID(kPiMinus);
150 // reader->AddParticleCut(readerpartcut);//read this particle type with this cut
152 analysis->SetReader(reader);
154 AliHBTPairCut *paircut = new AliHBTPairCut();
155 Float_t qinvmin = 0.0;
156 Float_t qinvmax = 0.05;//50MeV
157 paircut->SetQInvRange(qinvmin,qinvmax);
159 // paircut->SetAvSeparationRange(10.); //AntiMerging
160 paircut->SetClusterOverlapRange(-.5,1.0);//AntiSplitting
162 // AliHBTParticleCut* partcut= new AliHBTParticleCut();
163 // partcut->SetPID(kPiPlus);
164 // paircut->SetFirstPartCut(partcut);
165 // partcut->SetPID(kPiMinus);
166 // paircut->SetSecondPartCut(partcut);
168 analysis->SetGlobalPairCut(paircut);
170 /***********************************************************/
171 /***** W E I G H T S ******************************/
172 /***********************************************************/
174 AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
175 AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
176 AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
178 AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
179 //AliHBTLLWeights::Instance()->SetR1dw(8.0);
181 AliHBTLLWeights::Instance()->Init();
182 AliHBTLLWeights::Instance()->Set();
185 AliHBTWeightTheorQInvFctn *wqinvcfP = new AliHBTWeightTheorQInvFctn(100,qinvmax);
186 wqinvcfP->SetNumberOfBinsToScale(30);
187 wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function");
188 analysis->AddParticleFunction(wqinvcfP);
190 /************************************************************/
191 /**** Q INV Correlation Function ************************/
192 /************************************************************/
195 AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
196 AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
198 analysis->AddTrackFunction(qinvcfT);
199 analysis->AddParticleFunction(qinvcfP);
202 qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
203 qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
205 /************************************************************/
206 /**** Q OUT Correlation Function ************************/
207 /************************************************************/
209 AliHBTQOutLCMSCorrelFctn* qoutP = new AliHBTQOutLCMSCorrelFctn();
210 qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
211 AliHBTQOutLCMSCorrelFctn* qoutT = new AliHBTQOutLCMSCorrelFctn();
212 qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
214 AliHBTPairCut *outPairCut = new AliHBTPairCut();
215 outPairCut->SetQOutCMSLRange(0.0,0.15);
216 outPairCut->SetQSideCMSLRange(0.0,0.02);
217 outPairCut->SetQLongCMSLRange(0.0,0.02);
218 qoutP->SetPairCut(outPairCut);
219 qoutT->SetPairCut(outPairCut);
221 //analysis->AddTrackFunction(qoutP);
222 //analysis->AddParticleFunction(qoutT);
224 /************************************************************/
225 /**** Q SIDE Correlation Function ***********************/
226 /************************************************************/
228 AliHBTQSideLCMSCorrelFctn* qsideP = new AliHBTQSideLCMSCorrelFctn(100,qinvmax);
229 qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
230 AliHBTQSideLCMSCorrelFctn* qsideT = new AliHBTQSideLCMSCorrelFctn(100,qinvmax);
231 qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
233 AliHBTPairCut *sidePairCut = new AliHBTPairCut();
234 sidePairCut->SetQOutCMSLRange(0.0,0.02);
235 sidePairCut->SetQSideCMSLRange(0.0,0.15);
236 sidePairCut->SetQLongCMSLRange(0.0,0.02);
237 qsideP->SetPairCut(sidePairCut);
238 qsideT->SetPairCut(sidePairCut);
240 //analysis->AddTrackFunction(qsideP);
241 //analysis->AddParticleFunction(qsideT);
243 /************************************************************/
244 /**** Q LONG Correlation Function ***********************/
245 /************************************************************/
247 AliHBTQLongLCMSCorrelFctn* qlongP = new AliHBTQLongLCMSCorrelFctn(100,qinvmax);
248 qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
249 AliHBTQLongLCMSCorrelFctn* qlongT = new AliHBTQLongLCMSCorrelFctn(100,qinvmax);
250 qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
252 AliHBTPairCut *longPairCut = new AliHBTPairCut();
253 longPairCut->SetQOutCMSLRange(0.0,0.02);
254 longPairCut->SetQSideCMSLRange(0.0,0.02);
255 longPairCut->SetQLongCMSLRange(0.0,0.15);
256 qlongP->SetPairCut(longPairCut);
257 qlongT->SetPairCut(longPairCut);
259 //analysis->AddTrackFunction(qlongT);
260 //analysis->AddParticleFunction(qlongT);
263 /************************************************************/
264 /************************************************************/
265 /************************************************************/
266 /**** R E S O L U T I O N S ***********************/
267 /************************************************************/
268 /************************************************************/
269 /************************************************************/
273 AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015); //QInvLCMS Res Vs Kt
274 if (threeDcuts) qinvVsktF.SetPairCut(paircut);
275 //analysis->AddResolutionFunction(&qinvVsktF);
277 AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05); //QOutLCMS Res Vs Kt
278 if (threeDcuts) qoutVsktF.SetPairCut(outPairCut);
279 //analysis->AddResolutionFunction(&qoutVsktF);
281 AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015); //QSideLCMS Res Vs Kt
282 if (threeDcuts) qsideVsktF.SetPairCut(sidePairCut);
283 //analysis->AddResolutionFunction(&qsideVsktF);
285 AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015); //QLongLCMS Res Vs Kt
286 if (threeDcuts) qlongVsktF.SetPairCut(longPairCut);
287 //analysis->AddResolutionFunction(&qlongVsktF);
289 AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvLCMS Res Vs QInvLCMS
290 //analysis->AddResolutionFunction(&qInvVsqinvF);
292 AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05); //QOutLCMS Res Vs QInvLCMS
293 if (threeDcuts) qoutVsqinvF.SetPairCut(outPairCut);
294 //analysis->AddResolutionFunction(&qoutVsqinvF);
296 AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideLCMS Res Vs QInvLCMS
297 if (threeDcuts) qsideVsqinvF.SetPairCut(sidePairCut);
298 //analysis->AddResolutionFunction(&qsideVsqinvF);
300 AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongLCMS Res Vs QInvLCMS
301 if (threeDcuts) qlongVsqinvF.SetPairCut(longPairCut);
302 //analysis->AddResolutionFunction(&qlongVsqinvF);
305 AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
306 if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
307 //analysis->AddResolutionFunction(&pairThetaVsqinv);
309 AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
310 if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
311 //analysis->AddResolutionFunction(&pairThetaVsKt);
314 phipc.SetQLongCMSLRange(0.0,0.02);
316 AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
317 if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
318 //analysis->AddResolutionFunction(&pairPhiVsqinv);
320 AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
321 if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
322 //analysis->AddResolutionFunction(&pairPhiVsKt);
325 AliHBTQOutResolVsQOutFctn qoutVsqoutF; //QOutLCMS Res Vs QOut
326 if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
327 //analysis->AddResolutionFunction(&qoutVsqoutF);
329 AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideLCMS Res Vs QSide
330 if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
331 //analysis->AddResolutionFunction(&qsideVsqsideF);
334 AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongLCMS Res Vs QLong
335 if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
336 //analysis->AddResolutionFunction(&qlongVsqlongF);
339 /************************************************************/
340 /************************************************************/
341 /************************************************************/
342 /**** D I S T R B U T I O N S ***********************/
343 /************************************************************/
344 /************************************************************/
345 /************************************************************/
347 /************************************************************/
348 /**** OUT VS KT ***********************/
349 /************************************************************/
350 AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
351 qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
352 AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
353 qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
354 if (threeDcuts) qoutVsKtDistP.SetPairCut(outPairCut);
355 if (threeDcuts) qoutVsKtDistT.SetPairCut(outPairCut);
356 //analysis->AddParticleFunction(&qoutVsKtDistP);
357 //analysis->AddTrackFunction(&qoutVsKtDistT);
358 /************************************************************/
359 /**** Side VS KT ***********************/
360 /************************************************************/
361 AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
362 qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
363 AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
364 qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
365 if (threeDcuts) qSideVsKtDistP.SetPairCut(sidePairCut);
366 if (threeDcuts) qSideVsKtDistT.SetPairCut(sidePairCut);
367 //analysis->AddParticleFunction(&qSideVsKtDistP);
368 //analysis->AddTrackFunction(&qSideVsKtDistT);
369 /************************************************************/
370 /**** Long VS KT ***********************/
371 /************************************************************/
372 AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
373 qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
374 AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
375 qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
376 if (threeDcuts) qLongVsKtDistP.SetPairCut(longPairCut);
377 if (threeDcuts) qLongVsKtDistT.SetPairCut(longPairCut);
378 //analysis->AddParticleFunction(&qLongVsKtDistP);
379 //analysis->AddTrackFunction(&qLongVsKtDistT);
381 /************************************************************/
382 /************************************************************/
383 /************************************************************/
384 /**** M O N I T O R R E S O L U T I O N S *************/
385 /************************************************************/
386 /************************************************************/
387 /************************************************************/
389 AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
390 //analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
392 AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
393 //analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
395 AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
396 //analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);
398 AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
399 //analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);
401 AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
402 analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);
404 AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
405 analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);
407 AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
408 analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
413 /************************************************************/
414 /**** P R O C E S S ***********************/
415 /************************************************************/
416 analysis->SetBufferSize(2);
417 analysis->Process(processopt);
419 TFile histoOutput(outfile,"recreate");
420 analysis->WriteFunctions();
426 delete readerpartcut;