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