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