2 // Script to test the dN/dEta analysis using the dNdEtaAnalysis and
3 // dNdEtaCorrection classes. Note that there is a cut on the events,
4 // so the measurement will be biassed.
7 testAnalysis(Char_t* dataDir, Int_t nRuns=20) {
11 gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
12 gSystem->Load("libdNdEta.so");
14 // ########################################################
15 // selection of esd tracks
16 ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();
17 esdTrackCuts->DefineHistograms(1);
19 esdTrackCuts->SetMinNClustersTPC(50);
20 esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
21 esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
22 esdTrackCuts->SetRequireTPCRefit(kTRUE);
24 esdTrackCuts->SetMinNsigmaToVertex(3);
25 esdTrackCuts->SetAcceptKingDaughters(kFALSE);
27 AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
29 // ########################################################
30 // definition of dNdEta objects
32 dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
34 dNdEtaMap->LoadHistograms("correction_map.root","dndeta_correction");
35 dNdEtaMap->RemoveEdges(2,0,2);
37 dNdEtaAnalysis* analyse = new dNdEtaAnalysis("dndeta");
39 // ########################################################
42 sprintf(execDir,gSystem->pwd());
43 TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
44 TList* dirList = baseDir->GetListOfFiles();
45 Int_t nDirs = dirList->GetEntries();
46 // go back to the dir where this script is executed
50 // ########################################################
51 // definition of used pointers
58 // ########################################################
60 Int_t nRunCounter = 0;
61 for (Int_t r=1; r<=nDirs; r++) {
63 TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
64 if (!presentDir->IsDirectory())
66 // check that the files are there
67 sprintf(str,"%s/%s",dataDir, presentDir->GetName());
68 if ((!gSystem->Which(str,"galice.root")) ||
69 (!gSystem->Which(str,"AliESDs.root")))
72 if (nRunCounter++ >= nRuns)
75 cout << "run #" << nRunCounter << endl;
77 // #########################################################
78 // setup galice and runloader
80 delete gAlice->GetRunLoader();
85 sprintf(str,"%s/run%d/galice.root",dataDir,r);
86 AliRunLoader* runLoader = AliRunLoader::Open(str);
88 runLoader->LoadgAlice();
89 gAlice = runLoader->GetAliRun();
90 runLoader->LoadKinematics();
91 runLoader->LoadHeader();
93 // #########################################################
94 // open esd file and get the tree
96 sprintf(str,"%s/run%d/AliESDs.root",dataDir,r);
97 // close it first to avoid memory leak
99 if (esdFile->IsOpen())
102 esdFile = TFile::Open(str);
103 esdTree = (TTree*)esdFile->Get("esdTree");
106 esdBranch = esdTree->GetBranch("ESD");
107 esdBranch->SetAddress(&esd);
111 // ########################################################
113 AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
114 //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
116 // ########################################################
117 // getting number of events
118 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
119 Int_t nESDEvents = esdBranch->GetEntries();
121 if (nEvents!=nESDEvents) {
122 cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
125 // ########################################################
126 // loop over number of events
127 cout << " looping over events..." << endl;
128 for(Int_t i=1; i<nEvents; i++) {
130 esdBranch->GetEntry(i);
131 runLoader->GetEvent(i);
134 // ########################################################
135 // get the EDS vertex
136 AliESDVertex* vtxESD = esd->GetVertex();
142 vtx_res[0] = vtxESD->GetXRes();
143 vtx_res[1] = vtxESD->GetYRes();
144 vtx_res[2] = vtxESD->GetZRes();
146 // we need a good vertex
147 // => there will be a bias on the measurement, since this cuts away some events
148 if (strcmp(vtxESD->GetName(),"default")==0)
150 if (vtx_res[2]==0 || vtx_res[2]>0.1)
153 // ########################################################
154 // loop over esd tracks
155 Int_t nTracks = esd->GetNumberOfTracks();
156 for (Int_t t=0; t<nTracks; t++) {
157 AliESDtrack* esdTrack = esd->GetTrack(t);
159 if (!esdTrackCuts->AcceptTrack(esdTrack))
162 AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
164 if (tpcTrack->GetAlpha()==0) {
165 cout << " Warning esd track alpha = 0" << endl;
169 Float_t theta = tpcTrack->Theta();
170 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
171 Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta);
173 dNdEtaMap->FillMeas(vtx[2], eta);
175 analyse ->FillTrack(vtx[2], eta, correction);
177 } // end of track loop
178 analyse->FillEvent(vtx[2]);
180 } // end of event loop
185 TFile* fout = new TFile("out.root","RECREATE");
187 esdTrackCuts->SaveHistograms("esd_tracks_cuts");
188 dNdEtaMap->SaveHistograms();
189 analyse ->SaveHistograms();