[u/mrichter/AliRoot.git] / HLT / misc / MakePileup.cxx

#ifndef __CINT__
#include "AliL3FileHandler.h"
#include "AliL3DigitData.h"
#include "AliL3Transform.h"
#include "MakePileup.h"
#include <stdio.h>
#include <iostream.h>

void QSort(AliL3DigitData **a,Int_t first,Int_t last);
Int_t CompareDigits(AliL3DigitData *a,AliL3DigitData *b);
#endif

void MakePileup(Int_t triggerevent)
{
  Int_t srow[2] = {0,175};
  UInt_t size;
  const Int_t NEvents = 25;
  //Int_t Offset[25] = {0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240};
  //Int_t Offset[25] = {0,20,40,60,80,100,120,140,160,180,200,220,240,-20,-40,-60,-80,-100,-120,-140,-160,-180,-200,-220,-240};
  Int_t Offset[25]   = {30,60,90,120,150,180,210,240,270,300,330,360,-30,-60,-90,-120,-150,-180,-210,-240,-270,-300,-330,-360,-390};
  Offset[triggerevent]=0;
  AliL3FileHandler *hand[NEvents];
  AliL3DigitRowData *data[NEvents];
  AliL3DigitData *rowData[NEvents];
  AliL3DigitData *test;
  Char_t Carry[255];
  //Int_t event = 0;
  Int_t DigitsTot = 0;
  Int_t MeVDigit;
  Int_t DigitsPerRow = 0;
  //AliL3DigitRowData *NData;
  Byte_t *NData;
  Int_t tot_dig=0;
  
  AliL3Transform *transform = new AliL3Transform();
  const Int_t maxdigits=50000;
  AliL3DigitData *temp[maxdigits];
  AliL3DigitData *temp2[maxdigits];
  
  
  for(Int_t slice=0 ; slice<36 ; slice++){
    for(Int_t event = 0; event < NEvents ; event++){
      hand[event] = new AliL3FileHandler();
      sprintf(Carry,"/prog/alice/data/Rawdata/1_patch/pp/event_%d/digits_%d_0.raw",event,slice);
      if(!hand[event]->SetBinaryInput(Carry))
	cerr<<" Error opening file :"<<Carry<<endl;
      cout << "Open File: " << Carry;
      hand[event]->Init(slice,0,srow);
      data[event] = hand[event]->CompBinary2Memory(size);
      cout<<" done"<<endl;
      hand[event]->CloseBinaryInput();
    }
    DigitsTot=0;
    Int_t sign = slice < 18 ? -1 : 1;
    for(Int_t row = 0 ; row < 176 ; row++){
      
      for(Int_t event = 0; event < NEvents ; event++){
	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
	for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
	  {
	    Int_t time = rowData[event][dig].fTime + sign*Offset[event];
	    if(time < transform->GetNTimeBins() && time >= 0)
	      DigitsTot++;
	  }
	//DigitsTot += data[event]->fNDigit;
	hand[event]->UpdateRowPointer(data[event]);
      }
      
    }
    //cout << "Try to allocate : " << DigitsTot << " Digits in row : " << row << endl;
    
    
    //NData = (AliL3DigitRowData*) malloc(AllDigitsSize);// Allocate !
    
    Int_t AllDigitsSize = sizeof(AliL3DigitData) * DigitsTot + sizeof(AliL3DigitRowData) * 176;
    NData = new Byte_t[AllDigitsSize];
    memset(NData,0,AllDigitsSize);
    AliL3DigitRowData *AllRowData = (AliL3DigitRowData*)NData;
    
    //cout << "Allocated " << endl;
    
    //Reset the data pointers, because they changed when doing UpdateRowPointer.
    for(Int_t event=0; event<NEvents; event++)
      data[event] = (AliL3DigitRowData*)hand[event]->GetDataPointer(size);
    
    tot_dig=0;
    for(Int_t row = 0 ; row < 176 ; row++){
      DigitsPerRow = 0;
      for(Int_t event = 0; event < NEvents ; event++){
	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
	for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
	  {
	    Int_t time = rowData[event][dig].fTime + sign*Offset[event];
	    if(time < transform->GetNTimeBins() && time >= 0)
	      DigitsPerRow++;
	  }
	//DigitsPerRow += data[event]->fNDigit;// + DigitsPerRow; 
      }
      //tot_dig += DigitsPerRow;
      //cout << "Found : " << DigitsPerRow << " in all Events, Row :" << row << endl;
      
      /*
	AllRowData->fRow = row;
	AllRowData->fNDigit = DigitsPerRow;
	test = (AliL3DigitData *) AllRowData->fDigitData;
      */
      
      //cout << "Got Row Pointer " << endl;
      MeVDigit = 0;
      for(Int_t event = 0; event < NEvents ; event++){
	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
	for(UInt_t digit = 0; digit < data[event]->fNDigit ; digit++){
	  
	  //Copy the digits to a temporary storage, because we have to sort them before storing them :
	  Int_t time = rowData[event][digit].fTime + sign*Offset[event];
	  if(time >= transform->GetNTimeBins() || time < 0)
	    continue;
	  
	  temp[MeVDigit] = &rowData[event][digit];
	  temp[MeVDigit]->fTime = temp[MeVDigit]->fTime + sign*Offset[event];
	  
	  //cout << MeVDigit << " : " << digit << endl;
	  MeVDigit++;
	}
	hand[event]->UpdateRowPointer(data[event]);
      }
      //Sort the digits:
      QSort(temp,0,MeVDigit);
      int final_count=0;
      
      //Now, loop and check for overlaps:
      for(Int_t c=0; c<MeVDigit; c++)
	{
	  
	  if(c>0 && temp[c]->fTime == temp[c-1]->fTime && temp[c]->fPad == temp[c-1]->fPad)
	    {//this one is overlapping with the previous one:
	      temp2[final_count-1]->fCharge += temp[c+1]->fCharge;
	      if(temp2[final_count-1]->fCharge>=1024)
		temp2[final_count-1]->fCharge=1023;//Saturation
	    }
	  else
	    {
	      temp2[final_count] = temp[c];
	      final_count++;
	    }
	}
      tot_dig+=final_count;
      AllRowData->fRow = row;
      AllRowData->fNDigit = final_count;
      test = (AliL3DigitData *) AllRowData->fDigitData;
      
      for(Int_t c=0; c<final_count; c++)//and copy them to the new event:
	{
	  test[c].fCharge = temp2[c]->fCharge;
	  test[c].fPad = temp2[c]->fPad;
	  test[c].fTime = temp2[c]->fTime;
	  if(c>0 && test[c].fTime == test[c-1].fTime && test[c].fPad == test[c-1].fPad)
	    cout<<"Overlap at row "<<row<<" pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<endl;
	  //  cout<<"Copying back, pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<" charge "<<(int)test[c].fCharge<<endl;
	}
      
      if(MeVDigit!=DigitsPerRow)
	cerr<<endl<<"Error: "<<MeVDigit<<" "<<DigitsPerRow<<endl;
      if(final_count > MeVDigit)
	cerr<<"Error; final_count "<<final_count<<" MeVDigit "<<MeVDigit<<endl;
      Byte_t *tmp = (Byte_t *)AllRowData;
      Int_t UpdateSize = sizeof(AliL3DigitRowData) + sizeof(AliL3DigitData)*final_count;//DigitsPerRow;
      tmp += UpdateSize;
      AllRowData = (AliL3DigitRowData *) tmp; 
      
    }
    if(tot_dig>DigitsTot)
      cerr<<endl<<"Mismatching digitcount "<<tot_dig<<" "<<DigitsTot<<endl;
    
    for(Int_t event = 0; event < NEvents ; event++)
      {
	data[event]=0;
	delete hand[event];
      }
    
    AliL3FileHandler *Out = new AliL3FileHandler();
    //sprintf(Carry,"/prog/alice/data/Rawdata/PileUp/digits_1_0.raw");
    AllRowData = (AliL3DigitRowData*)NData;
    sprintf(Carry,"/prog/alice/data/Rawdata/1_patch/pp/pileups/event_%d/digits_%d_0.raw",triggerevent,slice);
    Out->SetBinaryOutput(Carry);
    cout << "Open File: " << Carry << " to wite PileUp Event" << endl;
    Out->Memory2CompBinary(176,AllRowData);
    Out->CloseBinaryOutput();
    delete Out;
    delete [] NData;
  }
  
  delete transform;
}

void QSort(AliL3DigitData **a,Int_t first,Int_t last)
{
  //General sorting routine. only sorting the pointers.

  AliL3DigitData *tmp;
  int i; 
  int j;
  
  while (last - first > 1) {
    i = first;
    j = last;
    for (;;) {
      while (++i < last && CompareDigits(a[i], a[first]) < 0)
	;
      while (--j > first && CompareDigits(a[j], a[first]) > 0)
            ;
         if (i >= j)
            break;

         tmp  = a[i];
         a[i] = a[j];
         a[j] = tmp;
    }
      if (j == first) {
         ++first;
         continue;
      }
      tmp = a[first];
      a[first] = a[j];
      a[j] = tmp;
      if (j - first < last - (j + 1)) {
         QSort(a, first, j);
         first = j + 1;   // QSort(j + 1, last);
      } else {
         QSort(a, j + 1, last);
         last = j;        // QSort(first, j);
      }
   }
  
  
}

Int_t CompareDigits(AliL3DigitData *a,AliL3DigitData *b)
{

  if(a->fPad==b->fPad && a->fTime == b->fTime) return 0;
  
  if(a->fPad<b->fPad) return -1;
  if(a->fPad==b->fPad && a->fTime<b->fTime) return -1;
  
  return 1;
  
}
Commit	Line	Data
686cf001	1	#ifndef __CINT__
	2	#include "AliL3FileHandler.h"
	3	#include "AliL3DigitData.h"
	4	#include "AliL3Transform.h"
	5	#include "MakePileup.h"
	6	#include <stdio.h>
	7	#include <iostream.h>
	8
	9	void QSort(AliL3DigitData **a,Int_t first,Int_t last);
	10	Int_t CompareDigits(AliL3DigitData a,AliL3DigitData b);
	11	#endif
	12
	13	void MakePileup(Int_t triggerevent)
	14	{
	15	Int_t srow[2] = {0,175};
	16	UInt_t size;
	17	const Int_t NEvents = 25;
	18	//Int_t Offset[25] = {0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240};
	19	//Int_t Offset[25] = {0,20,40,60,80,100,120,140,160,180,200,220,240,-20,-40,-60,-80,-100,-120,-140,-160,-180,-200,-220,-240};
	20	Int_t Offset[25] = {30,60,90,120,150,180,210,240,270,300,330,360,-30,-60,-90,-120,-150,-180,-210,-240,-270,-300,-330,-360,-390};
	21	Offset[triggerevent]=0;
	22	AliL3FileHandler *hand[NEvents];
	23	AliL3DigitRowData *data[NEvents];
	24	AliL3DigitData *rowData[NEvents];
	25	AliL3DigitData *test;
	26	Char_t Carry[255];
	27	//Int_t event = 0;
	28	Int_t DigitsTot = 0;
	29	Int_t MeVDigit;
	30	Int_t DigitsPerRow = 0;
	31	//AliL3DigitRowData *NData;
	32	Byte_t *NData;
	33	Int_t tot_dig=0;
	34
	35	AliL3Transform *transform = new AliL3Transform();
	36	const Int_t maxdigits=50000;
	37	AliL3DigitData *temp[maxdigits];
	38	AliL3DigitData *temp2[maxdigits];
	39
	40
	41	for(Int_t slice=0 ; slice<36 ; slice++){
	42	for(Int_t event = 0; event < NEvents ; event++){
	43	hand[event] = new AliL3FileHandler();
	44	sprintf(Carry,"/prog/alice/data/Rawdata/1_patch/pp/event_%d/digits_%d_0.raw",event,slice);
	45	if(!hand[event]->SetBinaryInput(Carry))
	46	cerr<<" Error opening file :"<<Carry<<endl;
	47	cout << "Open File: " << Carry;
	48	hand[event]->Init(slice,0,srow);
	49	data[event] = hand[event]->CompBinary2Memory(size);
	50	cout<<" done"<<endl;
	51	hand[event]->CloseBinaryInput();
	52	}
	53	DigitsTot=0;
	54	Int_t sign = slice < 18 ? -1 : 1;
	55	for(Int_t row = 0 ; row < 176 ; row++){
	56
	57	for(Int_t event = 0; event < NEvents ; event++){
	58	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
	59	for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
	60	{
	61	Int_t time = rowData[event][dig].fTime + sign*Offset[event];
	62	if(time < transform->GetNTimeBins() && time >= 0)
	63	DigitsTot++;
	64	}
65	//DigitsTot += data[event]->fNDigit;
66	hand[event]->UpdateRowPointer(data[event]);
67	}
68
69	}
70	//cout << "Try to allocate : " << DigitsTot << " Digits in row : " << row << endl;
71
72
73	//NData = (AliL3DigitRowData*) malloc(AllDigitsSize);// Allocate !
74
75	Int_t AllDigitsSize = sizeof(AliL3DigitData) * DigitsTot + sizeof(AliL3DigitRowData) * 176;
76	NData = new Byte_t[AllDigitsSize];
77	memset(NData,0,AllDigitsSize);
78	AliL3DigitRowData AllRowData = (AliL3DigitRowData)NData;
79
80	//cout << "Allocated " << endl;
81
82	//Reset the data pointers, because they changed when doing UpdateRowPointer.
83	for(Int_t event=0; event<NEvents; event++)
84	data[event] = (AliL3DigitRowData*)hand[event]->GetDataPointer(size);
85
86	tot_dig=0;
87	for(Int_t row = 0 ; row < 176 ; row++){
88	DigitsPerRow = 0;
89	for(Int_t event = 0; event < NEvents ; event++){
90	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
91	for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
92	{
93	Int_t time = rowData[event][dig].fTime + sign*Offset[event];
94	if(time < transform->GetNTimeBins() && time >= 0)
95	DigitsPerRow++;
96	}
97	//DigitsPerRow += data[event]->fNDigit;// + DigitsPerRow;
98	}
99	//tot_dig += DigitsPerRow;
100	//cout << "Found : " << DigitsPerRow << " in all Events, Row :" << row << endl;
101
102	/*
103	AllRowData->fRow = row;
104	AllRowData->fNDigit = DigitsPerRow;
105	test = (AliL3DigitData *) AllRowData->fDigitData;
106	*/
107
108	//cout << "Got Row Pointer " << endl;
109	MeVDigit = 0;
110	for(Int_t event = 0; event < NEvents ; event++){
111	rowData[event] = (AliL3DigitData *)data[event]->fDigitData;
112	for(UInt_t digit = 0; digit < data[event]->fNDigit ; digit++){
113
114	//Copy the digits to a temporary storage, because we have to sort them before storing them :
115	Int_t time = rowData[event][digit].fTime + sign*Offset[event];
116	if(time >= transform->GetNTimeBins() \|\| time < 0)
117	continue;
118
119	temp[MeVDigit] = &rowData[event][digit];
120	temp[MeVDigit]->fTime = temp[MeVDigit]->fTime + sign*Offset[event];
121
122	//cout << MeVDigit << " : " << digit << endl;
123	MeVDigit++;
124	}
125	hand[event]->UpdateRowPointer(data[event]);
126	}
127	//Sort the digits:
128	QSort(temp,0,MeVDigit);
129	int final_count=0;
130
131	//Now, loop and check for overlaps:
132	for(Int_t c=0; c<MeVDigit; c++)
133	{
134
135	if(c>0 && temp[c]->fTime == temp[c-1]->fTime && temp[c]->fPad == temp[c-1]->fPad)
136	{//this one is overlapping with the previous one:
137	temp2[final_count-1]->fCharge += temp[c+1]->fCharge;
138	if(temp2[final_count-1]->fCharge>=1024)
139	temp2[final_count-1]->fCharge=1023;//Saturation
140	}
141	else
142	{
143	temp2[final_count] = temp[c];
144	final_count++;
145	}
146	}
147	tot_dig+=final_count;
148	AllRowData->fRow = row;
149	AllRowData->fNDigit = final_count;
150	test = (AliL3DigitData *) AllRowData->fDigitData;
151
152	for(Int_t c=0; c<final_count; c++)//and copy them to the new event:
153	{
154	test[c].fCharge = temp2[c]->fCharge;
155	test[c].fPad = temp2[c]->fPad;
156	test[c].fTime = temp2[c]->fTime;
157	if(c>0 && test[c].fTime == test[c-1].fTime && test[c].fPad == test[c-1].fPad)
158	cout<<"Overlap at row "<<row<<" pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<endl;
159	// cout<<"Copying back, pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<" charge "<<(int)test[c].fCharge<<endl;
160	}
161
162	if(MeVDigit!=DigitsPerRow)
163	cerr<<endl<<"Error: "<<MeVDigit<<" "<<DigitsPerRow<<endl;
164	if(final_count > MeVDigit)
165	cerr<<"Error; final_count "<<final_count<<" MeVDigit "<<MeVDigit<<endl;
166	Byte_t tmp = (Byte_t )AllRowData;
167	Int_t UpdateSize = sizeof(AliL3DigitRowData) + sizeof(AliL3DigitData)*final_count;//DigitsPerRow;
168	tmp += UpdateSize;
169	AllRowData = (AliL3DigitRowData *) tmp;
170
171	}
172	if(tot_dig>DigitsTot)
173	cerr<<endl<<"Mismatching digitcount "<<tot_dig<<" "<<DigitsTot<<endl;
174
175	for(Int_t event = 0; event < NEvents ; event++)
176	{
177	data[event]=0;
178	delete hand[event];
179	}
180
181	AliL3FileHandler *Out = new AliL3FileHandler();
182	//sprintf(Carry,"/prog/alice/data/Rawdata/PileUp/digits_1_0.raw");
183	AllRowData = (AliL3DigitRowData*)NData;
184	sprintf(Carry,"/prog/alice/data/Rawdata/1_patch/pp/pileups/event_%d/digits_%d_0.raw",triggerevent,slice);
185	Out->SetBinaryOutput(Carry);
186	cout << "Open File: " << Carry << " to wite PileUp Event" << endl;
187	Out->Memory2CompBinary(176,AllRowData);
188	Out->CloseBinaryOutput();
189	delete Out;
190	delete [] NData;
191	}
192
193	delete transform;
194	}
195
196	void QSort(AliL3DigitData **a,Int_t first,Int_t last)
197	{
198	//General sorting routine. only sorting the pointers.
199
200	AliL3DigitData *tmp;
201	int i;
202	int j;
203
204	while (last - first > 1) {
205	i = first;
206	j = last;
207	for (;;) {
208	while (++i < last && CompareDigits(a[i], a[first]) < 0)
209	;
210	while (--j > first && CompareDigits(a[j], a[first]) > 0)
211	;
212	if (i >= j)
213	break;
214
215	tmp = a[i];
216	a[i] = a[j];
217	a[j] = tmp;
218	}
219	if (j == first) {
220	++first;
221	continue;
222	}
223	tmp = a[first];
224	a[first] = a[j];
225	a[j] = tmp;
226	if (j - first < last - (j + 1)) {
227	QSort(a, first, j);
228	first = j + 1; // QSort(j + 1, last);
229	} else {
230	QSort(a, j + 1, last);
231	last = j; // QSort(first, j);
232	}
233	}
234
235
236	}
237
238	Int_t CompareDigits(AliL3DigitData a,AliL3DigitData b)
239	{
240
241	if(a->fPad==b->fPad && a->fTime == b->fTime) return 0;
242
243	if(a->fPad<b->fPad) return -1;
244	if(a->fPad==b->fPad && a->fTime<b->fTime) return -1;
245
246	return 1;
247
248	}