(martin) pt vs eta correction matrix calculation macro using CorrectionMatrix2D class.
[u/mrichter/AliRoot.git] / PWG0 / ptSpectra / makeCorrectionPtEta.C
1 //
2 // Script to calculate PtvsEta correction map using the CorrectionMatrix2D class.
3 // 
4
5 makeCorrectionPtEta(Char_t* dataDir, Int_t nRuns=20) {
6
7   Char_t str[256];
8
9   //gSystem->Load("/home/poghos/CERN2006/pp/esdTrackCuts/libESDtrackQuality.so");
10   gSystem->Load("libPWG0base");
11   gSystem->Load("libCorrectionMatrix2D.so");
12   // ########################################################
13   // selection of esd tracks
14   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();    
15   esdTrackCuts->DefineHistograms(1);
16   esdTrackCuts->SetMinNClustersTPC(50);
17   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
18   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
19   esdTrackCuts->SetRequireTPCRefit(kTRUE);
20   esdTrackCuts->SetMinNsigmaToVertex(3);
21   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
22   AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
23   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
24
25   // ########################################################
26   // definition of PtEta correction object
27   CorrectionMatrix2D* PtEtaMap = new CorrectionMatrix2D("CorrectionMatrix2");
28   PtEtaMap->SetHist("PtEta",80,0.,10.,120,-2.,2.);
29   //Float_t x[]={-20., -15., -5., 0., 6., 20.};
30   //Float_t y[]={-6. , -3.,  -1., 2., 3., 6. };
31   //PtEtaMap->SetHist("PtEta",5,x,5,y);
32   PtEtaMap->SetHistTitle("p_{t}","#eta");
33
34
35   // ########################################################
36   // get the data dir  
37   Char_t execDir[256];
38   sprintf(execDir,gSystem->pwd());
39   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
40   TList* dirList            = baseDir->GetListOfFiles();
41   Int_t nDirs               = dirList->GetEntries();
42   // go back to the dir where this script is executed
43   gSystem->cd(execDir);
44
45   // ########################################################
46   // definition of used pointers
47   TFile* esdFile;
48   TTree* esdTree;
49   TBranch* esdBranch;
50
51   AliESD* esd =0;
52   // ########################################################
53   // loop over runs
54   Int_t nRunCounter = 0;
55   for (Int_t r=0; r<nDirs; r++) {
56
57     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
58     if (!presentDir || !presentDir->IsDirectory())
59       continue;
60     // check that the files are there
61     TString currentDataDir;
62     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
63     cout << "Processing directory " << currentDataDir.Data() << endl;
64     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
65           (!gSystem->Which(currentDataDir,"AliESDs.root")))
66       continue;
67
68     if (nRunCounter++ >= nRuns)
69       break;
70
71     cout << "run #" << nRunCounter << endl;
72
73     // #########################################################
74     // setup galice and runloader
75     if (gAlice) {
76       delete gAlice->GetRunLoader();
77       delete gAlice;
78       gAlice=0;
79     }
80
81     sprintf(str,"%s/galice.root",currentDataDir.Data());
82     AliRunLoader* runLoader = AliRunLoader::Open(str);
83
84     runLoader->LoadgAlice();
85     gAlice = runLoader->GetAliRun();
86     runLoader->LoadKinematics();
87     runLoader->LoadHeader();
88
89    
90     // #########################################################
91     // open esd file and get the tree
92
93     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
94     // close it first to avoid memory leak
95     if (esdFile)
96       if (esdFile->IsOpen())
97         esdFile->Close();
98
99     esdFile = TFile::Open(str);
100     esdTree = (TTree*)esdFile->Get("esdTree");
101     if (!esdTree)
102       continue;
103     esdBranch = esdTree->GetBranch("ESD");
104     esdBranch->SetAddress(&esd);
105     if (!esdBranch)
106       continue;
107
108     // ########################################################
109     // Magnetic field
110     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
111     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
112
113     // ########################################################
114     // getting number of events
115     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
116     Int_t nESDEvents = esdBranch->GetEntries();
117     
118     if (nEvents!=nESDEvents) {
119       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
120       return;
121     }
122     // ########################################################
123     // loop over number of events
124     cout << " looping over events..." << endl;
125     for(Int_t i=0; i<nEvents; i++) {
126
127       esdBranch->GetEntry(i);
128       runLoader->GetEvent(i);              
129
130       // ########################################################
131       // get the EDS vertex
132       AliESDVertex* vtxESD = esd->GetVertex();
133
134       Double_t vtx[3];
135       Double_t vtx_res[3];
136       vtxESD->GetXYZ(vtx);          
137     
138       vtx_res[0] = vtxESD->GetXRes();
139       vtx_res[1] = vtxESD->GetYRes();
140       vtx_res[2] = vtxESD->GetZRes();
141
142       // the vertex should be reconstructed
143       if (strcmp(vtxESD->GetName(),"default")==0) 
144         continue;
145
146       // the resolution should be reasonable???
147       if (vtx_res[2]==0 || vtx_res[2]>0.1)
148         continue;
149
150       // ########################################################
151       // loop over mc particles
152       AliStack* particleStack = runLoader->Stack();
153       Int_t nPrim    = particleStack->GetNprimary();
154
155       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
156       
157         TParticle* part = particleStack->Particle(i_mc);      
158         if (!part || strcmp(part->GetName(),"XXX")==0) 
159           continue;
160       
161         TParticlePDG* pdgPart = part->GetPDG();
162
163         Bool_t prim = kFALSE;
164         // check if it's a primary particle - is there a better way ???
165         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
166           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
167             prim = kTRUE;
168         }
169         if (!prim)
170           continue;
171
172         if (pdgPart->Charge()==0)
173           continue;
174
175         if (prim)
176           PtEtaMap->FillGene(part->Pt(), part->Eta());  
177         
178       }// end of mc particle
179
180       // ########################################################
181       // loop over esd tracks      
182       Int_t nTracks = esd->GetNumberOfTracks();            
183       for (Int_t t=0; t<nTracks; t++) {
184         AliESDtrack* esdTrack = esd->GetTrack(t);      
185         
186         // cut the esd track?
187         if (!esdTrackCuts->AcceptTrack(esdTrack))
188           continue;
189         
190       Double_t Ptr0[3];
191       if(!esdTrack->GetPxPyPz(Ptr0)) continue;
192       fTrP= new TVector3(Ptr0);
193         
194         PtEtaMap->FillMeas(fTrP->Pt(), fTrP->Eta());    
195
196       } // end of track loop
197     } // end  of event loop
198   } // end of run loop
199
200   PtEtaMap->Finish();  
201
202   TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
203   
204   esdTrackCuts->SaveHistograms("EsdTrackCuts");
205   PtEtaMap->SaveHistograms();
206   
207   fout->Write();
208   fout->Close();
209
210 }