HLT/exa/MakePileup.C

   1 //$Id$
   2
   3 /**
   4   Macro for making pileup events in pp. Remember to have the MC
   5   optione switched on and that the resulting digit files will
   6   not be RLE so use the PileUp option when tracking.
   7
   8   Compile and run it by:
   9   root>.L MakePileup.C++
  10   root>MakePileup(path,triggerevent);
  11
  12   Remember to have a "digitfile" link to the root
  13   file containing the pp events in "path".
  14
  15   Authors: Roland Bramm & Anders Vestbo & Constantin Loizides
  16 */
  17
  18 #ifndef __CINT__
  19 #include "AliHLTLogger.h"
  20 #include "AliHLTFileHandler.h"
  21 #include "AliHLTDigitData.h"
  22 #include "AliHLTTransform.h"
  23 #include <TNtuple.h>
  24 #include <TRandom.h>
  25 #include <TSystem.h>
  26 #include <stdio.h>
  27 #include <iostream.h>
  28 #include <time.h>
  29
  30 void QSort(void **dPt,Int_t first,Int_t last);
  31 Int_t CompareDigits(void *dPt1,void *dPt2);
  32 #endif
  33
  34 /* Nice script I am using for the generation of pileup events:
  35
  36 --------------------------------------------------------
  37 #!/bin/bash
  38
  39 from=$1
  40 to=$2
  41
  42 dir=/tmp/pp
  43
  44 mkdir -p $dir
  45 ln -s -f /data1/AliRoot/PythiaPP_1000Events_TPConly_ev0-375_digits.root $dir/digitfile.root
  46 cd ~/l3/level3code/exa
  47
  48 for i in `seq $from $to`; do
  49  echo $i;
  50  pdir="pileup-100-$i";
  51
  52  aliroot -b runMakePileup.C\(\"$dir\",0,375,\"$pdir\",101\);
  53 done
  54 --------------------------------------------------------
  55
  56   Path should point to the path in which there is a link to the root digitsfile.
  57   Startev and Endev mark the beginning and ending of valid event numbers in that file.
  58   Pileupdir is the directory in which the binary pileuped data will be written
  59   Npiles specifies the number of pileups wanted.
  60   TriggerEvent specifies the triggerevent in the linear mode or if -1 is given the random mode.
  61   Good_tpc_ev_filename is the rootfile name where the offline good-tracks are stored (see fill.C)
  62 */
  63
  64 void MakePileup(Char_t *path,Int_t startev,Int_t endev,Char_t *pileupdir="pileup",Int_t npiles=25,Int_t triggerevent=-1,Char_t *good_tpc_ev_filename=0)
  65 {
  66   Int_t ntotalev=endev-startev+1;
  67
  68   if(ntotalev<npiles){
  69     cerr << "Total number of events in root file must be at least equal to the number of pp events to be piledup." <<endl;
  70     return;
  71   }
  72   if(npiles<1 || npiles % 2 != 1 ){
  73     cerr << "npiles must be at least one and a odd number." <<endl;
  74     return;
  75   }
  76
  77   Int_t const minTriggerEvents=5;
  78   //Int_t RefOffset[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};
  79   Int_t tshift=30*25/npiles;
  80   Int_t *Offset=new Int_t[npiles];
  81   Offset[0]=0;
  82   for(Int_t i=1;i<=(npiles-1)/2;i++)
  83     Offset[i] = i*tshift;
  84   for(Int_t i=1;i<=(npiles-1)/2;i++)
  85     Offset[i+(npiles-1)/2] = -i*tshift;
  86
  87   Int_t *EventNumbers=new Int_t[npiles];
  88   if(triggerevent<0){
  89     Int_t evdiff=endev-startev;
  90     time_t t;
  91     Int_t iseed = (abs((12*1000+1)*time(&t)))%900000000;
  92     TRandom *rand=new TRandom(iseed);
  93
  94     while(triggerevent<0){
  95       Double_t dummy=evdiff*rand->Rndm();
  96
  97       if(good_tpc_ev_filename){
  98         TFile goodfile=TFile(good_tpc_ev_filename,"READ");
  99         if(!goodfile.IsOpen()){
 100           cerr << "Could not open good tpc tracks file "<<good_tpc_ev_filename<<endl;
 101           return;
 102         }
 103         Char_t gooddummy[1000];
 104         sprintf(gooddummy,"good-pptracks-%d",(Int_t)dummy+startev);
 105         TNtuple *ngtuppel = (TNtuple*)goodfile.Get(gooddummy);
 106         if(!ngtuppel) continue;
 107         Int_t ngoodevents=(Int_t)ngtuppel->GetEntries();
 108         if(ngoodevents>=minTriggerEvents) triggerevent=(Int_t)dummy+startev;
 109         goodfile.Close();
 110       } else
 111         triggerevent=(Int_t)dummy+startev;
 112     } //found a trigger event
 113     for(Int_t i=1;i<npiles;i++){
 114       Int_t ev=triggerevent;
 115       while(ev==triggerevent){
 116         Double_t dummy=evdiff*rand->Rndm()+startev;
 117         ev=(Int_t)dummy;
 118       }
 119       EventNumbers[i]=ev;
 120     }
 121     delete rand;
 122   } else {
 123     Int_t ev=startev;
 124     for(Int_t i=1;i<npiles;i++){
 125       if(ev==triggerevent) ev++;
 126       EventNumbers[i]=ev;
 127       ev++;
 128     }
 129   }
 130   EventNumbers[0]=triggerevent;
 131   cout << "Events: ";for(Int_t i=0;i<npiles;i++) cout << EventNumbers[i] << " "; cout << endl;
 132   cout << "Offsets: ";for(Int_t i=0;i<npiles;i++) cout << Offset[i] << " ";cout << endl;
 133
 134   /////////////
 135   //exit(1);
 136   /////////////
 137
 138   UInt_t size;
 139   Int_t patch = -1; //All of slice.
 140
 141   Int_t NEvents = npiles; //for compatibility with old version
 142   AliHLTFileHandler **hand=new AliHLTFileHandler*[NEvents];
 143   AliHLTDigitRowData **data=new AliHLTDigitRowData*[NEvents];
 144   AliHLTDigitData **rowData=new AliHLTDigitData*[NEvents];
 145   AliHLTDigitData *test;
 146   Char_t Carry[1024];
 147   Int_t DigitsTot = 0;
 148   Int_t tot_dig=0;
 149   Int_t MeVDigit;
 150   Int_t DigitsPerRow = 0;
 151   Byte_t *NData;
 152
 153   const Int_t maxdigits=50000;
 154   AliHLTDigitData *temp[maxdigits];
 155   AliHLTDigitData *temp2[maxdigits];
 156
 157   //Create 1 filehander per event:
 158   sprintf(Carry,"%s/digitfile.root",path);
 159   for(Int_t event=0; event<NEvents; event++)
 160     {
 161       hand[event] = new AliHLTFileHandler();
 162       if(!hand[event]->SetAliInput(Carry))
 163         cerr<<" Error opening file :"<<Carry<<endl;
 164     }
 165
 166   //Looping over slices and the over events
 167   for(Int_t slice=0 ; slice<36 ; slice++){
 168     cout<<"\nLoading slice "<<slice<<" "<<endl;
 169     for(Int_t event = 0; event < NEvents ; event++){
 170       Int_t eventno=EventNumbers[event];
 171       hand[event]->Init(slice,patch);  //give merge option to FileHandler
 172       data[event] = hand[event]->AliAltroDigits2Memory(size,eventno,kTRUE);
 173       cout << "."<<flush;
 174     }
 175     cout<<" done"<<endl;
 176
 177     //Find out how many digits to allocate per patch
 178     DigitsTot=0;
 179     Int_t sign = slice < 18 ? -1 : 1;
 180     for(Int_t row = 0 ; row < AliHLTTransform::GetNRows(patch) ; row++){
 181       for(Int_t event = 0; event < NEvents ; event++){
 182         rowData[event] = (AliHLTDigitData *)data[event]->fDigitData;
 183         for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
 184           {
 185             Int_t time = rowData[event][dig].fTime + sign*Offset[event];
 186             if(time < AliHLTTransform::GetNTimeBins() && time >= 0)
 187               DigitsTot++;
 188           }
 189
 190         hand[event]->UpdateRowPointer(data[event]);
 191       }
 192     }
 193
 194     //cout << "Try to allocate : " << DigitsTot << endl;
 195     Int_t AllDigitsSize = sizeof(AliHLTDigitData) * DigitsTot + sizeof(AliHLTDigitRowData) * AliHLTTransform::GetNRows(patch);
 196     NData = new Byte_t[AllDigitsSize];
 197     memset(NData,0,AllDigitsSize);
 198     AliHLTDigitRowData *AllRowData = (AliHLTDigitRowData*)NData;
 199     //cout << "Allocated " << endl;
 200
 201     //Reset the data pointers, because they changed when doing UpdateRowPointer.
 202     for(Int_t event=0; event<NEvents; event++)
 203       data[event] = (AliHLTDigitRowData*)hand[event]->GetDataPointer(size);
 204
 205     //Pileup the event
 206     tot_dig=0;
 207     for(Int_t row = 0 ; row < AliHLTTransform::GetNRows(patch) ; row++){
 208       DigitsPerRow = 0;
 209       for(Int_t event = 0; event < NEvents ; event++){
 210         rowData[event] = (AliHLTDigitData *)data[event]->fDigitData;
 211         for(UInt_t dig=0; dig<data[event]->fNDigit; dig++)
 212           {
 213             Int_t time = rowData[event][dig].fTime + sign*Offset[event];
 214             if(time < AliHLTTransform::GetNTimeBins() && time >= 0)
 215               DigitsPerRow++;
 216           }
 217       }
 218       //cout << "Found : " << DigitsPerRow << " in all Events, Row :" << row << endl;
 219
 220       //Copy the digits to a temporary storage,
 221       //because we have to sort them before storing them :
 222       MeVDigit = 0;
 223       for(Int_t event = 0; event < NEvents ; event++){
 224         rowData[event] = (AliHLTDigitData *)data[event]->fDigitData;
 225         for(UInt_t digit = 0; digit < data[event]->fNDigit ; digit++){
 226
 227           Int_t time = rowData[event][digit].fTime + sign*Offset[event];
 228           if(time >= AliHLTTransform::GetNTimeBins() || time < 0)
 229             continue;
 230
 231           temp[MeVDigit] = &rowData[event][digit];
 232           temp[MeVDigit]->fTime = temp[MeVDigit]->fTime + sign*Offset[event];
 233           //cout << MeVDigit << " : " << digit << endl;
 234           MeVDigit++;
 235         }
 236         hand[event]->UpdateRowPointer(data[event]);
 237       }
 238
 239       //Sort the digits:
 240       QSort((void**)temp,0,MeVDigit);
 241
 242       //Now, loop and check for overlaps:
 243       Int_t final_count=0;
 244       for(Int_t c=0; c<MeVDigit; c++)
 245         {
 246
 247           if(final_count>0 && temp[c]->fTime == temp[c-1]->fTime && temp[c]->fPad == temp[c-1]->fPad)
 248             {
 249               //this one is overlapping with the previous one:
 250               temp2[final_count-1]->fCharge += temp[c]->fCharge;
 251               if(temp2[final_count-1]->fCharge>=1024)
 252                 temp2[final_count-1]->fCharge=1023; //Saturation
 253
 254               Int_t testev=( (temp[c]->fTrackID[0])>>22 ) & 0x3ff;
 255               if( testev == triggerevent)
 256                 {
 257                   //This digit comes from trigger event, so store the mcid
 258                   temp2[final_count-1]->fTrackID[0] = temp[c]->fTrackID[0];
 259                   temp2[final_count-1]->fTrackID[1] = temp[c]->fTrackID[1];
 260                   temp2[final_count-1]->fTrackID[2] = temp[c]->fTrackID[2];
 261                 }
 262             }
 263           else
 264             {
 265               temp2[final_count] = temp[c];
 266               final_count++;
 267             }
 268         }
 269
 270       tot_dig+=final_count;
 271       AllRowData->fRow = row;
 272       AllRowData->fNDigit = final_count;
 273       test = (AliHLTDigitData *) AllRowData->fDigitData;
 274
 275       //and copy them to the new event:
 276       for(Int_t c=0; c<final_count; c++)
 277         {
 278           test[c].fCharge = temp2[c]->fCharge;
 279           test[c].fPad = temp2[c]->fPad;
 280           test[c].fTime = temp2[c]->fTime;
 281
 282
 283           Int_t testev=( (temp2[c]->fTrackID[0])>>22 ) & 0x3ff;
 284
 285           //Save the MCID
 286           //Only store the mc corresponding to the trigger event:
 287           if ( testev == triggerevent)
 288             {
 289               test[c].fTrackID[0] = ( (temp2[c]->fTrackID[0]) & 0x3fffff ) - 128;
 290               test[c].fTrackID[1] = ( (temp2[c]->fTrackID[1]) & 0x3fffff ) - 128;
 291               test[c].fTrackID[2] = ( (temp2[c]->fTrackID[2]) & 0x3fffff ) - 128;
 292             }
 293           else //Digit not corresponding to the triggerevent, so mark it as noise.
 294             {
 295               test[c].fTrackID[0]=-1;
 296               test[c].fTrackID[1]=-1;
 297               test[c].fTrackID[2]=-1;
 298             }
 299           if(c>0 && test[c].fTime == test[c-1].fTime && test[c].fPad == test[c-1].fPad)
 300             cerr<<"Overlap at row "<<row<<" pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<endl;
 301           //cout<<"Copying back, pad "<<(int)test[c].fPad<<" time "<<(int)test[c].fTime<<" charge "<<(int)test[c].fCharge<<endl;
 302         }
 303
 304       if(MeVDigit!=DigitsPerRow)
 305         cerr<<endl<<"Error: "<<MeVDigit<<" "<<DigitsPerRow<<endl;
 306       if(final_count > MeVDigit)
 307         cerr<<"Error; final_count "<<final_count<<" MeVDigit "<<MeVDigit<<endl;
 308
 309       Byte_t *tmp = (Byte_t *)AllRowData;
 310       Int_t UpdateSize = sizeof(AliHLTDigitRowData) + sizeof(AliHLTDigitData)*final_count;
 311       tmp += UpdateSize;
 312       AllRowData = (AliHLTDigitRowData *) tmp;
 313     }//end looping over row
 314
 315     if(tot_dig>DigitsTot)
 316       cerr<<endl<<"Mismatching digitcount "<<tot_dig<<" "<<DigitsTot<<endl;
 317
 318     //cout << "Clearing data structures. " << endl;
 319     for(Int_t event = 0; event < NEvents ; event++)
 320       {
 321         data[event]=0;
 322         hand[event]->Free();
 323         hand[event]->FreeDigitsTree();
 324       }
 325
 326     AliHLTFileHandler *Out = new AliHLTFileHandler();
 327     AllRowData = (AliHLTDigitRowData*)NData;
 328     if(slice==0){
 329       if(pileupdir[0]=='/')
 330         sprintf(Carry,"mkdir -p %s/",pileupdir);
 331       else
 332         sprintf(Carry,"mkdir -p %s/%s/",path,pileupdir);
 333       gSystem->Exec(Carry);
 334     }
 335     if(pileupdir[0]=='/')
 336       sprintf(Carry,"%s/digits_%d_%d_%d.raw",pileupdir,triggerevent,slice,patch);
 337     else
 338       sprintf(Carry,"%s/%s/digits_%d_%d_%d.raw",path,pileupdir,triggerevent,slice,patch);
 339
 340
 341     if(!Out->SetBinaryOutput(Carry))
 342       {
 343         cerr<<"Cannot open file: "<<Carry<<endl;
 344         return;
 345       }
 346     cout << "Writing to file: " << Carry <<endl;
 347     Out->Memory2Binary(AliHLTTransform::GetNRows(patch),AllRowData);
 348     Out->CloseBinaryOutput();
 349     delete Out;
 350     delete [] NData;
 351   }
 352
 353   for(Int_t event = 0; event < NEvents ; event++)
 354     delete hand[event];
 355
 356   delete[] hand;
 357   delete[] rowData;
 358   delete[] data;
 359 }
 360
 361 void QSort(void **dPt,Int_t first,Int_t last)
 362 {
 363   //General sorting routine. only sorting the pointers.
 364   AliHLTDigitData **a = (AliHLTDigitData**)dPt;
 365
 366   AliHLTDigitData *tmp;
 367   int i;
 368   int j;
 369
 370   while (last - first > 1) {
 371     i = first;
 372     j = last;
 373     for (;;) {
 374       while (++i < last && CompareDigits((void*)a[i], (void*)a[first]) < 0)
 375         ;
 376       while (--j > first && CompareDigits((void*)a[j], (void*)a[first]) > 0)
 377         ;
 378       if (i >= j)
 379         break;
 380
 381       tmp  = a[i];
 382       a[i] = a[j];
 383       a[j] = tmp;
 384     }
 385     if (j == first) {
 386       ++first;
 387       continue;
 388     }
 389     tmp = a[first];
 390     a[first] = a[j];
 391     a[j] = tmp;
 392     if (j - first < last - (j + 1)) {
 393       QSort((void**)a, first, j);
 394       first = j + 1;   // QSort(j + 1, last);
 395     } else {
 396       QSort((void**)a, j + 1, last);
 397       last = j;        // QSort(first, j);
 398     }
 399   }
 400 }
 401
 402 Int_t CompareDigits(void *dPt1,void *dPt2)
 403 {
 404   AliHLTDigitData *a = (AliHLTDigitData*)dPt1;
 405   AliHLTDigitData *b = (AliHLTDigitData*)dPt2;
 406   if(a->fPad==b->fPad && a->fTime == b->fTime) return 0;
 407
 408   if(a->fPad<b->fPad) return -1;
 409   if(a->fPad==b->fPad && a->fTime<b->fTime) return -1;
 410
 411   return 1;
 412 }