Support files added
[u/mrichter/AliRoot.git] / MFT / MergeClustersMFT.C
1 #include "AliMFTCluster.h"
2 #include "TFile.h"
3 #include "AliMFTCluster.h"
4 #include "TObjArray.h"
5 #include "TClonesArray.h"
6 #include "TTree.h"
7 #include "AliMFTConstants.h"
8 #include "TBranch.h"
9 #include "TRandom.h"
10 #include "TMath.h"
11
12 TFile *fileInPileUp=0, *fileInSignal, *fileInUnderlying, *fileOutPileUp=0, *fileOutTotal;
13
14 enum {kMergePileUpClusters_pp, kMergePileUpClusters_PbPb, kMergeAllClusters_pp, kMergeAllClusters_PbPb};
15 enum {k_pp, k_PbPb};
16
17 const Char_t *fPileUpDir[2] = {"$UnderlyingDirPP/spreadVtx",      // local directory with pile-up event p-p
18                                "$UnderlyingDirPbPb/spreadVtx"};   // local directory with pile-up event Pb-Pb
19
20 const Int_t fNRunsPileUp_pp   = 100;
21 const Int_t fNRunsPileUp_PbPb =   3;
22
23 const Int_t fNRunsPileUpAvailable[2] = {500, 500};    // pp, PbPb
24
25 TClonesArray *fRecPointsPerPlaneIn[AliMFTConstants::fNMaxPlanes] = {0};
26 TClonesArray *fRecPointsPerPlaneOut[AliMFTConstants::fNMaxPlanes] = {0};
27
28 void MergePileUpClusters(Int_t *runs, Int_t nRunsPileUp, const Char_t *pileUpDir, Int_t nEvents);
29 void MergeAllClusters();
30                                
31 //====================================================================================================================================================
32
33 void MergeClustersMFT(Int_t option = 9999,
34                       Int_t nEvents = 9999,
35                       Int_t maxPileUp = 9999) {
36   
37   if (option==kMergePileUpClusters_pp) {
38
39     Int_t runs_pp[fNRunsPileUp_pp] = {0};
40     Int_t nFiles_pp = 0;
41     while (nFiles_pp<TMath::Min(fNRunsPileUp_pp,maxPileUp)) {
42       Int_t run = gRandom->Integer(fNRunsPileUpAvailable[k_pp]);
43       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",fPileUpDir[k_pp],run));
44       if (!fileInPileUp->IsOpen() || !fileInPileUp->cd("Event0")) {
45         delete fileInPileUp;
46         continue;
47       }
48       Bool_t runAlreadyLoaded=kFALSE;
49       for (Int_t iAddedRun=0; iAddedRun<nFiles_pp; iAddedRun++) {
50         if (run==runs_pp[iAddedRun]) {
51           runAlreadyLoaded=kTRUE;
52           break;
53         }
54       }
55       if (!runAlreadyLoaded) runs_pp[nFiles_pp++] = run;
56       delete fileInPileUp;
57     }
58     MergePileUpClusters(runs_pp, TMath::Min(fNRunsPileUp_pp,maxPileUp), fPileUpDir[k_pp], nEvents);
59   }
60     
61   else if (option==kMergePileUpClusters_PbPb) {
62     Int_t runs_PbPb[fNRunsPileUp_PbPb] = {0};
63     Int_t nFiles_PbPb = 0;
64     while (nFiles_PbPb<TMath::Min(fNRunsPileUp_PbPb,maxPileUp)) {
65       Int_t run = gRandom->Integer(fNRunsPileUpAvailable[k_PbPb]);
66       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",fPileUpDir[k_PbPb],run));
67       if (!fileInPileUp->IsOpen() || !fileInPileUp->cd("Event0")) {
68         delete fileInPileUp;
69         continue;
70       }
71       Bool_t runAlreadyLoaded=kFALSE;
72       for (Int_t iAddedRun=0; iAddedRun<nFiles_PbPb; iAddedRun++) {
73         if (run==runs_PbPb[iAddedRun]) {
74           runAlreadyLoaded=kTRUE;
75           break;
76         }
77       }
78       if (!runAlreadyLoaded) runs_PbPb[nFiles_PbPb++] = run;
79       delete fileInPileUp;
80     }
81     MergePileUpClusters(runs_PbPb, TMath::Min(fNRunsPileUp_PbPb,maxPileUp), fPileUpDir[k_PbPb], nEvents);
82   }
83
84   if (option>=kMergeAllClusters_pp) MergeAllClusters();
85
86 }
87
88 //====================================================================================================================================================
89
90 void MergePileUpClusters(Int_t *runs, Int_t nRunsPileUp, const Char_t *pileUpDir, Int_t nEvents) {
91
92   fileOutPileUp = new TFile("./MFT.RecPoints.PileUp.root", "recreate");
93     
94   Bool_t eventExist = kTRUE;
95   Int_t iEv=0;
96
97   while (eventExist && iEv<nEvents) {
98
99     for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
100       fRecPointsPerPlaneOut[iPlane] = new TClonesArray("AliMFTCluster");
101     }
102
103     if (iEv) fileOutPileUp = new TFile("./MFT.RecPoints.PileUp.root", "update");
104
105     printf("Merging Event %d\n",iEv);
106     
107     fileOutPileUp-> cd();
108     TTree *treeOutPileUp = new TTree("TreeR_OutPileUp", "Reconstructed Points Container");
109
110     Int_t nFilesAdded = 0;
111     while (nFilesAdded<nRunsPileUp) {
112       fileInPileUp = new TFile(Form("%s/run_%d/MFT.RecPoints.MCShifted.root",pileUpDir,runs[nFilesAdded++]));
113       printf("Merging file %s\n", fileInPileUp->GetName());
114       if (!fileInPileUp->cd(Form("Event%d",iEv))) {
115         eventExist = kFALSE;
116         break;
117       }
118       for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
119         fRecPointsPerPlaneIn[iPlane] = new TClonesArray("AliMFTCluster");
120       }
121       TTree *treeInPileUp = (TTree*) gDirectory->Get("TreeR");
122       treeInPileUp -> SetName("TreeR_InPileUp");
123       Int_t iPlane=0;
124       while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
125         //      printf("Plane %02d\n",iPlane);
126         treeInPileUp ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
127         TBranch *branch = treeOutPileUp->GetBranch(Form("Plane_%02d",iPlane));
128         if (!branch) treeOutPileUp->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane])); 
129         iPlane++;
130       }
131       iPlane=0;
132       treeInPileUp -> GetEntry(0);
133       while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
134         Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
135         for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
136           //      printf("Cluster %4d\n",iCluster);
137           AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
138           new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
139         }
140         iPlane++;
141       }
142       for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
143         delete fRecPointsPerPlaneIn[jPlane];
144       }
145       delete treeInPileUp;
146       fileInPileUp -> Close();
147       delete fileInPileUp;
148     }
149     if (eventExist) {
150       treeOutPileUp -> Fill();
151       fileOutPileUp -> mkdir(Form("Event%d",iEv));
152       fileOutPileUp -> cd(Form("Event%d",iEv));
153       treeOutPileUp -> SetName("TreeR");
154       treeOutPileUp -> Write();
155       delete treeOutPileUp;
156       fileOutPileUp -> Close();
157       delete fileOutPileUp;
158       for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
159         delete fRecPointsPerPlaneOut[iPlane];
160       }
161       iEv++;
162     }
163   }  
164
165 }
166
167 //====================================================================================================================================================
168
169 void MergeAllClusters() {
170
171   fileOutTotal = new TFile("./MFT.RecPoints.root", "recreate");
172     
173   Bool_t eventExist = kTRUE;
174   Int_t iEv=0;
175
176   while (eventExist) {
177
178     for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
179       fRecPointsPerPlaneOut[iPlane] = new TClonesArray("AliMFTCluster");
180     }
181
182     if (iEv) fileOutTotal = new TFile("./MFT.RecPoints.root", "update");
183
184     printf("Merging Event %d\n",iEv);
185     
186     fileOutTotal-> cd();
187     TTree *treeOutTotal = new TTree("TreeR_Out", "Reconstructed Points Container");
188
189     // Adding Pile-Up Clusters
190
191     fileInPileUp = new TFile("./MFT.RecPoints.PileUp.root");
192     printf("Merging file %s\n", fileInPileUp->GetName());
193     if (!fileInPileUp->cd(Form("Event%d",iEv))) {
194       eventExist = kFALSE;
195       break;
196     }
197     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
198       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
199     }
200     TTree *treeInPileUp = (TTree*) gDirectory->Get("TreeR");
201     treeInPileUp -> SetName("TreeR_InPileUp");
202     Int_t iPlane=0;
203     while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
204       //        printf("Plane %02d\n",iPlane);
205       treeInPileUp ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
206       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
207       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane])); 
208       iPlane++;
209     }
210     iPlane=0;
211     treeInPileUp -> GetEntry(0);
212     while (treeInPileUp->GetBranch(Form("Plane_%02d",iPlane))) {
213       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
214       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
215         //        printf("Cluster %4d\n",iCluster);
216         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
217         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
218       }
219       iPlane++;
220     }
221     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
222       delete fRecPointsPerPlaneIn[jPlane];
223     }
224     delete treeInPileUp;
225     fileInPileUp -> Close();
226     delete fileInPileUp;
227
228     // Adding Underlying Event Clusters
229
230     fileInUnderlying = new TFile("./MFT.RecPoints.Underlying.root");
231     printf("Merging file %s\n", fileInUnderlying->GetName());
232     if (!fileInUnderlying->cd(Form("Event%d",iEv))) {
233       eventExist = kFALSE;
234       break;
235     }
236     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
237       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
238     }
239     TTree *treeInUnderlying = (TTree*) gDirectory->Get("TreeR");
240     treeInUnderlying -> SetName("TreeR_InUnderlying");
241     iPlane=0;
242     while (treeInUnderlying->GetBranch(Form("Plane_%02d",iPlane))) {
243       //        printf("Plane %02d\n",iPlane);
244       treeInUnderlying ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
245       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
246       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane])); 
247       iPlane++;
248     }
249     iPlane=0;
250     treeInUnderlying -> GetEntry(0);
251     while (treeInUnderlying->GetBranch(Form("Plane_%02d",iPlane))) {
252       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
253       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
254         //        printf("Cluster %4d\n",iCluster);
255         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
256         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
257       }
258       iPlane++;
259     }
260     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
261       delete fRecPointsPerPlaneIn[jPlane];
262     }
263     delete treeInUnderlying;
264     fileInUnderlying -> Close();
265     delete fileInUnderlying;
266
267     // Adding Signal Clusters
268
269     fileInSignal = new TFile("./MFT.RecPoints.Signal.root");
270     printf("Merging file %s\n", fileInSignal->GetName());
271     if (!fileInSignal->cd(Form("Event%d",iEv))) {
272       eventExist = kFALSE;
273       break;
274     }
275     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
276       fRecPointsPerPlaneIn[jPlane] = new TClonesArray("AliMFTCluster");
277     }
278     TTree *treeInSignal = (TTree*) gDirectory->Get("TreeR");
279     treeInSignal -> SetName("TreeR_InSignal");
280     iPlane=0;
281     while (treeInSignal->GetBranch(Form("Plane_%02d",iPlane))) {
282       //        printf("Plane %02d\n",iPlane);
283       treeInSignal ->SetBranchAddress(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneIn[iPlane]));
284       TBranch *branch = treeOutTotal->GetBranch(Form("Plane_%02d",iPlane));
285       if (!branch) treeOutTotal->Branch(Form("Plane_%02d",iPlane), &(fRecPointsPerPlaneOut[iPlane])); 
286       iPlane++;
287     }
288     iPlane=0;
289     treeInSignal -> GetEntry(0);
290     while (treeInSignal->GetBranch(Form("Plane_%02d",iPlane))) {
291       Int_t nClusters = fRecPointsPerPlaneIn[iPlane]->GetEntries();
292       for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
293         //        printf("Cluster %4d\n",iCluster);
294         AliMFTCluster *newCluster = (AliMFTCluster*) fRecPointsPerPlaneIn[iPlane]->At(iCluster);
295         new ((*fRecPointsPerPlaneOut[iPlane])[fRecPointsPerPlaneOut[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
296       }
297       iPlane++;
298     }
299     for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
300       delete fRecPointsPerPlaneIn[jPlane];
301     }
302     delete treeInSignal;
303     fileInSignal -> Close();
304     delete fileInSignal;
305
306     if (eventExist) {
307       treeOutTotal -> Fill();
308       fileOutTotal -> mkdir(Form("Event%d",iEv));
309       fileOutTotal -> cd(Form("Event%d",iEv));
310       treeOutTotal -> SetName("TreeR");
311       treeOutTotal -> Write();
312       delete treeOutTotal;
313       fileOutTotal -> Close();
314       delete fileOutTotal;
315       for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) {
316         delete fRecPointsPerPlaneOut[jPlane];
317       }
318       iEv++;
319     }
320
321   }   
322
323 }
324
325 //====================================================================================================================================================
326
327 //  LocalWords:  nFiles